1. Introducción

El agua subterránea es una fuente clave para la vida humana, los ecosistemas y la economía en muchos territorios. En México, donde existe un reconocido contraste espacio-temporal de la disponibilidad del recurso además de una concentración espacial de la población y la economía (^{Nava y Medrano-Pérez, 2019}; ^{Andrade-Velázquez y Medrano-Pérez, 2020}), el agua subterránea es considerada la fuente principal de abastecimiento y el recurso de mayor estrés en muchos territorios (^{Nava y Medrano-Pérez, 2019}; ^{Saiz-Rodríguez et al., 2019}; ^{Medrano-Pérez et al., 2022}). A diferencia del agua superficial, existe menor conocimiento sobre la dinámica del agua subterránea, lo que impacta en la gestión del recurso (^{Hatch-Kuri, 2017}; ^{López-Morales, 2017}; ^{Arreguín et al., 2018}). En este sentido, los métodos geofísicos pueden contribuir a la comprensión, gestión y protección del recurso agua en general (^{CONAGUA, 2007a}; ^{Gonzáles-Amaya et al., 2019}; ^{Saiz-Rodríguez et al., 2019}; ^{Capelo et al., 2021}).

De acuerdo con ^{Capelo et al. (2021)} las técnicas geofísicas son claves en la caracterización de cuencas y acuíferos (calidad y cantidad), en el conocimiento de las zonas de recargas e interacciones con el agua superficial, la intrusión salina, y en la comprensión de los eventos extremos y los impactos del cambio climático. Además, se consideran relevantes para el estudio de la dinámica y los procesos hidrológicos e hidrogeológicos, y en la prospección de aguas subterráneas (^{Brocca et al., 2017}; ^{McLachlan et al., 2017}; ^{McAnallen et al., 2018}; ^{Romero-Ruiz et al. 2018}; ^{Busato et al., 2019}; ^{Florsch et al., 2019}; ^{Pivetta et al., 2021}).

En México, los métodos magnéticos y eléctricos se han usado en la localización de agua subterránea en ambientes volcánicos áridos de San Luis Potosí (^{López-Loera, 2014}); en la evaluación del acuífero del Valle de Mezquital mediante métodos geológicos, geoeléctricos, petrofísicos e hidráulicos (^{Hernández-Espriú et al., 2016}); en la caracterización de acuíferos en la cuenca de México (^{Arango-Galván et al., 2016}); y en la caracterización de los lineamientos estructurales mediante la integración de datos potenciales para conocer las implicaciones hidrogeológicas en la cuenca de Villa Hidalgo, San Luis Potosí (^{Pérez-Corona et al., 2017}). Además, destacan los estudios de interpretación magnetométrica (^{Alatorre-Zamora et al., 2015}; ^{Salazar-Castro y López-Loera, 2017}; ^{Escudero-Durán et al., 2021}), y con métodos potenciales en la prospección hidrogeológica y construcción de obras hidráulicas (^{CONAGUA, 2007a}).

El uso extendido de los métodos eléctricos por encima de los métodos potenciales (magnéticos y gravimétricos) es reconocido en diferentes estudios (^{López-Loera y Tristán-González, 2013}; ^{Guevara-Mansilla et al., 2017}; ^{Escudero-Durán et al., 2021}). Sin embargo, factores como el costo, la cobertura territorial, el nivel de detalle, y la facilidad tecnológica y operativa, se consideran determinantes en la selección del método a utilizar (^{Escudero-Durán et al., 2021}). En este sentido, los datos magnéticos y gravimétricos son útiles para estudiar áreas con estructuras geológicas profundas y complejas, y la integración entre ambos contribuye al análisis de áreas extensas en múltiples temas hidrológicos (^{Pérez-Corona et al., 2017}; ^{Djebbi y Gabtni, 2018}).

Las cuencas Grijalva y Usumacinta se han caracterizado por distintas perturbaciones ambientales de origen antrópico que han impactado en la dinámica natural (^{Andrade-Velázquez y Medrano-Pérez, 2020}), y en la cantidad y calidad del agua (^{PNUD México-INECC, 2018}). A su vez, ambas cuencas son consideradas en su zona baja como una de las mayores reservas de agua subterránea (^{Toledo-Sánchez, 201l}). Sin embargo, el estudio sobre el funcionamiento del agua subterránea en esta zona es limitado, aunque se han analizado con técnicas magnéticas algunas porciones entre los estados de Chiapas y Tabasco (^{Escudero-Durán et al., 2021}).

Desde esta perspectiva, el objetivo de esta investigación es interpretar las anomalías magnéticas y gravimétricas producto del tratado de los datos, con el fin identificar y delimitar las estructuras rocosas del subsuelo con potencial acuífero de la porción oeste en la cuenca Grijalva, México. En este contexto, la principal contribución de esta investigación consiste en demostrar cómo el uso de ambos métodos geofísicos puede contribuir a la identificación de potenciales fuentes de agua en áreas extensas, favoreciendo la evaluación preliminar para estudios de mayor detalle en áreas específicas identificadas. Adicionalmente, esta investigación presenta un valor adicional al analizar una región con ecosistemas húmedos y áreas naturales protegidas interdependientes de las aguas subterráneas. Por ello, se concluye que el estudio realizado puede favorecer no sólo a la comprensión y funcionamiento de los sistemas hidrogeológicos, sino también a la toma de decisión en torno a la gestión sostenible de las aguas subterráneas y los ecosistemas naturales en el país.

