Introducción
Las condiciones fisiográficas, la distribución de la población y el desarrollo de las actividades económicas generan regiones contrastantes en cuanto a disponibilidad de recursos hídricos. En México, los climas áridos y semiáridos cubren más de la mitad del territorio. En estas zonas, el agua subterránea representa un recurso indispensable para satisfacer la demanda de agua. Sin embargo, la rápida expansión urbana e industrial, y las necesidades de agua para la agricultura han propiciado un incremento considerable en la explotación de este elemento, situación que pone en peligro el desarrollo sostenible de los acuíferos (Navarro, Garfias, & Mahlknecht, 2005).
El monitoreo permanente de los niveles del agua subterránea es esencial para un seguimiento adecuado de las variaciones en la dinámica del flujo de agua en la zona saturada del subsuelo. En muchos acuíferos de México se tienen establecidas redes de monitoreo con un largo historial de información donde se realizan mediciones manuales con sonda una o dos veces al año. Aunque la información generada es importante, se reconoce que esta frecuencia de monitoreo es limitada para entender a mayor detalle los procesos que rigen el funcionamiento del agua subterránea. Los sensores de medición automática ofrecen un método de medición alternativo que evita muchos de los problemas asociados con las mediciones manuales, como el costo de los recorridos frecuentes que deben realizarse para la toma de datos y la falta de precisión al usar sonda para medir el nivel del agua en pozos. Los sensores suministran los datos necesarios tanto para resolver problemas inmediatos como para la planificación a largo plazo (Keeland, Dowd, & Hardegree, 1997). Las primeras aplicaciones de sensores de medición en los estudios de aguas subterráneas se remontan a más de 40 años, pero la amplia disponibilidad de instrumentos comerciales ha llevado a un aumento en su uso en la última década debido a que permiten la medición rentable del nivel del agua con una alta resolución temporal (Post & Von Asmuth, 2013). En años recientes se ha buscado complementar las redes existentes con pozos instrumentados con sensores que permitan tomar los niveles del agua de forma automatizada con una frecuencia diaria.
Es deseable que la elección de los sitios a monitorear sea óptima, permitiendo obtener información representativa del acuífero a través del tiempo a bajo costo (Dhar & Datta, 2009; Dhar & Patil, 2012). Esto se logra mediante el monitoreo de un menor número de pozos o a través del monitoreo de la misma serie de pozos menos veces o por una combinación de estas dos estrategias (Herrera & Pinder, 2005). Una densidad muy alta de medición aumentará el tiempo y el costo y, por otro lado, una densidad muy baja puede resultar en la pérdida de información vital (Faisal, Shakeel, Benoit, & Jean, 2007).
Herrera (1998); Nunes, Cunha y Ribeiro (2004); Kumar, Sondhi y Phogat (2005); Faisal, Shakeel, Benoit y Jean (2007); Triki, Zairi y Dhia (2013), y Júnez-Ferreira y Herrera (2013), entre otros, han desarrollado métodos para el diseño óptimo de redes de monitoreo, en los que se busca obtener un nivel aceptable de información a través de éstas, evitando monitorear en sitios que aporten información redundante; esto, mediante el uso de una función objetivo y algún método de optimización. La naturaleza del fenómeno estudiado, así como el número de posibles diseños que pueden ser seleccionados, exigen la implementación de métodos de optimización que encuentren la red más apropiada para minimizar el error en la estimación de la(s) variable(s) de interés a monitorear. Considerando que el nivel del agua subterránea es una variable que depende no sólo del espacio sino también del tiempo es necesario el diseño de redes de monitoreo espacio-temporales (Júnez-Ferreira & Herrera, 2013).
El objetivo principal de este trabajo es proponer una metodología geoestadística para el diseño óptimo de redes de monitoreo automatizado de los niveles del agua subterránea que considere su correlación espacio-temporal. Los objetivos específicos son los siguientes: a) adaptar una metodología geoestadística espacio-temporal existente para el diseño óptimo de redes de monitoreo automatizado de los niveles del agua en un acuífero; b) probar la metodología en el acuífero Loreto del estado de Zacatecas, México.
Revisión del estado del arte para el diseño de redes de monitoreo
El diseño de redes de monitoreo ha considerado el marco hidrológico, estadístico y de modelación matemática. En el primero, la red y su programa de muestreo se definen tomando las condiciones hidrológicas del sitio; en el segundo, la red de monitoreo se define con base en inferencias estadísticas obtenidas de los datos; mientras que en el tercero se utilizan modelos numéricos de flujo del agua subterránea para determinar los sitios y frecuencias de monitoreo (Herrera & Pinder, 2005).
Marco hidrológico
Para realizar trabajos mediante este enfoque se necesita información hidrogeológica, cualitativa y cuantitativa, determinada con base en cálculos y el juicio de un hidrogeólogo. Por lo general se utiliza en programas de monitoreo del agua subterránea a largo plazo, sobre todo en la detección y cumplimiento de monitoreo de la calidad del agua subterránea y en aquellos encaminados a verificar que la calidad del agua cumpla con los reglamentos vigentes (Júnez-Ferreira, 2005). La ventaja de este enfoque es que se considera totalmente la información física de los sistemas hidrogeológicos y no se pueden despreciar los sitios clave de monitoreo. Sin embargo, la desventaja es que no existe algún criterio cuantitativo para determinar cuántos pozos se requieren (Triki, Zairi, & Dhia, 2013).
Marco estadístico
Dentro del enfoque estadístico se han realizado trabajos que consideran la correlación entre los datos para la selección de las posiciones óptimas de monitoreo en la búsqueda de objetivos específicos. Por ejemplo, Prakash y Singh (1998) utilizaron kriging para la selección de sitios óptimos, adicionales a los de la red existente, para el monitoreo del nivel del agua subterránea. En el análisis utilizaron los datos del nivel del agua de 32 pozos de observación, removieron la tendencia al semivariograma muestral calculado y ajustaron un modelo al semivariograma de los residuos obtenidos. Para la determinación de los errores utilizaron un procedimiento de validación cruzada, calcularon la varianza del error en la estimación y con base en la distribución espacial identificaron zonas potenciales para la selección de los sitios prioritarios en función de la ubicación de los valores máximos de varianza, los cuales se localizaron en las fronteras del acuífero.
Con el objetivo de optimizar recursos para el monitoreo a largo plazo, Cameron y Hunter (2000) propusieron una metodología para redes de monitoreo de la calidad del agua subterránea, con la finalidad de reducir la redundancia espacial y temporal a través de un algoritmo de optimización espacial y uno temporal. Para el análisis se utilizaron datos de dos plumas contaminantes del agua subterránea. En el algoritmo espacial se obtiene la estimación de la pluma contaminante a partir de la aplicación de kriging. El algoritmo temporal combina series de tiempo de datos de varios pozos para construir un semivariograma temporal, el cual se usa para definir frecuencias de muestreo que proporcionen datos temporales sin correlación. La redundancia temporal refleja una excesiva frecuencia en el monitoreo, mientras que la redundancia espacial indica una innecesaria densidad del número de pozos muestreados, de tal forma que los pozos que proporcionan información redundante podrían ser eliminados sin afectar la calidad de la información.
Nunes, Cunha y Ribeiro (2004) propusieron la optimización de redes de monitoreo del agua subterránea, considerando una reducción en la redundancia espacial y/o temporal, y proponen tres modelos de optimización para seleccionar el mejor subconjunto de sitios de una red de monitoreo del agua subterránea. El primero maximiza la certidumbre espacial, el segundo minimiza la redundancia temporal, y el tercero es un modelo que maximiza la certidumbre espacial y minimiza la redundancia temporal.
Faisal, Shakeel, Benoit y Jean (2007), a través de la validación cruzada, optimizaron una red piezométrica enfocada a la estimación de un balance de aguas subterráneas, y realizaron análisis geoestadísticos con los datos de nivel de agua para determinar la prioridad y/o redundancia de cada punto de medida, utilizando el método de validación cruzada. Se obtuvo una red de monitoreo óptima, la cual se usó para recalcular los componentes del balance hídrico.
Nejadkoorki, Nicholson y Hadad (2011) propusieron una metodología que consiste en la selección de las posiciones donde se obtiene la máxima varianza de la estimación de contaminantes medidos mientras es minimizada la autocorrelación espacio-temporal entre los sitios seleccionados. El método de estimación no incluye el espacio-tiempo de forma conjunta y no se determina la frecuencia de muestreo, sólo se indican las posiciones que deben ser monitoreadas.
Triki, Zairi y Dhia (2013) utilizaron geoestadística para la optimización de una red de monitoreo de la carga hidráulica en un acuífero al sureste de Túnez. El enfoque geoestadístico empleado en este trabajo se basa en el análisis de la varianza de kriging universal para la optimización de la red de monitoreo actual, combinado con la prueba de validación cruzada para identificar los pozos a ser excluidos de la red de monitoreo debido a su menor aportación en interpretaciones del nivel freático por la invariabilidad en la varianza de la estimación. Por último, añadieron pozos en zonas con alta varianza, con el fin de mejorar la cobertura espacial de la red de monitoreo.
Por otro lado, se han realizado trabajos buscando definir las frecuencias óptimas de monitoreo. Al respecto, Reza, Abrishamchi y Tajrishy (2011) desarrollaron una metodología basada en la teoría de la transinformación para la optimización de redes de monitoreo de calidad de agua subterránea en la que se determina la frecuencia de muestreo en cada sitio seleccionado. Recomiendan que las frecuencias de muestreo consideren factores como la magnitud de la concentración del contaminante, dirección del cambio, correlación con pozos vecinos, e incertidumbre de la tendencia de concentración derivada de datos históricos representativos.
Júnez-Ferreira y Herrera (2013) modifican el método propuesto por Herrera (1998) para el diseño óptimo de redes de monitoreo de la carga hidráulica empleando geoestadística espacio-temporal; en esta metodología se seleccionan posiciones y frecuencias en que debe ser monitoreada la variable de interés usando un semivariograma espacio-tiempo obtenido a partir de un análisis geoestadístico y la aplicación del filtro de Kalman. El método de optimización es de inclusiones sucesivas y la evaluación espacio-temporal de las posibles posiciones de monitoreo es implementada en tiempo real. La diferencia de la mayoría de los trabajos basados en geoestadística para el diseño de redes de monitoreo radica en que esta metodología emplea una estimación espacio-temporal.
En el presente trabajo se modifica el método propuesto por Júnez-Ferreira y Herrera (2013) para evaluar sitios disponibles para ser instrumentados con sensores de medición automatizada de los niveles del agua subterránea; para lograrlo, dentro de la optimización se considera una frecuencia de monitoreo fija en los sitios evaluados. La metodología se prueba diseñando la red de monitoreo automatizada para el acuífero Loreto en Zacatecas, México.
Marco de modelación
Loaiciga (1989) fue el primer autor que propuso un método para el diseño espacio-temporal de redes de monitoreo de la calidad del agua subterránea, y analizó en conjunto la redundancia espacial y temporal de una red de monitoreo. Utilizó el método de kriging generalizado para incluir el dominio en espacio y tiempo, y de esa forma evaluar cómo una muestra tomada de un pozo en un tiempo dado reduce la incertidumbre de todas las posiciones y tiempos en los que se obtendrían estimaciones.
Dentro del marco de modelación resulta relevante mencionar el trabajo presentado por Herrera (1998), que propuso un método para el diseño espacio-temporal de redes de monitoreo de la calidad del agua subterránea, combinando un filtro de Kalman, un modelo numérico de flujo y transporte estocástico, y un método de optimización de inclusiones sucesivas. Esta metodología permite seleccionar las posiciones y tiempos de monitoreo en forma óptima. El objetivo de la red fue minimizar la varianza total del error en la estimación en todas las posiciones y tiempos por una malla que cubre el área de estudio.
Wu (2003) aplicó el filtro de Kalman para el diseño de una red de monitoreo de aguas subterráneas; consideró la desviación estándar del error en la estimación obtenida a partir de un algoritmo de simulación que acopla el filtro de Kalman y el método de los elementos finitos. Se incluyeron pozos adicionales a los seleccionados inicialmente para el monitoreo, considerando la distribución espacial de la desviación estándar del error en la estimación y seleccionando zonas donde la desviación estándar alcanza valores máximos.
Zhang, Pinder y Herrera (2005) diseñaron una red óptima de monitoreo de la calidad del agua subterránea, combinando un filtro de Kalman y un algoritmo genético, a fin de reducir al máximo el coeficiente de variación en posiciones y tiempos establecidos. Utilizaron un método de generación de campo aleatorio, una técnica de optimización, y un simulador de flujo y transporte de las aguas subterráneas.
Briseño-Ruiz, Herrera y Júnez-Ferreira (2011) modificaron la metodología propuesta por Herrera (1998) para optimizar redes de monitoreo de los niveles del agua subterránea en el que seleccionan las posiciones de los pozos y tiempos que minimicen la incertidumbre de la carga hidráulica, empleando un modelo numérico de flujo estocástico. La metodología se dividió en dos pasos: el primero, la estimación de la carga hidráulica y la incertidumbre del error de esta estimación cuando se tienen datos de la carga hidráulica en diferentes pozos y tiempos de monitoreo; el segundo, un método que selecciona las posiciones de los pozos y tiempos de monitoreo que minimicen la predicción de la incertidumbre de la carga hidráulica obtenida del paso uno, con los cuales se definen la red de monitoreo y su programa de muestreo.
Materiales y métodos
Se modificó la metodología propuesta por Júnez-Ferreira y Herrera (2013), empleando como criterio de optimización reducir la varianza del error de estimación en ciertas posiciones espacio-temporales; la modificación consiste en evaluar la función objetivo en cada posición disponible para el monitoreo dada una frecuencia fija. De esta forma, durante la optimización se considera la correlación entre los datos generados por un sensor de medición automática de los niveles del agua subterránea.
Análisis geoestadístico espacio-temporal
Con el propósito de determinar la estructura de correlación espacio-temporal de los niveles del agua subterránea se realizó un análisis geoestadístico.
Análisis exploratorio
Este análisis se efectúa con el objetivo de caracterizar la muestra, tratando de obtener la mayor información posible a partir de los datos disponibles. Es una exploración preliminar que se lleva a cabo para identificar datos atípicos en la base de datos de la variable analizada; esto es, si las observaciones o datos obtenidos son muy diferentes unos con otros, y se verifica que la distribución es normal o al menos simétrica, que no presente tendencia y la distribución sea homogénea (Díaz, 2002; Webster & Oliver, 2007).
Análisis estructural
En este análisis se estima el semivariograma muestral espacio-temporal de los niveles del agua subterránea. Se obtiene un modelo geoestadístico para la función aleatoria que se estudia (Díaz, 2002).
Semivariograma muestral
El semivariograma es la herramienta central de la geoestadística. En geoestadística clásica, la covarianza y el variograma son funciones que describen la estructura espacial de la propiedad que está siendo estudiada. En la conceptualización espacio-temporal, esta definición se extiende para considerar también la correlación temporal. El semivariograma muestral para los incrementos Δx = (Δx, Δy) y Δt se basa en la siguiente expresión:
Donde N (Δx, Δt) es el número de pares Z (xk, ti), Z (xk + Δx, ti + Δt) separados por los incrementos Δx y Δt (Mendoza-Cáceres & Herrera, 2010).
Una vez que se cuenta con el semivariograma muestral se elige el semivariograma teórico, el cual debe reflejar la dependencia espacio-temporal que muestra el fenómeno.
Tendencia
El valor esperado de una función aleatoria puede ser constante o depender de las coordenadas de la posición. Para que una función aleatoria sea estacionaria de segundo orden o satisfaga la hipótesis intrínseca, no debe existir tendencia (Júnez-Ferreira & Herrera, 2013). Cuando el valor esperado de los datos depende de la posición espacial, el tiempo o de ambos, la función aleatoria no satisface las hipótesis de estacionariedad de segundo orden ni la intrínseca. En este caso es posible trabajar con una nueva variable llamada residuo. A continuación se menciona cómo proceder en este caso.
Una función aleatoria espacio-temporal Z(x, t) puede escribirse como:
Donde m(x, t) es una función determinista (conocida como tendencia) de coordenadas (x, t) y R(x, t) (conocida como residuo) es una función aleatoria espacio-temporal estacionaria de media cero que modela la fluctuación espacio-temporal alrededor de m(x, t). De esta manera se obtienen los datos estacionarios de media cero (residuos), los cuales son usados para calcular el semivariograma espacio-temporal (Júnez-Ferreira & Herrera, 2013).
Modelo producto-suma
El producto-suma es un modelo geoestadístico muy general de covarianza espacio-temporal y se expresa de la siguiente manera:
Para semivariogramas espacio-temporales:
Donde Cs y Ct son las funciones de covarianza espacial y temporal, respectivamente; γs y γt son sus correspondientes funciones de semivariograma; Cst (0, 0) es la varianza de γst; Cs (0) es la varianza de γs, y Ct (0) es la varianza de γt. Por definición, γst (0, 0) = γs (0) = γt (0) = 0.
La siguiente condición está implícita en la transformación covarianza-semivariograma:
Derivado de la ecuación (3), se deben cumplir las siguientes condiciones para obtener una covarianza positiva definida: k1 > 0, k2 ≥ 0 y k3 ≥ 0 (De Cesare, Myers, & Posa, 2001; De Iaco, Myers, & Posa, 2001).
Validación cruzada
El procedimiento de la validación cruzada utiliza el método denominado leave one out, el cual consiste en sacar un punto de la muestra y estimar con kriging el valor en ese punto utilizando el modelo de variograma obtenido. Por medio de la validación cruzada se midieron los errores entre los valores estimados y los medidos.
Derivación de la matriz de covarianza espacio-temporal
Una vez seleccionado el modelo de semivariograma espacio-temporal, los elementos de la matriz de covarianza Cst (hs, ht) a priori requerida por el método se calculan usando:
Este modelo permite estimar valores en ubicaciones en las que no se tienen datos, así como la predicción en tiempos futuros, por lo que resulta de gran utilidad para datos de agua subterránea (De Iaco, Myers, & Posa, 2003).
Filtro de Kalman estático
El filtro de Kalman es un conjunto de ecuaciones matemáticas que provee una estimación lineal insesgada de varianza mínima para el estado de un sistema, utilizando datos con ruido (Jazwinski, 1970). La ecuación lineal que relaciona el vector estado de la variable c en las posiciones y tiempos de interés, con datos medidos en z es:
Donde zj, j = 1, 2,… es una secuencia de mediciones del nivel del agua; Hj es la matriz de muestreo de 1xN, que es distinta de cero sólo en las posiciones de espacio-tiempo correspondientes a la entrada de c, donde se toma la muestra j-ésima y N es la dimensión del vector c; c = {cip} es el vector espacio-tiempo con los valores de los niveles del agua subterránea en las posiciones y tiempos de interés. Los vectores {vj j = 1, 2,…} son los errores de medición (Júnez-Ferreira & Herrera, 2013). La matriz de covarianza del error de estimación es:
Donde ĉn = E{(c/z1, z2,…zn)} es el valor esperado de c, dadas las mediciones z1, z2,…zn. El superíndice identifica el número n de vectores de medición usados para obtener una estimación.
Para implementar el filtro, se requiere proponer una estimación a priori de c (llamada ĉ0) y de su matriz de covarianza espacio-temporal. De estas estimaciones a priori, la estimación lineal de varianza mínima para c puede obtenerse secuencialmente mediante las siguientes fórmulas:
La matriz de covarianza del error de estimación a priori (P0) se obtiene de un análisis geoestadístico mediante el ajuste de un modelo de semivariograma espacio-temporal al semivariograma muestral de c. Una vez seleccionado el modelo, los elementos de la matriz de covarianza se calculan con la ecuación (6). Para la elección de los puntos y tiempos de monitoreo se utiliza un método de optimización en el que se minimiza la varianza del error de estimación; esta función es la suma de la varianza del error en la estimación sobre todas las posiciones y tiempos, a la que llamamos la varianza total del error de estimación y se calcula:
Donde
Selección de las posiciones óptimas de monitoreo
En la formulación de Júnez-Ferreira & Herrera (2013), el método evalúa las posiciones espacio-temporales empleando las ecuaciones del filtro de Kalman estático. En la aplicación del filtro de Kalman se estima c en una malla conformada por la unión de dos conjuntos de nodos definidos en el método de optimización; el primer conjunto está integrado por las posibles posiciones de monitoreo:
Donde M es el conjunto de posibles posiciones de monitoreo espacio-temporales, Npm es el número de posiciones de monitoreo y Ntm es el número de tiempos de monitoreo, el segundo por los puntos de estimación:
Donde E es el conjunto de puntos espacio-temporales de estimación; Npe, el número de posiciones de estimación, y Nte es el número de tiempos de estimación. A través del filtro de Kalman se asigna un orden de prioridad a cada uno de los elementos de M, en función de la forma en que minimizan la varianza del error en la estimación en los puntos de estimación espacio-temporales.
En la metodología presentada en este trabajo se evalúa cada posición de monitoreo (xi), considerando que será monitoreada para todos los tiempos p = 1, 2,…Ntm, pues se pretende instalar sensores que recolectarán de manera automática la información para todos los tiempos. Esto es, para evaluar un pozo ubicado en la posición xi, la matriz de covarianza será actualizada por el filtro de Kalman, empleando las posiciones espacio-temporales {(xi, tp), p = 1,…,Ntm}.
Para seleccionar las posiciones de monitoreo se utiliza un algoritmo iterativo que escoge una posición en el espacio a la vez. Cada posición que se selecciona minimiza la varianza total de la estimación, dadas las decisiones anteriores. Esto es, si las posiciones x1, x2,…, xn se han seleccionado, la posición xn+1 se elige de la siguiente manera: usando el filtro de Kalman, calcula la varianza de la estimación en tiempo real que se obtendrá al agregar cada una de las posibles posiciones de monitoreo para la frecuencia fijada en los sensores; después escoge la posición que proporciona la varianza total mínima, y actualiza la matriz de covarianza en tiempo real, i.e., la matriz es actualizada sólo para los tiempos presentes y futuros con respecto a la posición espacio-temporal monitoreada por el sensor.
Determinación del número de posiciones de monitoreo
Para definir el número de posiciones de la red se propone usar el valor
Puede demostrarse que:
El criterio consiste en seleccionar las posiciones de monitoreo espacio-temporales donde se alcanza un 90% del valor de la máxima reducción posible (MRP):
Descripción de la zona de estudio
El acuífero Loreto se localiza en la porción sureste del estado de Zacatecas, entre las coordenadas 102° 51’-102° 08’ longitud oeste y 22° 12’-22° 37’ latitud norte; comprende una superficie aproximada de 743 km2. Limita al sur con el acuífero Valle de Chicalote del estado de Aguascalientes y al sureste con el acuífero Villa García; al noroeste, con el acuífero La Blanca; al oeste, con el de Ojocaliente, y al este con el acuífero Villa Hidalgo; estos últimos del estado de Zacatecas (figura 1).
Debido a la naturaleza y condiciones climáticas de la zona, las aguas subterráneas representan la única fuente segura de suministro para los sectores productivos, incluyendo el abastecimiento de agua potable.
El clima en la zona corresponde a semiseco, con una precipitación media anual de 415.4 mm, la cual ocurre principalmente en verano; la temperatura media anual registrada es de 16.6 °C y la evaporación potencial es de 2 019 mm por año.
Forma parte de la provincia fisiográfica de la Mesa Central, la cual se ubica entre dos grandes sierras de México, la Oriental y la Occidental; asimismo, queda dentro de la Subprovincia Llanos de Ojuelos-Aguascalientes.
Los arroyos existentes son de carácter torrencial, sólo tienen escurrimientos durante la época de lluvias y descargan en pequeñas lagunas, presas y bordos dispersados en toda el área.
De acuerdo con la información geológica, las rocas que afloran en la zona tienen una edad que va del Triásico Superior hasta depósitos del periodo reciente. Las más antiguas pertenecen a una secuencia volcanosedimentaria constituida por tobas andesíticas, cuarcitas y arcillas, afectadas por metamorfismo de bajo grado y parecen ser correlacionables con la formación Zacatecas del Triásico Superior marino. Sobreyaciendo a éstas se encuentra una secuencia marina de tipo flysch.
El acuífero está alojado en una depresión donde se depositan arcillas, arenas y fragmentos líticos de origen volcánico; la erosión provocó el depósito de conglomerados poco consolidados. En la parte superior está formado por sedimentos aluviales, cuyo espesor se va incrementando de los bordes a la porción central de la zona. En conjunto, estas unidades forman un acuífero con régimen libre y espesores que varían de 150 m en los bordes a 530 m en la parte central.
El modelo conceptual de flujo establece que la recarga principal proviene de la precipitación pluvial que se infiltra sobre las unidades permeables constituidas por aluviones, pie de monte, gravas, arenas, riolitas y calizas fracturadas. El agua infiltrada drena hacia la parte central de la zona.
En menor proporción, otro volumen proviene de los retornos de riego por bombeo, además de la recarga vertical, producto del agua precipitada que escurre y se infiltra en las partes topográficamente bajas de los valles. La descarga se realiza de manera artificial por bombeo y de manera natural un volumen pequeño por flujo horizontal, que está siendo drenado hacia el Valle de Chicalote, en el estado de Aguascalientes. La dirección preferencial de flujo subterráneo es de norte a sur, recibiendo aportaciones laterales de las zonas topográficamente altas que limitan al valle.
En el acuífero existe información piezométrica para el periodo 1980-2015. Se cuenta con una red piezométrica de monitoreo manual compuesta por 29 pozos piloto, que incluye información en temporada de estiaje y de lluvias; sin embargo, tan sólo 15 pozos distribuidos en el acuífero muestran registro continuo; en la mayoría de los casos sucede que los pozos utilizados sólo cuentan con algunos años de registros (Conagua, 2009).
Resultados y discusiones
Análisis geoestadístico espacio-temporal
Análisis exploratorio
Se realizó el análisis exploratorio tanto espacial como temporal de los valores de los niveles del agua subterránea, y se evaluaron sus medidas estadísticas de tendencia central y dispersión, para determinar la distribución de los datos y evaluar la presencia de valores atípicos. En el análisis se observó la presencia de 662 valores atípicos en el periodo comprendido entre 1980 y 2015; estos datos no fueron usados en los análisis subsecuentes. Se calcularon los estadígrafos de la variable analizada (nivel del agua) y se apreció un comportamiento no simétrico, lo que implica desviaciones de la muestra con respecto a la distribución normal (cuadro 1). Con base en los resultados de la curtosis y el coeficiente de simetría existen sesgos importantes en cuanto a la distribución normal, lo cual tiene un efecto de incremento en el error al momento de estimar la variable; sin embargo, esto no representa una restricción para el ajuste de un semivariograma espacio-temporal.
Datos | Media (m) | Mediana (m) | Varianza (m2) | Desviación estándar (m) | Mínimo (m) | Máximo (m) | Simetría | Curtosis | Fecha |
17 | 2 032.50 | 2 025.80 | 1 531.08 | 39.13 | 1 985.70 | 2 083.70 | 0.25 | 1.41 | Oct/1980 |
10 | 2 028.30 | 2 121.50 | 1 224.02 | 34.99 | 1 988.60 | 2 093.90 | 0.97 | 2.72 | Jul/1981 |
25 | 2 014.90 | 2 000.10 | 983.70 | 31.36 | 1 988.50 | 2 082.20 | 1.14 | 2.81 | Oct/1981 |
43 | 2 020.60 | 2 000.00 | 1 365.01 | 36.95 | 1 985.60 | 2 093.30 | 0.92 | 2.21 | Dic/1981 |
37 | 2 018.90 | 1 999.10 | 1 388.68 | 37.27 | 1 984.60 | 2 093.90 | 0.97 | 2.27 | Abr/1982 |
15 | 2 011.80 | 1 994.30 | 1 391.66 | 37.31 | 1 984.80 | 2 084.20 | 1.36 | 3.02 | Mar/1983 |
26 | 2 026.50 | 2 016.60 | 1 541.11 | 39.26 | 1 983.00 | 2 094.70 | 0.60 | 1.67 | Jul/1983 |
28 | 2 019.20 | 2 002.80 | 1 286.15 | 35.86 | 1 985.60 | 2 095.10 | 1.02 | 2.53 | Dic/1983 |
3 | 1 989.70 | 1 986.10 | 58.33 | 7.64 | 1 984.60 | 1 998.50 | 0.68 | 1.50 | Ago/1984 |
8 | 2 028.40 | 2 007.50 | 2 406.39 | 49.06 | 1 981.40 | 2 093.50 | 0.41 | 1.33 | Sep/1984 |
20 | 2 008.20 | 1 999.00 | 980.88 | 31.32 | 1 982.50 | 2 083.10 | 1.49 | 3.99 | Dic/1984 |
10 | 1 988.80 | 1 989.80 | 777.91 | 27.89 | 1 981.60 | 2 075.30 | 2.34 | 7.04 | Mar/1985 |
12 | 2 011.60 | 2 004.40 | 833.82 | 28.88 | 1 978.60 | 2 073.30 | 1.06 | 3.09 | May/1985 |
17 | 2 026.00 | 2 012.50 | 1 470.80 | 38.35 | 1 982.40 | 2 082.10 | 0.42 | 1.48 | Ago/1985 |
21 | 2 027.30 | 2 011.80 | 1 533.35 | 39.16 | 1 982.60 | 2 093.80 | 0.58 | 1.67 | Dic/1985 |
9 | 2 012.00 | 1 999.10 | 1 383.24 | 37.19 | 1 981.50 | 2 081.20 | 1.16 | 2.66 | Jun/1986 |
16 | 2 016.90 | 2 005.00 | 1 178.21 | 34.33 | 1 981.20 | 2 082.90 | 0.99 | 2.40 | Mar/1987 |
13 | 2 016.90 | 1 996.10 | 1 662.44 | 40.77 | 1 979.00 | 2 081.60 | 0.78 | 1.72 | Ago/1987 |
30 | 2 015.50 | 2 055.40 | 1 419.63 | 37.68 | 1 979.60 | 2 082.90 | 0.87 | 2.11 | Dic/1987 |
31 | 2 019.10 | 2 002.60 | 1 304.73 | 36.12 | 1 978.50 | 2 082.50 | 0.66 | 1.86 | Ene/1989 |
11 | 2 015.70 | 1 933.50 | 1 457.64 | 38.18 | 1 983.20 | 2 078.50 | 0.80 | 1.98 | Ene/1990 |
22 | 2 003.00 | 1 996.10 | 764.02 | 27.64 | 1 978.10 | 2 079.30 | 1.55 | 4.39 | Dic/1990 |
32 | 2 011.30 | 1 996.20 | 1 246.65 | 35.31 | 1 977.20 | 2 079.30 | 1.05 | 2.53 | Mar/1992 |
34 | 2 009.80 | 1 997.50 | 1 206.80 | 34.74 | 1 976.30 | 2 077.90 | 1.11 | 2.69 | Feb/1993 |
30 | 2 000.10 | 1 992.50 | 841.46 | 29.01 | 1 972.70 | 2 082.80 | 1.76 | 5.05 | Dic/1996 |
5 | 2 010.40 | 1 996.90 | 1 121.25 | 33.49 | 1 982.30 | 2 063.40 | 0.84 | 2.21 | Dic/1998 |
31 | 2 004.80 | 1 990.90 | 1 170.26 | 34.21 | 1 972.50 | 2 081.50 | 1.28 | 3.15 | Ene/1999 |
3 | 2 015.30 | 1 984.00 | 3 073.93 | 55.44 | 1 982.50 | 2 079.30 | 0.71 | 1.50 | Mar/2006 |
3 | 1 994.70 | 1 955.20 | 5 369.08 | 73.27 | 1 949.70 | 2 079.30 | 0.70 | 1.50 | Jul/2009 |
6 | 1 964.60 | 1 971.10 | 437.65 | 20.92 | 1 933.00 | 1 990.70 | -0.43 | 1.98 | Jul/2012 |
4 | 1 948.60 | 1 949.30 | 77.99 | 8.83 | 1 937.20 | 1 958.80 | -026 | 2.01 | Nov/2013 |
4 | 1 944.70 | 1 942.80 | 115.41 | 10.74 | 1 935.40 | 1 958.00 | 0.33 | 1.45 | Dic/2013 |
5 | 1 946.10 | 1 949.40 | 92.18 | 9.60 | 1 936.30 | 1 958.90 | 0.11 | 1.61 | Ene/2014 |
5 | 1 945.90 | 1 949.40 | 98.56 | 9.93 | 1 935.10 | 1 958.90 | 0.06 | 1.58 | Feb/2014 |
5 | 1 945.30 | 1 949.00 | 100.38 | 10.02 | 1 933.80 | 1 958.00 | -0.03 | 1.55 | Mar/2014 |
4 | 1 941.30 | 1 942.20 | 76.32 | 8.74 | 1 931.70 | 1 949.20 | -0.12 | 1.17 | Abr/2014 |
4 | 1 940.80 | 1 941.80 | 80.70 | 8.98 | 1 930.70 | 1 948.90 | -0.15 | 1.22 | May/2014 |
4 | 1 940.70 | 1 941.60 | 78.71 | 8.87 | 1 930.90 | 1 948.80 | -0.13 | 1.19 | Jun/2014 |
4 | 1 940.90 | 1 941.60 | 72.55 | 8.52 | 1 931.80 | 1 948.50 | -0.09 | 1.13 | Jul/2014 |
4 | 1 940.30 | 1 941.30 | 78.67 | 8.87 | 1 930.50 | 1 948.30 | -0.13 | 1.19 | Ago/2014 |
4 | 1 940.30 | 1 941.20 | 78.18 | 8.84 | 1 930.50 | 1 948.10 | -0.13 | 1.18 | Sep/2014 |
4 | 1 940.00 | 1 941.00 | 79.91 | 8.94 | 1 930.00 | 1 948.00 | -014 | 1.19 | Oct/2014 |
13 | 1 969.70 | 1 947.90 | 2 567.85 | 50.67 | 1 928.80 | 2 077.80 | 1.44 | 3.58 | Nov/2014 |
6 | 1 953.00 | 1 947.60 | 636.81 | 25.24 | 1 930.40 | 2 000.50 | 1.22 | 3.25 | Dic/2014 |
12 | 1 955.20 | 1 947.80 | 550.00 | 23.45 | 1 929.50 | 2 000.60 | 0.72 | 2.29 | Ene/2015 |
Análisis estructural
Este análisis consistió en determinar la estructura de los datos espacio-temporales. La base de datos para el análisis geoestadístico espacio-temporal incluye información de 102 pozos para el periodo 1980-2015 (figura 2), sin embargo no existen registros anuales completos en ninguno de ellos. Como primer paso se graficó el semivariograma espacio-temporal muestral, con un total de 650 datos piezométricos para diferentes incrementos de distancia y tiempo; en el cuadro 2 se muestran los intervalos espacio-temporales considerados, así como la separación en el espacio y tiempo. Para la estimación de los semivariogramas se utilizó una rutina modificada del programa GAMV escrito en Fortran para modelación espacio-temporal (De Cesare, Myers, & Posa, 2002).
Acuífero | Intervalos espacio-temporales | Separación en el espacio | Separación en el tiempo |
Loreto | 20 | 826.575 | 1.8026 |
Al analizar visualmente el comportamiento de los semivariogramas muestrales se advirtió la presencia de tendencia (figura 3); en este punto se decidió removerla mediante una función polinomial. Así, se procedió a trabajar con los residuos para determinar la estructura de los semivariogramas y después calcular los parámetros que los caracterizan. La tendencia fue removida con un polinomio de primer grado, para lo cual se empleó el método de regresión por mínimos cuadrados.
Una vez obtenidos los residuos de primer grado se modelaron los semivariogramas experimentales o muestrales para diferentes incrementos de distancia en el eje del espacio y lo mismo para el tiempo.
Efectuado lo anterior, se ajustó un modelo producto suma al semivariograma espacio-temporal muestral (De Iaco, Myers, & Posa, 2002), donde se conoce la cima de la componente espacial y la temporal, considerando la información de los últimos 35 años (figura 4).
El ajuste de los variogramas se definió a prueba y error, se evaluaron varios modelos y se seleccionó el que arrojó los mejores resultados, empleando el método de la validación cruzada. Los parámetros del modelo seleccionado se presentan en el cuadro 3.
Los resultados de la validación cruzada error medio (EM), raíz del error cuadrático medio (ECM) y el error cuadrático medio estándar (ECME) se presentan en el cuadro 4. El valor obtenido de EM es próximo a cero; ello indica que los errores positivos y negativos se cancelan. Para el caso del ECM, que representa la magnitud promedio del error cuadrático, en la estimación resulta un valor pequeño. Respecto al ECME, un valor cercano a la unidad refiere que la varianza obtenida por el modelo proporciona una buena medida de los errores de la estimación; dicho de otra manera, representa un buen ajuste.
Variable | EM | ECME | Raíz ECM | Diferencia mínima | Diferencia máxima |
Nivel del agua subterránea | 0.04 | 1.03 | 2.69 | -12.18 | 28.51 |
Una vez concluida la etapa del análisis geoestadístico espacio-temporal, se procedió a trabajar en el diseño de la red óptima de monitoreo automatizado. Los sitios de monitoreo considerados para seleccionar la red óptima son pozos que se encuentran fuera de operación o inactivos debido a que en ellos se pretende instalar sensores de medición automática de forma permanente, que registrarán las variaciones de los niveles del agua en forma continua.
Para la selección de los pozos a instrumentar se definió una malla que cubre las zonas donde se requiere la estimación del nivel del agua subterránea (figura 5). Se eligió una malla de 195 nodos, que forman elementos cuadrados separados 1.5 km en las direcciones X y Y, que cubren una superficie envolvente de la zona donde se ubican los pozos empleados en el análisis. Al aplicar el método se incluyeron las 195 posiciones de la malla de estimación y 81 pozos aptos para ser instrumentados. Para evaluar la información que aportarían los pozos seleccionados, se consideró que cada uno de esos sitios sería monitoreado durante los siguientes diez años.
Matriz de covarianza
La matriz de covarianza espacio-temporal a priori que requiere el filtro de Kalman fue calculada del modelo de semivariograma espacio-temporal obtenido en el análisis geoestadístico; se construyó para un total de 304 posiciones de estimación/monitoreo y 10 tiempos, por lo cual su dimensión es de 3 040 columnas y 3 040 renglones.
Filtro de Kalman
Mediante el filtro de Kalman se calculó la varianza del error en la estimación con base en la posición de los sitios de monitoreo, sin necesidad de conocer el dato medido en campo. De esta manera, el filtro de Kalman fue utilizado para determinar n sitios de monitoreo, al permitir evaluar cómo afecta añadir uno adicional a la varianza del error de la estimación resultante.
Se obtuvo una reducción de la varianza total cada vez que se incluía una posición en la red de monitoreo (figura 6).
La varianza total se reduce de manera considerable al elegir las primeras posiciones a instrumentar. Durante el proceso, el método asigna un orden de prioridad para cada uno de los pozos. De este modo, el orden de selección es un indicador de la importancia de los datos obtenidos en esas posiciones.
Selección de las posiciones óptimas
Como se explicó antes, el número de posiciones fue determinado con base en el criterio de seleccionar las posiciones de monitoreo espacio-temporal donde se alcanza un 90% del valor de la máxima reducción posible.
El 90% de la máxima reducción posible se alcanza cuando se han seleccionado 29 pozos para ser instrumentados (figura 7). Los resultados muestran que la raíz cuadrada de la varianza del error promedio del error en la estimación alcanza un valor máximo de 15.81 m antes de elegir alguna posición de monitoreo, y un mínimo de 9.01 m, lo cual indica que si se monitorearan las 81 posiciones, la raíz cuadrada de la varianza promedio del error en la estimación se reduciría en 6.80 m; sin embargo, cuando se seleccionan las primeras 29 posiciones, se reduce de 15.81 a 9.65 m; es decir, 6.16 m, significando que las últimas 52 posiciones sólo la reducen 0.64 m más, por lo cual se consideran como redundantes.
Varianzas del error en la estimación
En la figura 8 se presentan las varianzas del error en la estimación obtenidas para la red de monitoreo propuesta. Las varianzas tienen valores entre 9.99x10-8 y 208.63 m2, esto es, la magnitud del error en la estimación varía entre 3.16x10-4 (mínimo) y 14.44 m (máximo). Los valores mínimos se sitúan alrededor del centro del acuífero, donde se ubica la mayor parte de los pozos, y las máximas varianzas se localizan en el margen del área, debido a que allí sólo se localizan pocos pozos. Se pudo comprobar que las varianzas tienen la misma distribución espacial para todos los tiempos de monitoreo. Este tipo de mapas de varianza, combinados con criterios hidrogeológicos, pueden ser de gran utilidad en la elección de nuevas posiciones de monitoreo en las zonas donde se sitúan las máximas varianzas. Del análisis se considera que utilizar la red de pozos propuesta es suficiente para tener una buena estimación en el acuífero.
Conclusiones
Dado que la geoestadística emplea la correlación de la variable de estudio, es posible obtener una medida de la incertidumbre en la estimación de ésta; en la metodología aplicada para el diseño óptimo de redes de monitoreo automatizado de los niveles del agua subterránea se incorpora la correlación espacio-temporal entre los datos, utilizando un modelo de semivariograma espacio-temporal obtenido de un análisis geoestadístico.
Con el análisis geoestadístico se obtuvo el semivariograma muestral, el cual quedó bien definido debido a la cantidad de datos espacio-temporales. Le fue ajustado un modelo teórico esférico, porque representaba mejor la estructura de los datos. Los valores resultantes en la validación cruzada muestran que el modelo producto suma fue ajustado de forma adecuada a los datos espacio-temporales de los niveles del agua subterránea.
El filtro de Kalman resulta eficiente para reducir la varianza del error en la estimación. Tal característica lo convierte en una herramienta muy útil para el diseño de redes de monitoreo.
La aplicación de la metodología permitió asignar un orden de prioridad a cada uno de los pozos disponibles para ser instrumentados en función del nivel de representatividad que tienen sobre las estimaciones de los niveles del agua subterránea en una malla espaciotemporal definida para las zonas de interés.
Fue posible evaluar el nivel de información aportado por un sensor de medición de los niveles del agua subterránea durante el periodo en que es programado, dada una frecuencia de monitoreo previamente establecida.