3.1. ENTORNO GEOLÓGICO E HIDROGEOLÓGICO

Los municipios de estudio se enmarcan dentro de la Región Hidrológica Administrativa número IV, conocida como "Balsas". Desde el punto de vista de conformación fisiográfica, la región se halla dividida por el cinturón volcánico TransMexicano, en su porción más septentrional y, en su porción sur, por la Sierra Madre del Sur, que corresponde a la mayor parte del área de estudio. Los rasgos morfológicos de la región están relacionados con esta última, específicamente con aquellos valles y sierras con altitudes medias y bajas (menores a los 2600 msnm). Los procesos geológicos determinados en ella establecen, sobre todo, la situación de los puntos considerados en esta investigación.

Desde el punto de vista geológico, ^{Esteller et al. (2019)} reconocen tres grupos principales y ofrecen una descripción más detallada de los mismos: (1) Rocas del precretácico, que involucran dos unidades: esquisto Taxco y Roca Verde Taxco Viejo. Estas rocas forman parte de un complejo metamórfico (Pt y TRtv), constituido por una alternancia de metalavas, metatobas, metareniscas y metalutitas. Esta unidad es la más vieja de la zona de estudio y se considera el basamento sobre el cual la secuencia cretácica evolucionó; (2) Rocas del Cretácico: son parte de las formaciones Xochicalco (Kix) y Morelos (Kim). La primera se compone de una serie densa de calizas, finamente laminadas, de color oscuro, que alcanza, en promedio, los 150 m. La segunda, por otra parte, consiste de calizas y dolomitas interestratificadas; y (3) Rocas del cenozoico, en las que se encuentran los siguientes afloramientos: (a) Granodiorita (Ti), su afloramiento principal se halla a 7 km al oriente de la población de Tonatico, (b) Formación Chontalcoatlán (Tpch), que consiste en una secuencia de depósitos clásticos derivados de rocas del Sistema Volcánico Zinacantépetl, (c) Travertino (Qtr), formada por un paquete de rocas evaporíticas carbonatadas que afloran en la zona de estudio y (d) Aluvión, constituida por depósitos recientes derivados del intemperismo y erosión de las unidades más antiguas. Su granulometría es variable, pudiendo conformarse de gravas y arenas, así como de arcillas en otros lugares; en general, son materiales de espesores reducidos. La litología de la zona está caracterizada por la presencia de lahares (período neógeno), meta-andesita-arenisca (período cretácico) y metavolcano sedimentario (período cretácico), aunque existen otros tipos, pero de menor presencia en la misma (Figura 1).

Desde el punto de vista hidrogeológico, la zona de estudio se enmarca en el acuífero de Tenancingo. Este acuífero, en su porción inferior, que corresponde al área de estudio, alberga materiales granulares y rocas con poca y discontinua fractura y alta compacidad, que forman acuíferos de importancia mínima. En esta región ^{Martínez-Florentino (2015)} señala la existencia de dos unidades hidrogeológicas de importancia: (1) las secuencias volcánicas del Plio-Cuaternario, formadas por andesitas y basaltos, las cuales constituyen las principales unidades hidrogeológicas con capacidad para almacenar y transmitir importantes volúmenes de agua subterránea y (2) las calizas del cretácico de la Formación Morelos, que funcionan como zonas de transmisión de agua a otros materiales. Un tercer grupo, de menor importancia, estaría determinado por materiales granulares no consolidados de, relativamente, poco espesor, que generan horizontes saturados y acuíferos locales, como en el caso del municipio de Ixtapan de la Sal. Estos materiales impiden o limitan, en buena medida, la explotación del agua a través de pozos, siendo aprovechada ésta en manantiales que presentan variabilidad de caudal en función de, principalmente, las características de precipitación. Situación contraria presenta el municipio de Tonatico, con espesores mayores que permiten la explotación del agua en pozos someros.

Para la zona se distingue un flujo subterráneo preferencial de recarga con dirección Norte-Sur, proveniente, sobre todo, del Nevado de Toluca (^{Esteller et al., 2019}). Regionalmente, dos sistemas de fallas se reconocen en la zona de estudio, la primera, con dirección NE-SO y, la segunda, NO-SE, la cual estaría determinando la circulación del agua en la región (^{Juan-Pérez, 1994}). Tales fallas determinarían, de igual forma, el flujo de aguas termales que emerge en la zona de Ixtapan de la Sal, las cuales explican las condiciones de salinidad en ellas.

4.1. CARACTERIZACIÓN HIDROGEOQUÍMICA

4.1.1. MUESTREO YMÉTODOS ANALÍTICOS

La caracterización hidrogeoquímica fue determinada, tanto espacial como temporalmente en la zona de estudio, con el análisis de dos conjuntos de datos: enero de 2015 (^{Martínez-Florentino, 2015}) y abril de 2019 (determinados en esta investigación). Cuatro tipos de fuentes fueron considerados: manantiales termales, manantiales no termales, pozos y fuentes superficiales. Dos tipos de muestras fueron considerados: (1) muestras filtradas para aniones mayores, (2) muestras filtradas y acidificadas para cationes y Arsénico. Dichas muestras se conservaron en refrigeración a 4°C.

La temperatura, el pH, la conductividad eléctrica (CE) y los sólidos disueltos totales (SDT), fueron determinados en campo con un equipo portátil Thermo Scientific Orion Star A329. La dureza total y la alcalinidad, con un equipo Digital Titrator modelo 16900, utilizando, como agentes titulantes, solución de EDTA y ácido sulfúrico, respectivamente. Los sulfatos y nitratos, en laboratorio, con metodología de colorimetría utilizando equipo portátil DR 850 y empleando como agente titulante Nitrato de Bario y Cadmio, respectivamente. Los cloruros se determinaron en laboratorio con nitrato de plata como agente de titulación. El arsénico total fue determinado por flama con generador de hidruros utilizando celda de cuarzo, con equipo Perkin Elmer AAnalyst 200, el cual tiene un límite de detección de 0.001 mg/l. El control de calidad utilizado consistió en solución elemental de 100 mg/l, Marca Environmental Express Cat#QCS-27. El contenido de dióxido de carbono (denominado EpCO₂), relativo a una muestra de agua pura, fue calculado a partir de la Ecuación 1 (^{Neal et al., 1998}):

donde: AlkμEq/l corresponde a la alcalinidad expresada en μEq/l, pH es la concentración de iones hidrógeno en el agua y pHendpoint, el valor de pH del punto final de la determinación de la concentración de la alcalinidad, el cual usualmente es 4.3, por lo que la expresión 106-pHendpoint, equivale a 31.6 μEq/l (^{Jeong, 2001}).

El error de electroneutralidad, para los iones mayores, fue calculado en las muestras como un indicador de la calidad de los análisis de hidroquímica realizados. Los errores se determinaron dentro del rango de ±10%, valor que coincide con el considerado adecuado por otros autores, tomando en cuenta los errores sistemáticos introducidos por las técnicas de titulación utilizadas en la investigación (^{Pacheco-Ávila et al., 2017}).

4.2. ANÁLISIS ESTADÍSTICO MULTIVARIADO

El conjunto de determinaciones fisicoquímicas, de los años 2015 y 2019, fue utilizado como base para el análisis de componentes principales (ACP). Fueron incluidas aquellas variables que representan los elementos mayoritarios disueltos, así como algunos parámetros fisicoquímicos de interés correlacionados con los primeros, con lo que se consideró que podría tenerse una visión más amplia de los procesos geoquímicos que tienen lugar en la zona. Las variables seleccionadas para la realización del ACP fueron: temperatura, pH, conductividad eléctrica, sólidos disueltos totales, Ca2+,k+,Mg2+,Na+,Cl-,NO3-,SO42-,HCO3-,As y EpCO2.