2. Área de estudio

La zona de interés se encuentra entre las coordenadas de latitud norte de 18° 04’ 30.5400’’ y 16° 54’ 05.3784’’, y de longitud oeste de 92° 25’ 53.1470’’ y 93° 52’ 49.9539’’. Abarca una extensión total aproximada de 20,000 km² con elevaciones que van desde -6 hasta 2,485 metros sobre el nivel del mar (m s.n.m.). Administrativamente, se encuentra en los estados de Tabasco, Chiapas, y una pequeña porción del estado de Veracruz al este, y se asientan un total de 4,501 localidades con una población de 2,160,913 habitantes. En la porción de estudio están presentes áreas naturales como la Reserva de la Biósfera Selva El Ocote, y una porción del Parque Nacional Cañón del Sumidero y la Reserva de la Biósfera Pantanos de Centla (Figura 1a).

Fuente. Elaboración propia basado en ^{ESRI (2019)}, ^{CGIAR-CSI (2022)}, ^{SINA (2022)} e ^{INEGI (2022)}.

Figura 1 Área de estudio: a) Ubicación, hidrografía y elevación; b) Subprovincias fisiográficas; c) Subcuencas; y, d) Acuíferos. 

2.3. Entorno geológico estructural e hidrogeológico

En la región SNC-AC encontramos elementos del denominado Arco Volcánico Chiapaneco (AVC) como lo son el volcán Chichonal, las estructura volcánica inactiva Catedral y el cuerpo intrusivo Santa Fe (Figuras 1a, 1b, 2a y 2b; ^{Mora et al., 2007}: ^{Mora et al., 2012}; ^{Plascencia et al., 2014}; ^{CONAGUA, 2014b}; ^{López-Loera et al., 2020}; ^{Escudero-Duran et al., 2021}).

Fuente. Elaboración propia basado en ^{SGM (2022)}, ^{REPDA (2021)}, ^{INEGI (2022)}, ^{ESRI (2019)} y ^{Marín y Torres-Ruata (1990)}.

Figura 2 Litología (a) y Unidades hidrogeológicas (b) con las estructuras geológicas y las concesiones de agua. 

Las unidades hidrogeológicas están ligadas al origen de las cuencas del Sureste. A principios del Paleógeno el Macizo de Chiapas continuaba emergido en su parte sur, mientras que sobre las rocas carbonatadas marinas (de plataforma en Chiapas y de cuenca en Campeche) de su parte norte, se depositaban gruesas secuencias de clásticos. Durante el Oligoceno continuó el depósito de clásticos en el sureste, aunque en el área de Macuspana se empezó a desarrollar un depocentro en donde se depositaron secuencias de arcillas y sobre la Cadena Chiapas-Reforma-Akal se depositaron espesores menores de arenas y arcillas, mientras en Comalcalco comenzó la movilización de grandes volúmenes de sal (^{Padilla y Sánchez, 2007}). En el Mioceno medio ocurrió la máxima etapa de deformación que plegó y cabalgó las rocas de la cadena de Chiapas-Reforma-Akal, a partir de un décollement formando pliegues orientados NW-SE y convergencia al noreste. El gran aporte de clásticos provenientes del Macizo de Chiapas durante el Plioceno-Pleistoceno causó el depósito de varios kilómetros de espesor de sedimentos, cuya sobrecarga empezó a generar fallas orientadas NE-SW con bloque caídos al norte, formando el pilar de Akal (horst de Villahermosa), que divide a las cuencas Comalcalco y Macuspana favoreciendo el ascenso de arcillas, razón por la cual es posible encontrar este tipo de rocas dentro de las secuencias terrígenas (^{Padilla y Sánchez, 2007}; ^{Saavedra-Guerrero et al., 2019}).

La litología de la provincia fisiográfica LPT corresponde a depósitos sedimentarios con facies carbonatadas y secuencias siliciclásticas con cambios de facies del jurásico al reciente. El relleno sedimentario terrígeno que compone esta región descansa sobre grandes depresiones generadas por subsidencia y evacuación salina con presencia de fallas normales del pilar de Akal (^{Padilla y Sánchez, 2007}; ^{Saavedra-Guerrero et al., 2019}; ^{Escudero-Duran et al., 2021}). Estas características de cuenca sedimentaria que genera una interestratificación compleja de arenas, limos y arcillas forman unidades hidrogeológicas con espesores que varían de 8 a más de 200 m, y permeabilidad media a alta (^{Velázquez-Aguirre y Ordaz-Ayala, 1993}).

Respecto a la SNC-AC, esta región comprende la porción sur de la zona de estudio y está compuesta por sierras escarpadas complejas, siendo la caliza la roca predominante, con un basamento de rocas metamórficas y graníticas (^{Velázquez-Aguirre y Ordaz-Ayala, 1993}; ^{Salazar-Juárez, 2008}; ^{CEIEG, s.f.}). También es posible encontrar rocas ígneas de composición andesítica y de brecha volcánica intermedia producto de la actividad volcánica del volcán Chichonal; además de depósitos volcanoclásticos relacionados a la emisión de flujos piroclásticos del AVC (Figuras 2a y 2b; ^{CEIEG, s.f.}; ^{Mora et al., 2007}; ^{López-Loera et al., 2020}). Dentro de las porciones con menor elevación (i.e.: hacia el norte, en dirección a las LPT), se tiene una litología mayormente compuesta de rocas detríticas del tipo arenisca, sedimentos aluviales y calizas arcillosas de permeabilidad baja (Figuras 2a y 2b; ^{Marín y Torres,1990}). En la sierra de Chiapas, constituida por un macizo granítico presenta condiciones desfavorables para la explotación acuífera, reduciéndose el aprovechamiento de pequeños manantiales originados en zonas de rocas fracturadas (^{Velázquez-Aguirre y Ordaz-Ayala, 1993}). La morfología tectónica de rocas calizas está representada por fallas, anticlinales y sinclinales con orientación preferencial NW-SE, formando sistemas kársticos que permitieron la formación de fosas, pilares y cuyo flujo de agua subterránea es controlado por el patrón de fracturamiento (Figuras 2a y 2b; ^{Marín y Torres, 1990}; ^{Velázquez-Aguirre y Ordaz-Ayala, 1993}; ^{Islas et al., 1996}).

En cuanto a las concesiones de agua, a febrero 2019, de acuerdo con los datos del Registro Público de Derechos de Agua (^{REPDA, 2021}) en la zona de estudio se tienen otorgadas 1,221 concesiones de agua con un volumen total de 204.4 hectómetros cúbicos/año, concentradas en los usos agrícolas con 54.3%, publico-urbano con 20.3%, e industria con 16.5%. Estas se concentran en las LPT, la cual se caracteriza por la presencia de materiales sedimentarios (gravas, arenas y limos) y depósitos aluviales y lacustres, teniendo una permeabilidad media alta generalizada, y en donde se concentran cerca del 94% de las concesiones de la porción estudiada. En esta provincia fisiográfica se tienen norias con profundidades del nivel estático de 10 m y pozos que van de 1 a 16 m (^{CONAGUA, 2020c}, ^2020d, ^2020e, ^2020f). En la unidad que corresponde a rocas sedimentarias marinas de permeabilidad baja a alta generalizada se localizan pozos con profundidades del nivel estático de hasta 60 m (^{CONAGUA, 2020c}) (Figuras 2a y 2b).

3. Bases de Datos y Procesamiento

3.3. Procesado de base de datos

Una etapa posterior a la implementación de los filtros propuestos fue el ajuste del modelo matemático. Este ajuste se llevó a cabo a través de un semivariograma, con el cual se buscó representar de manera eficiente el comportamiento espacial de los datos, reduciendo así la incertidumbre en su distribución. Como resultado del ajuste, se obtuvieron los parámetros (i.e.: nugget, slope y anisotropía) que mejor modelaron el comportamiento y distribución de los datos. La obtención de estos parámetros es fundamental, pues con ello es posible realizar un acondicionamiento del dato, permitiendo la optimización en la interpolación de los resultados. El modelo ajustado para los datos aeromagnéticos y sus parámetros se muestran en la Figura 3a, y en el caso de los datos gravimétricos satelitales se muestran en la Figura 3b.

Fuente. Elaboración propia basado en ^{SGM (2022)}, ^{REPDA (2021)}, ^{INEGI (2022)} y ^{Marín y Torres-Ruata (1990}).

Figura 3 Litología (a) y Unidades hidrogeológicas (b) con las estructuras geológicas y las concesiones de agua. 

Realizado lo anterior, se aplicó a los datos aero-magnéticos y gravimétricos satelitales un filtro compuesto sobre las configuraciones principales de las anomalías (i.e.: Anomalía Magnética del Campo Total -AMCT- y ABC), considerados en el tratamiento de los datos es la detección y delimitación de los contactos y bordes geológicos (i.e.: Señal Analítica -SA- y TDR -Derivada Tilt-; ^{Nabighian, 1972}, ¹⁹⁸⁴; ^{Roest et al., 1992}; ^{Miller & Singh, 1994}; ^{Reyes-Moreno et al., 2019}). Este ajuste se realizó para identificar los rasgos geológicos de interés hídrico asociados al marco estructural (i.e.: fallas, fracturas, contactos y/o bordes geológicos), aplicándose filtros que permitieran la identificación y delimitación de los contactos, límites y/o bordes de las fuentes magnéticas y gravimétricas. En particular, estos filtros resaltan el ruido de alta frecuencia en ambos tipos de datos y los artefactos para el caso de los datos aeromagnéticos (i.e.: la influencia de las líneas de vuelo). Para filtrar el ruido y artefactos indeseados simultáneamente, se utilizaron de manera conjunta los filtros correspondientes al Coseno Direccional y el filtro Butterworth, según lo establecido en ^{Fedi y Florio (2003)} y ^{Mogren et al. (2009)}.

Posterior a la aplicación conjunta de los filtros, para el caso de los datos aeromagnéticos, estos se redujeron al polo magnético mediante el algoritmo propuesto por ^{Baranov y Naudy (1964)} para eliminar la asimetría en las respuestas magnéticas y asegurar que las anomalías estén simétricamente centradas sobre las fuentes que las causan. La reducción en el polo se aplicó con una inclinación y declinación de 46.41° y 3.68°, respectivamente.