Previo a la realización del ACP, los resultados analíticos fueron transformados logarítmicamente y, posteriormente, normalizados (z-scores), tanto para mitigar los efectos de las escalas de medida como para asegurar que parámetros altamente variables no dominaran el análisis (^{McLeod et al., 2017}), logrando con ello que cada variable tuviera el mismo peso en el análisis estadístico. Para la normalización se utilizó la Ecuación 2:

donde: Z_i es la puntuación estándar de la muestra i, X_i, el valor de la muestra, μ, la media y S la desviación estándar.

El programa estadístico SPSS fue utilizado para llevar a cabo el ACP (^{Satyaji-Rao et al., 2010}). Se adoptaron dos criterios para determinar el número de componentes principales (CP) a seleccionar: (1) aquellos CP cuyos eigenvalores (valores propios) fueron mayores a la unidad y (2) el trazado y análisis del gráfico de sedimentación (^{Subyani y Al-Ahmadi, 2010}). Asimismo, se utilizaron dos modos de extracción: (1) el modo R, que realiza la clasificación de acuerdo con la interacción entre las variables, y (2) el modo Q, que clasifica las muestras de acuerdo a sus parámetros, es decir, establece las relaciones existentes entre ellas, permitiendo el establecimiento de grupos con similitudes o, en el sentido contrario, disimilitudes fisicoquímicas. Para facilitar la interpretación de los factores seleccionados se recurrió a la rotación ortogonal de los ejes factoriales a través del método Varimax con la normalización Kaiser (^{Subyani y Al-Ahmadi, 2010}). En este sentido, la explicación del fenómeno en estudio sigue un orden descendente, en la medida que, el primer CP, explica el máximo de la varianza de los procesos en el área de estudio, mientras que el segundo factor y los subsecuentes explican el máximo de la varianza común restante, en magnitudes (varianzas) cada vez menores.

4.3. GEOESTADÍSTICA PARA EL MAPEO DE LA CALIDAD DE AGUA

4.3.1. ANÁLISIS ESTRUCTURAL DE LOS DATOS: ANÁLISIS VARIOGRÁFICO

La geoestadística se basa en el análisis de la continuidad espacial y la variación temporal de las variables hidroquímicas, y tiene su fundamento en el cálculo del variograma experimental y teórico (análisis variográfico). Su base matemática se deriva de la Ecuación 3.

donde: n es el número de pares o puntos de muestreo separados por una distancia estándar denominada h, Z(x_i), el valor de la variable z en la localidad x_i, los cuales tienen valores medidos en la variable regionalizada Z(x) y x_i, un punto de muestreo con i=1...n (^{Webster y Oliver, 2014}).

A diferencia del enfoque adoptado por otros autores (^{Sánchez-Martos et al., 2001}), quienes utilizaron el modelo de variograma esférico sin un análisis estadístico previo de adecuación, esta investigación realizó un análisis variográfico, basado en la validación cruzada, para la selección del modelo de ajuste más adecuado para llevar a cabo la interpolación de las variables. Para ello se seleccionaron cinco índices de ajuste: el error máximo (ME), raíz cuadrada media del error (RMSE), coeficiente de determinación (CD), eficiencia del modelo (EF) y el coeficiente de masa residual (CRM). Sus expresiones matemáticas están descritas en literatura relacionada (^{Homaee et al., 2002}). Los valores esperados en un ajuste perfecto son los siguientes: ME=0, RMSE=0, CD=1, EF=1 y CRM=0.

4.3.2. INTERPOLACIÓN DE LAS VARIABLES: KRIGING ORDINARIO

Una vez realizado el análisis variográfico, es posible realizar la interpolación de las variables de interés considerando la Ecuación 4, en la que se fundamenta la técnica de kriging ordinario.

donde: z^ es el valor estimado de un atributo en un punto de interés x_o, z, el valor observado en el punto muestreado x_i, λ_i, el peso asignado al punto de muestreo y n representa el número de puntos muestreados usados para la estimación (^{Webster y Oliver, 2014}).

El kriging ordinario ha sido de los más ampliamente utilizados (^{Gorai y Kumar, 2013}). En tal sentido, ^{Ebrahimi et al. (2011)} compararon 8 métodos para interpolar la variable sodio en Kohpayeh-Segzi, Irán, para un total de 36 muestras de agua subterránea utilizada para consumo. Sus resultados muestran que, dentro de todas las metodologías de interpolación evaluadas, el kriging ordinario ofreció los mejores resultados con el menor valor del índice de ajuste RMSE.

La presente investigación, en contraposición a otros autores (^{Nas y Berktay, 2010}; ^{Gorai y Kumar, 2013}), quienes realizaron la interpolación individual de variables fisicoquímicas específicas de interés, se inclinó h acia la metodología desarrollada por ^{Sánchez-Martos et al. (2001)} y ^{McLeod et al. (2017)}, quienes optaron por la interpolación de las puntuaciones factoriales en sus zonas de estudio (resultantes de la aplicación del modo Q del ACP). El enfoque metodológico, retomado por estos últimos, posee algunas características que lo hacen particularmente atractivo para los intereses de esta investigación: (1) la superficie de predicción lograda lleva implícito el conjunto de parámetros fisicoquímicos que, en él (el componente principal), están involucrados, llevando a generalizar la influencia inherente que éstos conllevan en el agua de la zona de estudio y (2) se posibilita el análisis en un espacio ortogonal multivariado (las variables son linealmente independientes), el cual es más reducido que al realizar el análisis con las variables originalmente establecidas en la investigación.

Se utilizó el programa Jeostat para la realización del análisis variográfico, así como para la realización del kriging ordinario de las variables, y el programa Arcgis 10.;, para la generación de las superficies de predicción de las variables en función de los resultados obtenidos con el primero. Su principio de funcionamiento, al igual que el proceso de generación de los archivos derivados de su implementación, son descritos en la literatura relacionada (^{Ali-Mert y Dag, 2017}).

5.1. ANÁLISIS DE LA CARACTERIZACIÓN HIDROQUÍMICA

Las Figuras 2a y 2b muestran la composición química del agua en función del tipo de fuente muestreada: (1) manantiales termales, (2) manantiales no termales, (3) pozos y (4) fuentes superficiales. En tal sentido, para ambos períodos de muestreo, tres grupos químicos se identifican en las muestras. Para el realizado hacia el año 2019, los sitios de muestreo se distribuyen, en ellos, de la siguiente manera: tipo Ca2+-HCO3- (puntos 8, 9, 10, 13 y 14, Tabla 1), Ca²⁺-Mg²⁺-Cl^- (puntos 1, 2, 4 y 12) y Na⁺-Cl^- (puntos 3, 5, 6, 7, 11 y 15). El segundo grupo corresponde a sitios considerados no termales pero que poseen, en concordancia con los resultados obtenidos por ^{Esteller et al. (2019)}, características intermedias entre los pozos (tipo Ca2+-HCO3-)y las aguas termales (tipo Na⁺-Cl^-).

Figura 2 Diagrama de Piper mostrando la composición química de las muestras de agua para los dos períodos de muestreo establecidos en la zona de estudio: (a) 2015 y (b) 2019, (c) el diagrama que muestra la relación entre el EpCO₂ y la composición química del agua para el año 2019. 

Al revisar el parámetro SDT (Tabla 1), para los tres grupos químicos determinados (Figura 2b), en promedio, éste se cuantifica, respectivamente, en 303.5 mg/l (tipo Ca2+-HCO3-), 797.97 mg/l (tipo Ca ²⁺ -Mg ²⁺ -Cl^-) y 3214.68 mg/l (tipo Na ⁺-Cl^-), clasificando las aguas como dulces (<1000 mg/l en SDT) y salobres (>1000 mg/l en SDT). De entre estos tres tipos, las muestras de aguas termales, todas ubicadas dentro del tipo Na⁺ -Cl^-, son las que presentan las mayores concentraciones de SDT (Tabla 1).

Las concentraciones catiónicas y aniónicas predominantes en la zona están determinadas de la siguiente manera: Na⁺>Ca²⁺>K⁺>Mg²⁺, de las cuales, los primeros dos elementos representan el 58.32% y 30.94% de los cationes, respectivamente, mientras que en los aniones existe una relación descendente de HCO3->Cl->SO42->NO3-, donde los dos primeros representan el 40.27% y 40.05% de los mismos. En términos generales, las características del agua indican una naturaleza dura: 10 muestras de naturaleza muy dura (>300 mg/L CaCO₃) contra 5 de tipo duro (entre 150 a 300 mg/L CaCO₃).

Un amplio rango de las concentraciones de sulfatos caracteriza las muestras analizadas en el área de estudio (15.2-956.7 mg/l). El promedio (92.2 mg/l) de las fuentes con potencialidad de uso para consumo humano (pozos y manantiales no termales) no supera el límite máximo permisible de 250 mg/l de sulfatos (^{WHO, 2011}). Para este parámetro pueden ser clasificados tres rangos en las muestras: menores a 74.5 mg/l, entre 92.2 y 290.5 mg/l y mayores a 872.5 mg/l. Las muestras provenientes de los pozos se clasifican dentro del primer grupo. En el segundo grupo se ubican todos los manantiales no termales y dos fuentes superficiales, mientras que el rango más alto corresponde a las fuentes termales. Las concentraciones de bicarbonatos varían entre 248 y 1543.03 mg/l, mientras que el pH oscila desde 6.1 hasta 7.9. Las aguas con alto contenido de sulfatos se caracterizan por altos contenidos de bicarbonatos y, en términos generales, puede inferirse un origen común de ambos elementos por su alto nivel de correlación (R²=0.982).

Adicionalmente, y considerando que éste es un acuífero carbonatado, la disolución de dolomita y de calcita se consideran procesos importantes. En ese sentido, los sitios de muestreo con contenidos medios y altos de sulfatos presentarían mayores tasas de disolución de calcita tomando en cuenta la relación molar Mg²⁺/Ca²⁺, que, en promedio, corresponden a 0.18 y 0.23, respectivamente. Por el lado contrario, un incremento en dicha relación molar (0.51) indicaría, para las aguas con menores niveles de sulfato, una mayor disolución de dolomita, posiblemente causado por un mayor tiempo de residencia y, en consecuencia, una mayor interacción agua-roca (^{Garfias et al., 2010}). Esta disolución de dolomita es un proceso relevante en la zona de estudio, en la que existe una correlación directa entre el contenido de sulfatos y calcio (Figura 3g). En la zona se infieren procesos de dedolomitización, en los que se generan incrementos de calcio, magnesio y contenido de EpCO₂ y un descenso en la alcalinidad y el pH. Los resultados obtenidos en esta investigación coinciden con los reportados por ^{Cardenal et al. (1994)}, quienes establecen que este proceso (dedolomitización) es importante en su zona de estudio cuando el yeso o la anhidrita están presentes, situación que ha sido demostrada en los resultados aquí obtenidos (Figura 3g). En ese sentido, señalan, la reacción neta que involucra este proceso está definida por la reacción química 5, en la que "el proceso continúa con la disolución simultánea de yeso (o anhidrita) y dolomita, junto con la precipitación de la calcita".

Figura 3 Relaciones químicas entre aniones y cationes mayores para el período de muestreo del año 2019: (a) Cationes mayoritarios y cloruros, (b) (Ca2++Mg2+)-(HCO3-+SO42-) y (Na⁺+K⁺)-Cl^-, (c) Na⁺ y Cl^-, (d) (Ca²⁺+Mg²⁺) y SO42-+HCO3-, (e) (Ca²⁺+Mg²⁺) y HCO3-, (f) Ca²⁺ y HCO3-, (g) Ca²⁺ y SO42-, (h) Mg²⁺ y HCO3-, (i) Mg²⁺ y SO42-, (j) Cl^- y SO42- y (k) HCO3- y pH. 

5.2. ANÁLISIS HIDROQUÍMICO

Tomando en consideración que uno de los objetivos de esta investigación consistió en establecer el origen para algunos de los elementos químicos predominantes en el agua de la zona de estudio, se razonó necesario realizar la evaluación de las relaciones existentes entre las concentraciones de iones mayoritarios. Es por esto que la utilización de relaciones hidroquímicas que identifiquen procesos de disolución específicos se introdujo a partir de experiencias previas (^{Ma et al., 2014}), logrando con ello una aproximación hacia los principales procesos que controlan la evolución química del agua y que han sido señalados en las etapas anteriores. La Figura 3 muestra las relaciones establecidas.

5.2.1. INFLUENCIA DE LAS SALES DE CLORURO

La gráfica de iones mayoritarios fue utilizada para establecer algunas conclusiones respecto de la composición química del agua en la zona de estudio. Así se espera una relación 1:1 entre los cationes mayoritarios (Ca²⁺+Mg²⁺+Na⁺+K⁺) y el ion cloruro, estimando la disolución de las sales de cloruros (^{Ma et al., 2014}). En la Figura 3a puede apreciarse que todos los puntos se dispersan ampliamente de la equilínea, lo que indica que la disolución de sales de cloruro no es la fuente mayor de los cationes totales. Tal desviación podría estar asociada a otros procesos, como la disolución de minerales de silicatos y carbonatados, o ser consecuencia del intercambio iónico.

5.2.2. DISOLUCIÓN DE MINERALES DE SILICATO E INTERCAMBIO IÓNICO

La relación mostrada en la Figura 3b permite introducir la posible disolución de silicatos o los procesos de intercambio iónico en la zona de estudio. Si dichos procesos son significativos en la composición química del agua, la relación entre ambos parámetros utilizados Ca2++Mg2+-HCO3-+SO42-yNa++k+-Cl- debería de ser lineal y con una pendiente negativa equivalente a -1 (^{Jalali, 2005}; ^{Kim et al., 2004}). En la Figura 3b puede apreciarse que un 66% de las muestras se ubica alrededor de la equilínea mostrada, implicando que, al menos para ellas, los procesos de disolución de silicatos y el intercambio iónico podrían ser significativos en su composición química.

5.2.3. ORIGEN DEL SODIO

La relación entre sodio y cloruros ha sido utilizada para identificar los mecanismos de adquisición de salinidad. La Figura 3c, que relaciona estos elementos, muestra que el 100% de las muestras se ubica alrededor o sobre la equilínea 1:1, pues un incremento de sodio va asociado a uno de cloruros (R²= 0.9937), por lo que puede inferirse que existe un probable proceso de meteorización de las rocas (^{Meybeck, 1987}). Esta alta correlación podría indicar, adicionalmente, la disolución de halita en el medio, puesto que ésta aporta concentraciones iguales de sodio y cloruros, considerándose así que estos elementos tienen su origen en el mismo proceso.