Asimismo, luego del tratamiento de los datos, se continuó con la aplicación de los filtros matemáticos enfocados a delimitar las anomalías relacionadas con las estructuras geológicas de interés hídrico (i.e.: rasgos geológicos asociados al marco estructural del tipo falla, fractura, contacto y/o borde geológico) como la Señal Analítica (SA) y la Derivada Tilt (TDR, por sus siglas en inglés). En particular, si bien ambos filtros delimitan los límites de las fuentes causantes de las anomalías mediante el uso de los gradientes, las respuestas de la TDR arrojan resultados más precisos, pues a diferencia de la SA, responde tanto a fuentes someras como profundas.

Los filtros propuestos fueron aplicados sobre la Anomalía Magnética Reducida al Polo (AMRP) y la ABC. Posteriormente, se obtuvieron los resultados correspondientes a la AMCT, ABC, AMRP, Señal Analítica (^{Nabighian, 1972}, ¹⁹⁸⁴; ^{Roest et al., 1992}) y TDR (^{Miller y Singh, 1994}; ^{Reyes-Moreno et al., 2019}) utilizando el software Oasis Montaj v. 8.4 de Geosoft (^{GEOSOFT, 2019}). A su vez, se implementó la deconvolución de Euler (DE; ^{Reid et al., 1990}; ^{Thompson, 1982}) para los índices estructurales (IE) 0, 0.5, 1.0, 1.5, 2.0, 2.5 y 3.0 para el caso de los datos aeromagnéticos; y, 0, 0.5, 1.0, 1.5 y 2.0 para el caso de los datos gravimétricos satelitales, con la finalidad de construir el modelo de rocas asociadas al basamento. Las soluciones de Euler se calcularon sobre la configuración correspondiente a la AMCT y ABC. Las ventanas utilizadas para el conjunto de datos fueron de 9 x 9 datos con una tolerancia de 60% para el caso de los datos aeromagnéticos, y 9x9 datos con una tolerancia de 5% para los datos gravimétricos satelitales.

Posteriormente, se continuó con la construcción del modelo (Voxel View, [^{Golden Software, 2019}]) para ambos tipos de datos (i.e.: aeromagnéticos y gravimétricos satelitales). Para esto, las soluciones de Euler de cada uno de los índices fueron integradas en una sola base de datos para su posterior interpolación. En este sentido, ^{Escudero-Durán et al. (2021)} expresa que el uso de un índice estructural solo proyecta una parte de la estructura, así que al interpolar las soluciones conjuntas se genera una imagen que muestra la distribución completa de las fuentes magnéticas. Para realizar la interpolación de las soluciones integradas se utilizó el método de distancia inversa. Esta interpolación permitió integrar la respuesta de la fuente equivalente (i.e.: magnética y gravimétrica) asociada a la capa de roca, la cual constituye el basamento cristalino impermeable (^{Escudero-Durán et al., 2021}). Posterior a la interpolación de las soluciones integradas, se realizó un proceso de asociación de fuentes para diferenciar e identificar las respuestas de materiales asociados al basamento, respecto a las repuestas de zonas asociadas y de espesores de paquetes sedimentarios.

4. Resultados y Discusión

Los resultados obtenidos de la aplicación de los diferentes filtros se muestran a continuación. Se comienza mostrando y describiendo las configuraciones principales de las anomalías (Anomalía Magnética del Campo Total (AMCT), Anomalía Magnética Reducida al Polo (AMRP) y Anomalía de Bouger Completa (ABC); Figuras 4 y 5), seguido de la descripción de las configuraciones correspondientes a la Señal Analítica (SA) y Derivada Tilt (TDR, por sus siglas en inglés). En este sentido, las anomalías obtenidas de la TDR permitieron inferir la interpretación estructural en donde se marcan los lineamientos de acuerdo con las anomalías reflejadas y las estructuras geológicas reportadas; además al tratarse de un filtro que resalta los elementos geológicos someros (fallas, litología, fracturas, entre otros), permite inferir en cierta medida los rasgos litológicos presentes. En cambio, la SA resalta elementos geológicos profundos, destacándose: intrusivo Santa Fe, la Catedral y el complejo volcánico El Chichonal. Posteriormente, sobre la configuración de la TDR, se realiza una interpretación estructural en donde se marcan los lineamientos de acuerdo con las anomalías reflejadas y las estructuras geológicas reportadas. Asimismo, se muestran los resultados obtenidos después de aplicar la deconvolución de Euler y las visualizaciones correspondientes a los modelos elaborados mediante la asociación de las fuentes equivalentes correspondientes a las soluciones de Euler.

Figura 4 a) Configuración de anomalía magnética del campo total (AMCT); b) Configuración de la Anomalía de Completa Bouguer (ABC). 

Figura 5 a) Configuración de Anomalía Magnética del Campo Total (AMCT); b) Configuración de la Anomalía Magnética Reducida al Polo (AMRP). 