Es conveniente, en este punto, introducir otro aspecto que vincularía el origen del sodio con la disolución de silicatos antes revisada. La relación molar Na+/Cl^- está vinculada con la liberación de sodio a partir de la disolución de minerales de silicato. En ese sentido, puede afirmarse que este proceso tiene importancia para el aporte de este elemento en la zona, puesto que un 73% de los sitios presenta una relación superior a uno. Un aporte adicional de sodio estaría relacionado con la disolución de minerales de sulfato de magnesio, esto al considerar la alta correlación existente entre sodio y sulfato (R²=0.9799), lo que podría indicar una fuente común de origen hacia el medio acuífero.

5.2.4. ORIGEN DEL CALCIO Y MAGNESIO

^{Jalali (2005)} sugiere que si los elementos calcio, magnesio, sulfatos y bicarbonatos proceden de la disolución simple de calcita, dolomita y yeso, debería conducir a un balance de carga entre estos cationes y aniones. La Figura 3d muestra la relación existente entre tales elementos. Puede apreciarse que un 53% (8 puntos) se halla alrededor y sobre la equilínea, indicando la participación de la disolución de minerales de carbonato y sulfato en la hidroquímica, mientras que 7 sitios se hallan por debajo de ella, indicando la disolución de silicatos y el intercambio iónico, determinados previamente en la Figura 3b. Desde este punto de vista, la Figura 3d determinaría que las concentraciones de calcio y magnesio, de ese 53% de las muestras, son afectadas por la disolución de minerales carbonatados, como la calcita y la dolomita. La Figura 3e, que relaciona Ca²⁺+Mg²⁺ en función de HCO3- fue utilizada para profundizar en la participación de la disolución de dolomita sobre la concentración de los primeros dos elementos. Puede apreciarse que, para un 73% de las muestras, la disolución de este mineral es de importancia pues ellas se ubican alrededor de la equilínea.

De forma complementaria, el intercambio iónico, como fue descrito anteriormente, podría estar teniendo participación en el incremento en la concentración de sodio de algunas muestras, en las que este elemento podría intercambiarse con el calcio y magnesio adsorbido en las superficies de los minerales, dando como resultado una reducción de la concentración de estos últimos y un aumento del primero. Ello se expresaría en la relación existente de concentraciones catiónicas promedio para la zona: Na⁺>Ca²⁺>K⁺>Mg²⁺ (Tabla 1).

Para complementar, la Figura 3f muestra que existe una alta correlación entre el Ca²+ y HCO3- (R²=0.9291), indicando que una fuente adicional importante de calcio en la zona podría estar conformada por la disolución de la calcita. Por otro lado, existe una alta correlación entre sulfatos y calcio (R²=0.9225), con lo que se puede inferir la existencia de procesos de disolución evaporítica. La Figura 3g es utilizada para determinar la participación de la disolución de yeso en el aporte de calcio. Puede apreciarse que sólo tres sitios (20%) se hallan alrededor de la equilínea, lo cual indica que el yeso no es una fuente importante de calcio en la zona, pudiendo así derivarse éste de la disolución de otro mineral evaporítico tal como la anhidrita.

La Figura 3h, que relaciona el magnesio en función de los bicarbonatos, muestra que pueden existir otras fuentes de magnesio en la zona, ya que el 100% de los puntos se aleja de la equilínea. Todavía cabe señalar que otra fuente adicional de magnesio podría estar constituida por la disolución de minerales de sulfato de magnesio, puesto que existe una relación relativamente alta entre éstos (R²=0.951). Tal afirmación está relacionada con el comportamiento gráfico en la Figura 3i, que relaciona magnesio en función de sulfatos. En ella un 66% de las muestras se distribuye alrededor de la equilínea mostrada, señalando un aporte sustancial del proceso antes indicado.

De manera análoga a los resultados obtenidos por ^{Garfias et al. (2010)}, puede inferirse que para los sulfatos con respecto a calcio, magnesio y bicarbonatos existe, ya sea una fuente común de origen o similares procesos responsables de sus concentraciones en las muestras analizadas, esto a partir de la alta correlación de los sulfatos con ellos (R²=0.9225, 0.951 y 0.9823, respectivamente). En definitiva, puede concluirse que los procesos de disolución de minerales de silicatos, halita, sulfato de magnesio, calcita y dolomita, así como los procesos de intercambio iónico, son los principales involucrados en la composición química del agua en la zona de estudio.

Especial mención de análisis tiene la importancia de los procesos de dedolomitización en la composición química del agua, los cuales se evidencian temporalmente en la misma. En la dedolomitización, ante un incremento de sulfatos, incrementos de calcio, magnesio y cloruros (Figuras 3g, 3i y 3j), van asociados a reducciones en el pH. Este proceso es corroborado por la relación inversa existente entre los bicarbonatos y el pH (Figura 3k), generando una subsaturación respecto a dolomita, y provocando su disolución con un consecuente incremento en las concentraciones de magnesio en el agua (^{López-Chicano et al., 2001}). De manera análoga a los resultados obtenidos por ^{Garfias et al. (2010)}, los procesos de dedolomitización incrementan los contenidos de calcio y magnesio, mientras que hacen decrecer los de alcalinidad y pH, por lo que se considera un proceso fundamental en la zona.

Además de la disolución de minerales carbonatados, el enriquecimiento con HCO3- puede originarse de la disolución de CO₂ atmosférico y del suelo (^{Jeong, 2001}), derivado de la descomposición de la materia orgánica, la cual podría ser una fuente adicional de bicarbonatos en la zona de estudio.

En ese sentido, y como puede intuirse a partir de los resultados antes mostrados, el grado de mineralización de las aguas en un acuífero kárstico depende de la conjunción de los factores volumen de recarga, tiempo de residencia, las interacciones generadas en el binomio agua-roca en el acuífero y el componente gaseoso del subsuelo (CO₂). El dióxido de carbono es generado en los suelos por la degradación de la materia orgánica y la respiración radicular, por ello el uso de suelo y la actividad biológica están directamente relacionadas con él.

De esta forma, con el objetivo de profundizar la relación existente entre el contenido de EpCO₂ de las muestras y las concentraciones de elementos de interés, la Figura 4 muestra las relaciones de las variables pH, aniones, dureza y cationes en función del EpCO₂. La relación existente entre el pH y EpCO₂ es negativa para los datos analizados (Figura 4a), situación que concuerda con los resultados obtenidos por ^{Jeong (2001)}, lo que estaría relacionado con la provisión y consumo de dióxido de carbono en el proceso de evolución del agua. En este aspecto, los valores de pH varían en el rango de 6.1 a 7.9, lo que estaría generando la presencia de, aproximadamente, 30% de HCO3- junto a un 70% de ácido carbónico (H₂ CO₃), para el límite inferior (pH=6.1), contra un 95% de HCO3- y 5% de H₂ CO₃ , para el límite superior (pH=7.9), medido en las muestras analizadas para 2019. Esto estaría indicando que, ante un pH más básico, podrían esperarse concentraciones menores de EpCO₂, situación que puede ser corroborada en la relación establecida en la Figura 4a.

Figura 4 Gráfica de los niveles de: (a) pH, (b) Sulfatos, (c) Cloruros, (d) Nitratos, (e) dureza total y (f) sodio+potasio en función del EpCO₂ para las muestras de agua para el período de muestro del año 2019. 

En la Figura 4 puede verse que altos niveles de EpCO₂ se corresponden con altos niveles de sulfatos, cloruros, dureza total, y sodio+potasio para las fuentes termales y, en mucha menor proporción, para los otros tipos de fuentes (manantiales no termales, pozos y fuentes superficiales). Los resultados difieren de los obtenidos por ^{Jeong (2001)}, quien no asigna el alto contenido de CO₂ a actividad magmática en su zona de estudio; más bien, encausa su origen a la infiltración de aguas ricas en contaminantes orgánicos. En contraposición a este autor, los resultados registrados en la presente investigación coindicen con los obtenidos por ^{Cerón et al. (1998)}, quienes determinaron un enriquecimiento en bicarbonatos junto a concentraciones altas de CO₂ en aguas termales del sureste de las Cordilleras Béticas. Adicionalmente, este factor determinaba en dicha zona, junto a procesos de dedolomitización, la formación de travertinos que pueden hallarse en ella, algo que se evidencia en algunas áreas de la presente investigación.

La investigación retomó un enfoque adicional para evaluar el efecto que tiene el EpCO₂ en la conformación de las características fisicoquímicas del agua. Para ello la Figura 2c muestra la relación entre la concentración determinada de éste en el agua y su clasificación respecto del diagrama de Piper para 2019 (Figura 2b). En este sentido, y retomando la metodología desarrollada por ^{Jeong (2001)}, pueden establecerse dos rangos principales: (1) muestras cuyo contenido de EpCO₂ oscila entre 8 y 84, que contemplan la totalidad de pozos, tres muestras superficiales y dos manantiales no termales (puntos 3-6, 8-10 y 13-14, Tabla 1) y (2) aquellas muestras con un EpCO₂ arriba de los 270, en las que se incluyen tres manantiales considerados no termales y tres procedentes de fuentes termales, que pueden considerarse contaminadas por fuentes geogénicas de la zona y que sustentarían los procesos de mezcla entre dichas fuentes al compartir esta característica común de la alta concentración de EpCO₂.

La afirmación anterior se sustenta en el hecho de lo expresado por ^{Jeong (2001)} respecto de la implementación de esta metodología, pues afirma que el vínculo establecido entre el contenido de CO₂ de las muestras analizadas y su clasificación química en el diagrama de Piper (Figura 2b), puede ser utilizado como "una herramienta analítica para identificar la contaminación orgánica e inorgánica del agua subterránea". Por ello, las altas concentraciones de EpCO₂ determinadas en las "fuentes no termales" (puntos 1, 2 y 12, Tabla 1), y que son compartidas con las fuentes termales (puntos 7, 11 y 15, Tabla 1), pueden ser adjudicadas a procesos endógenos del sistema acuífero relacionados a la mezcla de dichas fuentes.

5.3. ANÁLISIS DE COMPONENTES PRINCIPALES (ACP)

5.3.1. PROCESOS HIDROGEOQUÍMICOS

Numerosos factores geológicos están involucrados en la composición química del agua en la zona de estudio, en la que la disolución de minerales carbonatados, como la calcita y dolomita, es predominante. Adicionalmente, factores antropogénicos derivados del uso del suelo también ejercen influencia, en grados diversos, en la contaminación del agua en la zona con elementos que representan un riesgo potencial a la salud de la población. Los principales procesos involucrados en la conformación de las características químicas del agua han sido destacados previamente a través del diagrama de Piper (Figura 2) y del análisis de relaciones iónicas de importancia en la zona (Figura 3 y 4).

Las técnicas de análisis multivariado, entre las que destaca el Análisis de Componentes Principales (ACP), se formulan como una alternativa adicional a las técnicas convencionales de análisis hidroquímico (diagramas de Piper, Schöeller-Berkalof y Stiif, entre otros) que ya han sido implementadas en la zona de estudio. De igual forma, se consolidan en una elección, como lo postula ^{Hair et al. (2014)}, para la reducción de un número grande original de variables fisicoquímicas (pH, SDT, Ca²⁺, Mg²⁺, entre otros), hacia un número nuevo, menor de ellas, denominadas componentes o factores (Tabla 2). Estos factores presentan una mínima pérdida de información y permiten la inferencia de los procesos hidrogeoquímicos que tienen lugar en la zona (^{Subyani y Al-Ahmadi, 2010}; ^{Subba Rao et al., 2007}; ^{McLeod et al., 2017}).

Como un primer paso para la definición de las nuevas variables, se aplicó el criterio de Kaiser para identificar el número de componentes principales (CP) a mantener para fines de esta investigación. En función del mismo, se conservaron solamente aquellos cuyo eigenvalor fue superior a la unidad, manteniendo así, para ambos períodos de muestreo, dos CP que explicaban el 91.93% (2015) y 87.96% (2019) de la varianza en los datos, los cuales se consideran adecuados en términos de representatividad estadística de los procesos hidroquímicos estudiados. La Tabla 2 muestra las cargas de las variables para ambos períodos como resultado de la rotación Varimax con la normalización Kaiser, aplicada para facilitar la interpretación de los factores seleccionados. En función del tamaño de la muestra utilizado para esta investigación, se consideraron como representativas, para motivos de análisis, aquellas variables cuyos cargas fueron mayores o iguales a 0.7 (^{Hair et al., 2014}), las cuales se muestran en negrita para fines de visualización en la tabla antes mencionada.

Para ambos períodos, el CP-I está relacionado positivamente con aquellos elementos vinculados a los procesos de salinidad, en la que se involucran procesos geogénicos como antropogénicos de composición química del agua: sólidos disueltos totales, potasio, sodio, cloruros, nitratos, sulfatos y arsénico. Este CP explica el 84.30% y 45.01% de la varianza de los datos para 2015 y 2019, respectivamente. Es importante mencionar que existe una fuerte correlación entre los SDT con calcio, sodio, magnesio y cloruros (R²=0.9537, 0.9953, 0.9385 y 0.9975, respectivamente), implicando que estas variables químicas seguirán una tendencia similar, por lo que se esperaría que este proceso de salinidad incremente a medida que lo haga la interacción entre el agua subterránea y las formaciones geológicas. Con ello puede deducirse que la alta contribución de la variable SDT a este CP se debe a la disolución iónica en el agua.

Igualmente, la participación de los iones mayoritarios (sodio, cloruros y potasio) está asociada a la disolución mineral y a la interacción agua-roca en el acuífero (^{Jiang et al, 2015}), así como a la influencia geotermal (^{Halim et al., 2010}), aspecto, este último, de importancia en la conformación de las características químicas del agua de la zona de estudio. Este hecho ya ha sido evidenciado en el diagrama de Piper (Figura 2), y se considera de relevancia en las concentraciones de arsénico registradas en el medio acuoso, situación que puede ser corroborada al verificar la fuerte carga que este elemento tiene en el CP-I (0.882).

A partir de ello, puede inferirse que este CP-I está influenciado por factores: (1) climáticos, puesto que en él intervienen parámetros fisicoquímicos relacionados a la disolución mineral (SDT, K⁺, Na⁺, Cl^-), (2) geogénicos, relacionados a la presencia de aguas termales en la zona (As presenta una carga superior a 0.8 para ambos períodos de muestreo (Tabla 2)) y (3) antropogénicos, proceso manifestado por la presencia de nitratos en este CP. Factores todos ellos que, en su conjunto, definen las características del agua en la zona de estudio.