En primer lugar, la configuración de la AMCT muestra una distribución de anomalías (Figura 5), con un rango de valores del campo magnético que varía entre los -437.54 nT hasta los 496.95 nT. Estas anomalías son mayormente de carácter regional, por lo cual muestra grandes longitudes de onda con valores que van desde los ~4 km hasta los ~35 km con gradientes horizontales de hasta 120 nT/m, y con distribuciones mayormente suaves que se extienden sobre toda el área de estudio y se pueden notar a primera vista. Respecto a las repuestas con valores de tendencia positiva (74.73 nT a 496.95 nT), como se muestra en las Figuras 4a y 5a, encontramos tres de carácter regional (con longitudes de onda desde los ~24 km hasta los ~35 km) dispuestas en una anomalía parcialmente continua en forma de “gancho” de ~196 km de longitud, la cual se extiende desde la parte central-sur, continuando hacia extremo occidental y terminado en la porción central-norte de la zona de estudio (dentro de estas ubicaciones es donde se encuentran los máximos positivos). Este alto magnético presente en la porción central-norte corresponde, al cuerpo intrusivo denominado Santa Fe e influencias del distrito minero Minita, de acuerdo con lo reportado por ^{López-Loera et al. (2020)}. Esta anomalía cruza toda el área de estudio desde el sur hasta el norte. Asimismo, se aprecia un alto magnético con dirección preferencial este-oeste y un valor máximo de 148.9 nT, bifurcado y alargado de menor dimensión (~28 km en la horizontal y ~6.5 km en la vertical) cercano a la porción central de la zona de estudio, el cual se extiende de manera horizontal hacia el este. En cuanto a los valores con tendencia negativa (-26.17 nT a -437.54 nT), estos se encuentran rodeando a los altos magnéticos mencionados anteriormente. Son sobresalientes dos bajos magnéticos distribuidos de manera aislada, uno al NE con dirección noroeste-sureste (NW-SE), y el segundo en la porción central de la zona de estudio. Este segundo bajo magnético forma parte de una respuesta dipolar correspondiente al intrusivo Santa Fe, como lo reporta ^{López-Loera et al. (2020)}. Asimismo, relacionado con este mismo evento, encontramos la repuesta magnética asociada al volcán Chichonal.

De acuerdo con ^{Carrillo-de la Cruz et al. (2020)} las respuestas magnéticas presentes en las inmediaciones de zonas con actividad volcánica reciente se ven fuertemente afectadas, debido a que las anomalías termales en profundidad influyen en la magnetización de las rocas, lo que a su vez influye en las anomalías magnéticas observadas en superficie. Además, en toda el área de estudio se muestran dispersas anomalías de carácter local superpuestas a las anomalías de carácter regional mencionadas anteriormente.

En el mapa de la ABC mostrado en la Figura 4b se observan anomalías con un comportamiento un poco más regional, en comparación a la AMCT, con valores que van desde los -137.29 mGal hasta 35.56 mGal, de longitudes de onda que van desde los ~5 km hasta anomalías de ~89 km, y gradientes del orden de 115 nT con distribuciones mayormente suaves. Este mismo comportamiento lo encontramos en las anomalías de las configuraciones de SA y TDR. Los valores con tendencia positiva (-29.87 a 35.56 mGal) se observan principalmente al norte y en la porción central (alto magnético que podría estar asociado al cuerpo intrusivo Santa Fe) y central-este de la zona de estudio. En la ubicación correspondiente al volcán Chichonal se observa una anomalía circular pequeña con un valor máximo de -59.5 mGal, la cual podría estar vinculada al mismo evento geológico. Respecto a los valores con tendencia negativa (-60.50 mGal a -137.29 mGal) estos se encuentran al sur distribuidos de manera aislada por valores medios que van de los -46.23 mGal a los -58.65 mGal. Al oeste encontramos valores medios (-46.23 mGal a -58.65 mGal) con la presencia de algunos altos magnéticos distribuidos de manera dispersa.

La Figura 5b muestra la configuración de la Anomalía Magnética Reducida al Polo (AMRP), en el cual se puede observar el cambio de posición de las anomalías, mostrando de esta manera la ubicación de los cuerpos que causan estas anomalías; i.e.: el cuerpo intrusivo Santa Fe puede apreciarse de manera eficiente de acuerdo con su posición. Sobre esta configuración se calcularon los distintitos filtros propuestos.

En cuanto a las configuraciones correspondientes a la SA (Figura 6), estas ayudaron a establecer los límites de las fuentes causantes de las anomalías, pues es una herramienta que delimita los gradientes de estas (^{Nabighian, 1972}; ¹⁹⁸⁴; ^{Roest et al., 1992}). Por la característica propia del filtro, las anomalías presentes en estas configuraciones muestran gradientes muy marcados con un mayor contraste de valores. En particular, cabe mencionar que las anomalías para los datos aeromagnéticos y gravimétricos satelitales tienen diferencias al tratarse de fenómenos físicos de distinta naturaleza. En el caso analizado las anomalías magnéticas responden principalmente a materiales que se encuentran sobre el basamento, mientras que las anomalías gravimétricas al ser de origen satelital ofrecen información de materiales más profundos.

Figura 6 a) Configuración de Señal Analítica de la Anomalía Completa del Campo Total (AMCT); b) Configuración de Señal Analítica de la Anomalía de Bouguer Completa (ABC); c) Configuración de la Derivada Tilt (TDR) de la AMCT; d) Configuración de la TDR de la ABC. 

Respecto a la configuración de la SA aplicada a la AMCT (Figura 6a), en los altos magnéticos (0.01187 nT/m a 0.329254 nT/m) encontramos una gran anomalía regional de forma ovalada hacia el centro de la zona de estudio que se extiende hacia el sur con longitudes de onda del orden de los 46 km, la cual está asociada al volcán Chichonal, al cuerpo intrusivo Santa Fe e influencias del distrito minero Minita (^{López-Loera et al., 2020}). Asimismo, en la porción norte se encuentra un trend parcialmente horizontal de anomalías que coincide con la porción de terreno con los valores altimétricos más bajos, correspondiente a las LPT. Hacia el oeste encontramos otro trend, el cual se extiende hacia el sur. Esta anomalía guarda relación con el alto magnético observado en el extremo occidental de la AMCT. En cuanto a los bajos magnéticos (0.00487 nT/m a 0.00021 nT/m), las anomalías guardan una distribución similar a los observados en la AMCT y AMRP. El alto magnético central se encuentra rodeado por una serie de anomalías de valores medios y bajos (0.00785 nT/m a 0.00021 nT/m), las cuales lo separan de los trends presentes al norte y al oeste. Similar al caso de la AMCT y AMRP, encontramos anomalías, tanto positivas como negativas, de carácter local superpuestas sobre las de carácter regional.

En cuanto a la SA aplicada a la ABC, en la Figura 6b se observa una distribución de anomalías un poco más sencilla en comparación a la SA aplicada a la AMCT. Los valores con tendencia positiva (0.01366 mGal/m a 0.02351 mGal/m) los encontramos a partir de la porción central hacia el sur, abarcando parcialmente la mitad de la zona de estudio. Dentro de esta zona encontramos de manera dispersa una serie de pequeños bajos magnéticos de forma redonda, dando un aspecto “granular”. Los valores con tendencias negativas (0.00125 mGal/m a 0.00016 mGal/m) se extienden desde el centro hasta el norte de la zona de estudio. Los bajos gravimétricos (0.00125 mGal/m a 0.00016 mGal/m) máximos se encuentran aislados por una serie de anomalías con valores medios.

En las configuraciones de la TDR para el caso de la AMCT y ABC mostradas en las Figuras 6c y 6d, se observa una mayor cantidad de fuentes magnéticas de menor tamaño asociadas a fuentes superficiales. Esta herramienta, al igual que la SA, ayuda en la detección de límites o contactos entre los cuerpos que originan las anomalías (^{Miller y Singh, 1994}; ^{Reyes-Moreno et al., 2019}).

Se observan que en la configuración de la SA aplicada a la AMCT (Figura 6a) se tiene que las direcciones preferenciales de las anomalías presentes en esta configuración son NW- SE y noreste-suroeste (NE-SW). La tendencia de las respuestas de la TDR (Figura 6c) guarda cierta correlación con las repuestas de la AMCT, puesto que es posible observar el trend de altos magnéticos en forma de “gancho” presente en la AMCT, aunque con menores extensiones y gradientes más pronunciados. En la zona estudiada se tienen complejos volcánicos que de acuerdo con ^{López-Loera y Tristán-González (2013)} el análisis del magnetismo de rocas volcánicas permite la localización de zonas permeables con capacidad de almacenar agua. Por tanto, se sugieren estudios geológicos-geofísicos in situ sobre el conjunto de estructuras volcánicas, eventos geológicos, y su influencia en la dinámica hidrológica de la zona. El trend de altos magnéticos (0.725 rad a 1.557 rad) observado en la porción oeste dentro de las configuraciones de la AMCT y SA también se percibe en la TDR. Este trend sobresaliente en todas las configuraciones podría estar asociado con anticlinales reportados en la zona (^{Martínez-Méndez, 2015}). El mismo caso se tiene con el alto magnético (0.725 rad a 1.557 rad) ubicado en la poción central-norte de la zona de estudio. En la porción central se encuentra otro trend de anomalías con tendencia positiva (0.725 rad a 1.557 rad) que se extienden hacia el este de la zona de estudio. Asimismo, las anomalías del volcán Chichonal y al intrusivo Santa Fe siguen apreciándose, aunque en altos magnéticos más pequeños y aislados (0.725 rad a 1.557 rad; Figura 6c). Para el caso de la TDR de la ABC (Figura 6d) se encuentra una particular distribución de repuestas pertinente a la morfología de la configuración. En la parte que va desde la coordenada en Y con valor UTM de 1950000 m hasta el sur de la zona de estudio se encuentra una serie intercalada de anomalías aisladas, tanto de tendencia positiva (-0.650 rad a 1.566 rad) como negativa (-0.898 rad a -1.569 rad), alargadas y paralelas ente sí, con longitudes de onda del orden de 3.6 km a 10.3 km y gradientes del orden de los 0.3 rad a los 3.0 rad. Similar a la configuración de la TDR para la AMCT (Figura 6c), la dirección preferencial que guardan estas anomalías es NW-SE y NE-SW. En la parte norte de la zona de estudio se observan anomalías de longitudes de onda mayores, en donde se encuentra un alto gravimétrico de gran extensión circundado por bajos gravimétricos.