Esta individualización se sustenta en el principio inherente a la técnica de análisis multivariado del ACP, el cual permite la inferencia de procesos individuales, independientes los unos de los otros. Esto se confirma por lo expresado por ^{Subba Rao et al. (2007)}: "Cada componente principal puede ser usado para interpretarse como un proceso hidrogeoquímico específico o como diferentes procesos a través de un examen de las cargas (de las variables fisicoquímicas registradas en los CP)". Este CP, y sus elementos positivamente relacionados (Tabla 2), afectan, sobre todo, los puntos de muestreo 1-6, 11 y 14-15 (Tabla 1, Figura 1), considerando que éstos presentan puntuaciones factoriales mayores a cero (Figuras 5c y d), por lo que se consideran los puntos con mayor índice de disolución iónica, entre los que cabe destacar la presencia de arsénico a distintas concentraciones.

Figura 5 Distribución bidimensional de las variables en las componentes principales (a) 2015 y (b) 2019 y de los sitios de muestreo en las puntuaciones factoriales en: (c) 2015 y (d) 2019. 

Adicionalmente, el CP-II está relacionado positivamente con aquellos elementos vinculados con los procesos de alcalinidad que tienen lugar en la zona: pH, magnesio, calcio, bicarbonatos, contenido de dióxido de carbono y temperatura. Debe destacarse de la Tabla 2 que, para ambos períodos, las variables pH (CP-II) y nitratos (CP-I) presentan una carga negativa, lo que implica que aquellos sitios con puntuaciones factoriales negativas (modo Q_ del ACP) tendrán tendencia a valores más altos de ellas y, relativamente menores, de las otras variables consideradas, es decir, serán aguas más alcalinas y de mayor contenido de nitratos. Las variables con carga positiva en el CP-II afectan, sobre todo, a los puntos 7, 10-12 y 15, ya que éstos presentan puntuaciones factoriales positivas en él.

Por otra parte, profundizando en algunas variables de interés en la zona, puede mencionarse que las altas concentraciones de HCO3- en las muestras pueden tener su origen, tanto a partir de la disolución de minerales carbonatados (Figuras 2 y 3) como de la degradación de la materia orgánica, por lo que el CP-II podría ser indicativo de procesos naturales donde la materia orgánica juega un papel importante. ^{Halim et al. (2010)} sugieren que la degradación de esta materia orgánica podría tener una función en los procesos de óxido reducción, potenciando con ello la movilización de arsénico en el medio acuífero, y considerándose así éste un mecanismo posible de enriquecimiento del metaloide en la zona. Tal situación puede sustentarse en la moderada correlación existente entre el arsénico y los bicarbonatos (R²=0.6887) y la carga de este último (0.640) en el CP-I, que puede considerarse representativa, desde el punto de vista estadístico, como factor de influencia e interacción con el arsénico, implicando con ello su movilización en un ambiente anóxico en la zona de estudio.

En este sentido, ^{Gupta y Gupta (2013)}, expresan que el potencial de óxido reducción (ORP), junto al pH, pueden ser considerados los dos factores principales que controlan la especiación de arsénico y determinan su solubilidad en el agua. De esta forma, para establecer la influencia de este factor (ORP) en la zona de estudio, fue aplicada la especiación de arsénico presente utilizando el programa PHREEQC (^{Avilés et al., 2013}) para los puntos 1 y 2, cuyo valor promedio de Eh es de 82.95 mV, pudiendo así determinarse un 97.2% de especies de As (V). Estos resultados de influencia coinciden con los obtenidos por ^{Rasool et al. (2016)}, quienes determinaron una predominancia de la misma especie química de arsénico en aguas bajo condiciones oxidantes en estudios de hidrogeoquímica en Tehsil Mailsi, Punjab, Pakistan.

De igual forma, ^{Zheng et al. (2004)}, citado por ^{Halim et al. (2010)}, especifican que la mayoría de los cationes mayores liberados son incrementados a partir del CO₂ derivado de la degradación, tanto óxica como anóxica de la materia orgánica. Por esta razón, explican que el proceso general para la liberación de las altas concentraciones de HCO3- en la zona puede derivarse de las siguientes reacciones químicas:

En ese sentido, esta posible generación de HCO3- a partir de la materia orgánica, podría estar relacionada con la liberación de arsénico, pues, como lo demuestran ^{DeVore et al. (2019)}, este compuesto podría tener relación directa con el mecanismo de liberación, vía desplazamiento competitivo de iones de As(V) que pudieran estar débilmente unidos a la superficie de los sólidos del subsuelo.

De manera complementaria, a partir de las relaciones establecidas entre los elementos en el plano bidimensional dictado por los componentes principales extraídos (Figuras 5a y 5b), pueden definirse tres grupos de variables, el primero está definido por EpCO₂ y magnesio, mientras que el segundo, por bicarbonatos y calcio y el tercero está conformado por otros elementos involucrados en la composición química del agua, que incluye a los sulfatos, arsénico, sodio, potasio, cloruros, sólidos disueltos totales y conductividad eléctrica. Estas agrupaciones pueden atribuirse a la disolución de los estratos minerales existentes en el acuífero estudiado, que han sido descritas en el apartado 5.2. Análisis hidroquímico: minerales evaporíticos (anhidrita y halita), así como de calcita, dolomita y aquellos de sulfato de magnesio. La relación cercana, determinada por el ACP entre el sodio y los cloruros, confirmaría su origen común en el medio acuífero (Figura 3c).

Un hecho importante que debe de resaltarse en este punto es que, para ambos períodos de muestreo, todos los elementos relacionados al enriquecimiento iónico del agua tienen una relación inversa con el pH (ubicado en una posición opuesta en el diagrama de ACP, Figuras 5a y 5b), implicando con ello que las aguas más mineralizadas, o al menos las más enriquecidas en elementos sensibles al potencial de óxido reducción, sulfatos y arsénico, tenderán a ser menos alcalinas (^{Garfias et al., 2010}). Así pues, puede considerarse éste (pH) como uno de los principales parámetros involucrados en la conformación química del agua, lo que es particularmente importante en términos de las concentraciones de arsénico en ella, considerando que es uno de los parámetros fisicoquímicos que participan tanto en su especiación como en su solubilidad en el medio acuífero (^{Ali et al., 2019}).

Similar consideración puede establecerse respecto a los nitratos, parámetro que ocupa una posición lejana y opuesta respecto de los otros elementos analizados, infiriéndose, de esta manera, un proceso distinto de origen, específicamente antropogénico, responsable de su ingreso al medio acuoso. Por otro lado, todos los elementos relacionados a la mineralización del agua se hallan en una posición cercana al EpCO₂, con lo que puede concluirse el involucramiento de este parámetro en los procesos naturales de disolución, aspecto que confirmaría las relaciones químicas visualizadas en la Figura 4.

Al comparar las gráficas de las cargas en el plano bidimensional establecido por los dos CP extraídos para los parámetros analizados en los años 2015 y 2019 (Figuras 5a y 5b), puede determinarse que las mayores variaciones se presentan para el calcio, la conductividad eléctrica, los nitratos y el pH. Esto pondría en evidencia los procesos de disolución/concentración de elementos en función de las condiciones de flujo, la precipitación y la influencia de las actividades antropogénicas/ agrícolas desarrolladas.

Estos resultados de tendencia estacional son reforzados con los obtenidos por ^{Satyaji-Rao et al. (2010)}, quienes determinaron una variación significativa de influencia, principalmente para el calcio, en función de los cambios de las condiciones climáticas derivadas de las épocas de muestreo de su investigación. Ellos establecieron, sin lugar a dudas, una relación inherente entre dichas condiciones, expresadas en términos de la precipitación, y las condiciones hidrogeológicas, las cuales, en conjunto, determinan la mayor o menor concentración de los elementos constitutivos del agua en las zonas de estudio y, señalan, las variaciones que tienen lugar en los períodos analizados en la presente investigación (Figura 5a y 5b). De este análisis se desprende que el modo R del Análisis de componentes principales corrobora los procesos ya definidos en la etapa previa de relaciones químicas mostradas. Su aporte radica en destacar, separar e individualizar entre los procesos antropogénicos y geogénicos que determinan la composición química del agua en la zona de estudio e identifican los elementos que en ellos intervienen.