Las configuraciones correspondientes a la SA y TDR de la ABC muestran una transición de valores con tendencia negativa a valores con tendencia positiva aproximadamente en el centro de la zona de estudio (Figura 6b y 6c). Este comportamiento podría estar asociado con una zona de transición entre materiales de cobertura sedimentaria (correspondiente a la porción norte; i.e.: LPT) y a estratos con componentes de rocas volcánicas (correspondiente a la porción sur; i.e.: SNC-AC) producto de la actividad volcánica del volcán Chichonal y de los distintos elementos correspondientes al AVC (i.e.: depósitos volcanoclásticos asociados a volcanes extintos; ^{Mora et al., 2007}; ^{López-Loera et al., 2020}; ^{CEIEG, s.f.}). Este comportamiento gravimétrico coincide con la transición de los elementos fisiográficos correspondientes a las LPT (valores altimétricos menores que 100 m s.n.m.) y la SNC-AC (valores altimétricos mayores a 600 m s.n.m.).

En comparación con las configuraciones anteriores la TDR es la que mejor representa las tendencias lineales correspondientes a lineamientos estructurales, las cuales representan el mayor interés hídrico, y están distribuidas dentro de la zona de estudio para ambos casos, i.e.: AMCT y ABC (Figuras 7a y 8a). Por ello, se decidió realzar la caracterización de lineamientos sobre estas configuraciones (Figuras 7b y 8b), utilizando como referencia las respuestas de la TDR y el mapa de estructuras regionales (fallas). Esto hizo posible inferir que los lineamientos siguen las orientaciones preferenciales a las fracturas y fallas reportadas en la zona, presentándose también estructuras secundarias más locales y/o superficiales asociadas.

Figura 7 a) Marco estructural regional del área de estudio y b) lineamientos interpretados sobre la configuración correspondiente a la TDR para los datos aeromagnéticos. Nota. Se puede observar, al costado derecho de la Figura 7a, la dirección preferencial de los lineamientos con orientación NW-SE y NE-SW. 

Figura 8 a) Marco estructural regional del área de estudio y b) lineamientos interpretados sobre la configuración correspondiente a la TDR para los datos gravimétricos satelitales. Nota. Se puede observar, al costado derecho de la Figura 8a, la dirección preferencial de los lineamientos con orientación NW-SE y NNE-SSW. 

Cabe mencionar que la porción sur de la zona de estudio representa el menor interés hídrico a pesar de que las repuestas de la TDR corroboran la presencia de fallas (Figuras 7 y 8). Esto se debe a que estas respuestas coinciden con el relieve más alto correspondiente a la SNC-AC y, en comparación a la zona de la planicie, esta cuenta con un potencial acuífero mucho menor, relacionado con acuíferos fracturados (^{Velázquez-Aguirre y Ordaz-Ayala, 1993}). Lo contrario sucede a partir de la porción central de la zona estudiada, pues es aquí donde la altimetría del relieve comienza a descender. Esto propicia que el área norte de la zona estudiada cuente con el mayor potencial acuífero, puesto que se encuentra en la porción de menores elevaciones de bajo de de 100 m s.n.m., lo que ocasiona que esta zona reciba escurrimientos de la cuenca alta. Respecto a la deconvolución de Euler (DE), esta se calculó en ambos tipos de datos, y las soluciones obtenidas se muestran en la Figura 9 para el dato aeromagnético (a) y el gravimétrico satelital (b).

Figura 9 Visualización de las soluciones de la deconvolución de Euler para los IE correspondientes a ambos sets de datos: a) Soluciones para el dato aeromagnético, siendo los IE 0 (negro), 0.5 (Azul marino), 1.0 (Azul celeste), 1.5 (Verde), 2.0 (Amarillo), 2.5 (Rojo) y 3.0 (purpura); b) Soluciones para el dato gravimétrico satelital, siendo los IE 0 (negro), 0.5 (Azul marino), 1.0 (Azul celeste), 1.5 (Verde), 2.0 (Amarillo). Las zonas en blanco en el área estudiada corresponden a las porciones en donde no se tienen soluciones de Euler. En la porción inferior izquierda se muestra la vista en planta de los resultados. Nota. Los índices estructurales son adimensionales. 

En las Figuras 10a y 10b se muestran la representación (Voxel view) de los modelos correspondientes a ambos tipos de datos. Con respecto a los colores, aquellas áreas representadas con tonos azules corresponden a las zonas en donde la interpolación no encontró asociación de fuentes magnéticas equivalentes dentro de las soluciones de Euler. Por dicha razón, estas zonas pueden ser asociadas a espesores conformados por paquetes sedimentarios, los cuales representan las zonas con mayores características para formaciones acuíferas. Los espesores máximos obtenidos en ambos modelos fueron de 3,500 m. Las zonas representadas con tonos grises pueden relacionarse con estructuras asociadas al basamento de la zona de estudio. Estos resultados muestran una fuerte correlación en cuanto a su distribución en la región fisiográfica denominada SNC-AC (Figura 1b), que corresponde a la parte sur de la zona de estudio. A primera vista, puede notarse que las zonas asociadas con los espesores sedimentarios se encuentran distribuidas mayormente desde el centro de la zona de estudio hasta el norte, zona que corresponde a las LPT. Esta zona es topográficamente más baja, donde escurren los flujos provenientes de las zonas más altas (i.e.: SNC-AC). Así, los resultados observados en ambos modelos y las características topográficas de la zona muestran que el sector central-norte representa las formaciones geológicas de mayor espesor y el de mayor potencial acuífero en la zona estudiada (Figuras 11a y 11b).