5.3.2. PUNTUACIONES FACTORIALES DE LOS COMPONENTES PRINCIPALES

Las Figuras 5c y 5d presentan la distribución de los sitios de muestreo respecto de las puntuaciones factoriales extraídas para ambos períodos que están siendo considerados (modo Q del ACP). Cuatro grupos de sitios con características similares pueden ser inferidos:

El grupo "A" reúne sitios con características termales, con altas tasas de mineralización, elevado contenido de sulfatos, CE, SDT, bicarbonatos y arsénico, manteniéndose, en el tiempo, dentro del polo correspondiente a los elementos de salinidad. Este grupo se corresponde con aguas del tipo Na⁺-Cl^- y presentan, predominantemente, puntuaciones positivas en el CP-I.

El grupo "B", para ambos períodos, coincide con los pozos, todos ubicados en el municipio de Tonatico, con aguas relativamente de baja mineralización, bajas concentraciones promedio de CE (538 μS/cm) y SDT (272 mg/l) y, en general, las menores concentraciones de cationes y aniones mayoritarios, lo cual se suma a las más bajas concentraciones de arsénico.

Este análisis confirma los resultados obtenidos en el diagrama de Piper (Figura 2), en el que se ubican dentro del grupo Ca2+-HCO3-, de aguas frescas, que han pasado por procesos de disolución de minerales carbonatados. El único parámetro químico que limitaría la aptitud de uso para consumo humano de estas fuentes estaría determinado por la concentración promedio de nitratos, que asciende a 57.44 mg/l, lo que supera el límite máximo permisible de 50 mg/l establecido desde ^{WHO (2011)}. En general, este grupo se caracteriza por tener puntuaciones negativas para ambos componentes principales.

El grupo "C", aunque presenta pequeñas variaciones, en general tiende a mantener su posición dentro del plano bidimensional establecido por las puntuaciones factoriales (Figuras 5c y 5d ). Este grupo es prioritario, pues incluye las fuentes MNT1 y MNT2 (Tabla 1), que son utilizadas para la provisión de aguas para consumo humano en el muncipio de Ixtapan de la Sal. En ellas, el arsénico se consolida como el principal elemento limitante de calidad, con una concentración promedio de 0.33 mg/l (Límite máximo permisible de 0.01 mg/l (^{WHO, 2011})). Para 2019 (Figura 5d), estas fuentes tienden a acercarse un poco más a las puntuaciones del grupo "A', infiriéndose, con ello, un incremento en los procesos de mezcla evidenciados con el diagrama de Piper (Figura 2c).

Muy posiblemente, pueden ser aguas que evolucionaron, por el proceso antes mencionado, desde el tipo Na⁺-Cl^- hacia Ca²⁺-Mg²⁺-Cl^-, por su paso por minerales carbonatados, así como por los procesos de intercambio iónico que tienen lugar en el acuífero. Su movimiento vertical en el eje del CP-II, relacionado a pH y EpCO₂, puede ser indicativo de la participación de estos parámetros fisicoquímicos en los procesos de gasificación-desgasificación y en el enriquecimiento con bicarbonatos y arsénico, principalmente. En general se caracterizan por ser aguas cuyas puntuaciones en ambos componentes principales no son mayores a 0.5.

El grupo denominado "D", que se corresponde con aguas de características Ca2+-HCO3-,Ca2+-Mg2+-Cl-yNa+-Cl-, presenta puntuaciones altas en el CP-I, relacionado a la salinidad, lo que determina que presenten concentraciones de sulfatos, arsénico, CE y SDT incluso mayores a los sitios ubicados en el grupo "C". Por el lado contrario, presentan puntuaciones negativas en el CP-II, determinando valores de pH más alcalinos.

Una conclusión adicional que puede desprenderse de las Figuras 5c y d, y que deriva de la capacidad del sistema de mantener las fuentes reunidas en grupos más o menos constantes en el tiempo con la realización del ACP, es el hecho de que éste puede considerarse un sistema capaz de homogenizar las aguas, lo que hidrogeológicamente implicaría una baja susceptibilidad de las fuentes muestreadas a la contaminación. Este aspecto hidrodinámico deberá ser evaluado a profundidad en el área de estudio en investigaciones posteriores para lograr una mejor comprensión del sistema acuífero.

La Tabla 3 muestra las concentraciones promedio de las variables involucradas en los grupos identificados en el modo Q del ACP Con ella pueden visualizarse los cambios que las variables de composición físico-química del agua van teniendo en función de la posición del grupo respecto de los ejes de puntuaciones factoriales. Puede notarse que para 2019, conforme disminuyen las puntuaciones en el eje de "Puntuaciones factoriales I" (Figura 5d), las variables CE, SDT, K⁺, Na⁺, Cl^-, Mg²⁺ y SO42- siguen una concentración en el sentido de los grupos A>D>C>B. De este modo, puede concluirse, en términos de la aplicación del modo Q del ACP, que esta etapa permitió el agrupamiento de los sitios con base en características comunes, ya que estas concentraciones revelan algunas tendencias entre los grupos identificados.

5.4. ANÁLISIS GEOESTADÍSTICO

5.4.1. DEFINICIÓN DE MODELOS DE VARIOGRAMAS

El uso de herramientas de geoestadística ha cobrado importancia en el estudio de la calidad del agua en diversas partes del mundo (^{Sánchez-Martos et al., 2001}; ^{Satyaji-Rao et al., 2010}; ^{McLeod et al., 2017}). Por tal motivo, estas técnicas fueron utilizadas para estudiar la continuidad espacial de las variables seleccionadas: las puntuaciones factoriales de CP-I y CP-II (modo Q del ACP) y las concentraciones de arsénico. La Tabla 4 resume algunos parámetros de importancia utilizados en el análisis estructural, así como los resultados de la validación cruzada, desarrollada para la selección del modelo de variograma que brindó el mejor ajuste. El modelo Gaussiano proveyó el mejor ajuste para las variables CP-I y CP-II, mientras que el modelo esférico es el que lo hace para la concentración de As.

Como una forma de establecer una comparación temporal, se aplicó el análisis variográfico a los resultados de la campaña de muestreo de enero de 2015 (Tabla no mostrada), encontrándose resultados distintos en los modelos más representativos para las variables utilizadas. Así pues, el modelo exponencial se ajustó mejor a las variables CP-I y CP-II, mientras que el Gaussiano lo hace para las concentraciones de As. Considerando que el variograma representa la variación espacial de una variable, es probable que este cambio de tendencia pueda explicarse, principalmente, por dos factores: (1) Fechas diferentes de muestreo (enero de 2015 y abril de 2019), que determinan condiciones de flujo distintas, sobre todo, en lo relacionado a las condiciones de precipitación en la zona y (2) diferencias en los puntos de muestreo establecidos para ambas investigaciones.

Esta característica de cambio del análisis variográfico es consistente con las investigaciones realizadas por otros autores, quienes reportan, de manera similar, diferencias significativas de los parámetros meseta, rango y efecto nugget, para distintas épocas de muestreo e interpolación, utilizando el variograma esférico para el mapeo de sus variables de interés (^{Sánchez-Martos et al., 2001}; ^{Satyaji-Rao et al., 2010}). Situación que lleva a concluir que existe un control regional del variograma representativo con base en factores climáticos, indiscutiblemente, así como geológicos e hidrogeológicos, determinados por el flujo subterráneo en las zonas de estudio.