Figura 10 Visualización del modelo del basamento obtenido a través de la interpolación de las soluciones de Euler para a) los datos aeromagnéticos y b) los datos gravimétricos satelitales. Las zonas de color azul corresponden a aquellas porciones de mayor espesor sedimentario asociadas con zonas de mayor potencial acuífero y los tonos grises con estructuras relacionadas con el basamento. Nota. En la porción inferior izquierda de ambas figuras se muestra la vista en planta de los resultados. 

Figura 11 Elevaciones superpuestas sobre los modelos mostrados en planta (a y b) (^{JAXA/METI ALOS PALSAR L, 2011}); y las concesiones de agua en el área de estudio (REPDA-CONAGUA) superpuestas sobre los modelos mostrados en planta (c y d). 

En adición a ello, cabe señalar que las concesiones se encuentran concentradas en el área media-baja de la zona de estudio, donde se identificaron espesores sedimentarios con una profundidad máxima de 3,500 m, de lo que se infiere que representa el límite del acuífero regional con el basamento (Figuras 11c y 11d). En estas zonas se han identificado secuencias de clásticos con espesores que varían de 2,500 m en la zona de Campeche, hasta cerca de 3,000 m en Chiapas (^{Padilla y Sánchez, 2007}). En el acuífero de Reforma se tienen aprovechamientos subterráneos con profundidad del nivel estático de 10 m en promedio para norias y los pozos están perforados a profundidades que varían entre 25 y 60 m, y caudales de 50 a 120 litros por segundo (^{CONAGUA, 2020c}). Asimismo, en los acuíferos La Chontalpa, Samaria-Cunduacán y Macuspana ubicados en la LPT se tienen pozos con profundidad media de nivel estático de entre 1 a 16 m (^{CONAGUA, 2020d}, ^2020e, ^2020f). Es decir, la identificación regional de potencial del acuífero se relaciona con unidades hidrogeológicas de buena permeabilidad.

Desde esta perspectiva, los resultados confirman la idoneidad de la presente metodología para el análisis de las formaciones geológicas de buena permeabilidad, así como la identificación preliminar de zonas y estructuras geológicas que representan direcciones preferenciales de flujo de agua subterránea a nivel cuenca. Por tanto, los resultados obtenidos se enmarcan como una herramienta preliminar de bajo coste que puede contribuir al análisis y comprensión del recurso subterráneo, además de favorecer la toma de decisión en torno al ordenamiento y manejo del territorio en esta zona, con la finalidad de proteger y garantizar eventuales o futuras fuentes de agua ante eventos de contaminación y de clima extremo.

5. Conclusión

Los resultados obtenidos al aplicar los filtros propuestos a los datos mostraron consistencia en la AMCT y la ABC, y en el caso de la configuración de la TDR esta resultó ser de mayor aporte descriptivo para ambos tipos de datos, resaltando los rasgos lineales asociados al marco estructural de la zona estudiada (i.e.: fallas, fracturas y contactos geológicos). De acuerdo con los lineamientos interceptados en la configuración TDR los rasgos resaltados en la porción central-norte representan las áreas de mayor potencial acuífero con características litológicas del tipo sedimentario, aluvial y lacustre con permeabilidad media alta generalizada. Los espesores sedimentarios observados en los modelos predominan en la parte central-norte del área de estudio con una profundidad máxima de 3,500 m. Las unidades hidrogeológicas de la parte central-norte se componen principalmente por grandes paquetes sedimentarios permeables y de alta infiltración. Los elementos identificados en las configuraciones de la TDR y los modelos (i.e.: las zonas con presencia de lineamientos estructurales (fallas, fracturas, contactos geológicos) que se relacionan con direcciones de flujos preferenciales de agua subterránea con mayor permeabilidad, así como las porciones vinculadas a espesores sedimentarios), pueden interpretarse unidades hidro-geológicas con mayor almacenamiento acuífero vinculadas a espesores sedimentarios.

En la parte centro-norte de la porción estudiada, con elevaciones topográficas menores a 100 m s.n.m., está caracterizada por materiales clásticos de mayor permeabilidad y representan el mayor potencial acuífero con espesores de hasta 3,500 m, que marcan el límite del basamento de la cuenca sedimentaria y representa una zona de almacenamiento de agua subterránea con aportes de escurrimientos superficiales provenientes de la SNC-AC. Además, en la parte centro-norte del área de estudio destaca la concentración de las concesiones de agua subterráneas en el área identificada con mayor potencial hídrico, teniendo una fuerte correlación principalmente con las zonas asociadas con espesores sedimentarios mostradas en el modelo.

La metodología de la investigación se muestra idónea para la identificación preliminar de áreas con potencial acuífero al combinar el uso de datos aeromagnéticos y gravimétricos satelitales en la identificación unidades hidrogeológicas y lineamientos estructurales que representan zonas preferenciales de flujo subterráneo. Se recomienda la aplicación de prospecciones geofísicas in situ, como sondeos y tomografías eléctricas, para la descripción y delimitación la geometría de las unidades hidrogeológicas a detalle para continuar con una mejor caracterización de los sistemas acuíferos que permite una mejor toma de decisiones relacionadas con disponibilidad de cuencas-acuíferos y ordenamiento territorial.