Retomando la información brindada en la Tabla 4, puede verse que algunos de los índices utilizados para la validación cruzada se alejan de su valor "ideal" esperado. Esta variación entre los valores obtenidos experimentalmente y su valor ideal podrían ser ocasionados por el tamaño de muestra utilizado. ^{Webster y Oliver (2014)}, explican que, ante un aumento del tamaño de muestra, desde 87 a 434 sitios seleccionados para la medición de la variable log_w K⁺ en una región agrícola, la suma de los cuadrados de los residuos osciló de 0.01127 a 0.00368 y, la varianza explicada, lo hizo de 43.8% a 99.4%. Con ello puede concluirse que el tamaño de la muestra se consolida como un factor determinante en los resultados de la predicción lograda con el kriging ordinario para esta investigación. De este modo, ellos identifican el hecho de que, mientras mayor sea el número de datos que se tenga mayor será la exactitud lograda en los índices seleccionados para la determinación del mejor modelo de ajuste para las variables de interés.

Especial mención debe hacerse respecto a la información relacionada a la validación cruzada, específicamente con los resultados relativos al índice EF (eficiencia del modelo), cuyo valor, para el CP-II, es negativo, indicando que existe una tendencia a que los valores promedio medidos brinden una mejor estimación que los valores simulados por el modelo. Adicionalmente, el coeficiente de masa residual (CRM), negativo para el CP-I, mostraría que la tendencia de los valores simulados es a la sobreestimación, generando, posiblemente, valores en la superficie de predicción que podrían ser mayores de lo que, en campo, serían medidos en una situación real (^{Homaee et al., 2002}). Como consecuencia de las anteriores aseveraciones, puede establecerse que el aumento de los sitios de muestreo para la zona de estudio podría mejorar los resultados de la validación cruzada, así como las superficies de predicción logradas por la metodología de kriging ordinario para las variables en estudio.

5.4.2. MAPEO BIDIMENSIONAL DE LA DISTRIBUCIÓN DE LAS NUEVAS VARIABLES EN LA ZONA DE ESTUDIO

Las superficies de predicción e influencia para cada una de las variables de interés se muestran en las Figuras 6 y 7. Puede verse que la tendencia de la predicción para el CP-I, para ambos períodos, es a situar las puntuaciones más altas en el área norte de la zona de estudio, en el municipio de Ixtapan de la Sal. Este componente (CP-I) está relacionado con la salinidad en el medio y se ve influenciado por el enriquecimiento derivado de la presencia de aguas termales, situación evidenciada en los diagramas de Piper (Figuras 2b y 2c) y el análisis de relaciones químicas.

Figura 6 Mapas de distribución espacial, establecidos a partir del kriging ordinario, para las variables (a) CP-I, (b) CP-II y (c) As (mg/l) para el período de muestreo correspondiente al año 2015. 

Figura 7 Mapas de distribución espacial, establecidos a partir del kriging ordinario, para las variables (a) CP-I, (b) CP-II y (c) As (mg/l) para el período de muestreo correspondiente al año 2019. 

Para 2019 (Figura 7a), puede apreciarse que existe una tendencia a la reducción del área de influencia de las puntuaciones altas en la zona norte del área de estudio para el CP-I, aunque sí debe enfatizarse que, por el lado contrario, hay un aumento de 0.07 unidades en el límite superior establecido. Para ambos períodos el comportamiento de predicción determinado es en franjas. Similar tendencia de distribución fue determinada por ^{Yescas-Coronado et al. (2018)}, en la región de Comarca Lagunera, para la variable de carbono inorgánico en suelos que reciben enmiendas orgánicas para el desarrollo de actividades agrícolas.

La variable nitrato, que tiene carga negativa en este componente (CP-I), tiende a encontrarse hacia la zona sur de la zona de estudio, muy probablemente, como fue detallado antes en la sección de ACP y de relaciones químicas, por el desarrollo de actividades agrícolas en ella. Sus mayores concentraciones se ubican en los pozos (Tabla 1), que se hallan en el municipio de Tonatico. Con el propósito de complementar, la relación SO42-/Cl- ha sido utilizada para medir el impacto de actividades agrícolas. Así pues, 33% de las muestras presenta una proporción superior a 0.25 (puntos 8, 9, 10, 13 y 14, Tabla 1), correspondientes a los pozos y a una muestra superficial, indicando la incidencia de la lixiviación de fertilizantes en ellas (^{Satyaji-Rao et al., 2010}). Especial análisis requiere el parámetro cloruros, también asociado positivamente en el CP-I, el cual comúnmente incrementa conforme lo hace la concentración mineral. Este elemento, además de factores geogénicos, puede tener un origen desde el punto de vista antropogénico, como la actividad agrícola que ha sido demostrada en la zona de estudio.

Como elemento de interés, el CP-I está relacionado con el Arsénico (con una carga positiva alta de 0.877). Los resultados presentados en la Tabla 1 demuestran que la presencia de aguas termales en la zona norte está determinando el comportamiento y presencia de arsénico en ella, lo que fue igualmente demostrado en el diagrama de Piper (Figura 2c), es decir, de los procesos de mezcla que tienen lugar en algunas fuentes consideradas no termales, cuya composición tiene características intermedias por el proceso antes mencionado. Consecuentemente, la distribución de este elemento tiene una tendencia decreciente en la dirección noreste-suroeste.

La tendencia del CP-II (Figuras 6b y 7b), relacionado a alcalinidad, está influenciada positivamente, para el año 2019, con los elementos bicarbonatos, EpCO₂, magnesio, calcio y, negativamente, con el pH. A diferencia del comportamiento notado en el CP-I, el CP-II tiende a un comportamiento localizado en algunos puntos de la zona de estudio, llegando a las cargas menores, para ambos períodos, hacia la zona central, con el nivel más bajo de -0.37 y -0.43, para 2015 y 2019, respectivamente.

En cuanto a la predicción de la tendencia de arsénico (Figuras 6c y 7c), puede decirse que ésta es similar a la encontrada para las variables anteriores, sin embargo, debe enfatizarse un cambio significativo mayor hacia la reducción del área de influencia de las puntuaciones altas, las cuales se concentran hacia el área noreste de la zona de estudio. No obstante, a pesar de lo anterior, existe una variación en la concentración límite máxima, pasando desde 0.90 mg/l, en 2015, hacia 1.24 mg/l, en 2019, lo que significaría un aumento en su impacto en la zona. Su tendencia de distribución en franjas es coincidente con el componente principal I. ^{Sánchez-Martos et al. (2001)}, de manera coincidente con esta investigación, refuerzan la influencia de los factores climáticos en la variación de las superficies de predicción debido a la realización de muestreos de manera previa y posterior a la época de recarga de su zona de estudio.

Teniendo en cuenta la información antes revisada, puede asegurarse que el análisis realizado al interpolar las nuevas variables, los componentes principales (CP-I y CP-II), permite un análisis más completo y sintético, ya que cada una de ellas resume, en sí misma, los elementos fisicoquímicos representativos de influencia en una zona determinada. De esta forma, puede afirmarse, para la zona de estudio, que la calidad del agua decrece en el sentido suroeste-noreste, teniendo al arsénico como principal elemento limitante de calidad en función de su toxicidad. Sin embargo, otros elementos adicionales están igualmente involucrados en la pérdida de la aptitud de uso del recurso para consumo humano.