Método para el diseño óptimo de redes de monitoreo de los niveles del agua subterránea

Method for the optimal design of networks to monitor groundwater levels

Jessica Vanessa Briseño-Ruiz, Graciela del Socorro Herrera-Zamarrón, Hugo Enrique Júnez-Ferreira
Universidad Nacional Autónoma de México

Dirección institucional de los autores

Universidad Nacional Autónoma de México Facultad de Ingeniería
Ubicado en el interior de las instalaciones del
Instituto Mexicano de Tecnología del Agua
Paseo Cuauhnáhuac 8532, Colonia Progreso
62550 Jiutepec, Morelos, México
Teléfono: +52 (777) 3293 600, extensiones 135 y 136
jevabrir@hotmail.com

Dra. Graciela del Socorro Herrera Zamarrón

Universidad Nacional Autónoma de México Instituto de Geofísica ]]> Ciudad Universitaria, Delegación Coyoacán
04510 México, D.F., México
Teléfonos: +52 (55) 5622 41 33 y 5622 4135
Fax: +52 (55) 5622 4198
ghz@geofisica.unam.mx

M.I. Hugo Enrique Júnez-Ferreira

Universidad Nacional Autónoma de México Facultad de Ingeniería
Ubicado en el interior de las instalaciones del
Instituto Mexicano de Tecnología del Agua
Paseo Cuauhnáhuac 8532, Colonia Progreso ]]> 62550 Jiutepec, Morelos, México
Teléfono: +52 (777) 3293 600, extensiones 135 y 136
hejunez@hotmail.com

Recibido: 27/04/10
Aprobado: 31/05/11

Resumen

Para estimar la disponibilidad del agua subterránea y su evolución es indispensable medir los niveles del agua subterránea por periodos largos de tiempo. Sin embargo, los recursos económicos son una limitante en la recopilación de suficientes datos que permitan estimar con certeza la evolución de dichos niveles. Por este motivo se requiere de métodos para el diseño de redes de monitoreo de los niveles del agua subterránea que faciliten utilizar los recursos económicos disponibles para el monitoreo de manera efectiva. El método de Herrera y Pinder (Herrera, 1998) para el diseño óptimo de redes de monitoreo de la calidad del agua subterránea, aplicado por Herrera et al. (2001), y también por Herrera y Pinder (2005), tiene varias ventajas. Por un lado, permite seleccionar tanto posiciones como tiempos de monitoreo en forma óptima; además, con la información obtenida de la red de monitoreo, permite hacer estimaciones en espacio y tiempo, y utiliza un modelo de flujo y transporte, el cual sintetiza el conocimiento de las características del acuífero, de la dinámica del flujo subterráneo y de la problemática de calidad del agua subterránea que se tenga, así como del conocimiento de las leyes de la física. En este trabajo, una modificación de esta metodología se prueba por primera vez para el diseño de redes de monitoreo de los niveles del agua subterránea. Se presenta un caso sintético en el que la red se diseña para un solo tiempo, suponiendo que el flujo está en estado estacionario. Los resultados presentados indican que la metodología propuesta por Herrera y Pinder para el diseño de redes de monitoreo tiene buenas perspectivas para aplicarse exitosamente al diseño de redes de monitoreo de los niveles del agua.

Palabras clave: redes de monitoreo, agua subterránea, diseño óptimo.

Abstract

It is essential to measure groundwater levels for long periods of time in order to estimate groundwater availability and its evolution. However, economic resources are a limitation to gathering enough data to estimate the evolution of these levels with certainty. Therefore, methods are required to design groundwater monitoring networks so as to effectively use the economic resources available for monitoring. The Herrera and Pinder method (Herrera, 1998) for the optimal design of groundwater quality sampling networks, applied by Herrera et al. (2001) and by Herrera and Pinder (2005), has several advantages: it allows for optimally selecting both monitoring positions and monitoring times; the information from the monitoring network enables obtaining estimates in space and time and; the method uses a flow and transport model that synthesizes knowledge about the characteristics of the aquifer, the groundwater flow dynamics and the groundwater quality problem, as well as knowledge about the laws of physics. This work tests, for the first time, a modification of the original method for the design of groundwater level monitoring networks. A synthetic case is presented in which the network is designed for a single time, assuming steady-state flow. The results indicate that the methodology proposed by Herrera and Pinder has good prospects for successful application to the optimal design of groundwater level monitoring networks.

Keywords: groundwater, monitoring networks, optimal design.

Introducción

El diseño de redes de monitoreo del agua subterránea consiste en establecer un programa de supervisión continua en pozos de muestreo para medir niveles del agua subterránea y tomar muestras para el análisis de variables físicas, químicas y biológicas previamente establecidas. Este tipo de redes de monitoreo piezométricas y de calidad del agua subterránea son indispensables para conocer el comportamiento de un acuífero, así como sus tendencias de calidad del agua. Esta información es básica para lograr una gestión eficaz y controlar los impactos ambientales negativos de la extracción del agua subterránea y de las cargas de contaminantes. Sin embargo, se requiere de un alto costo para medir los niveles del agua subterránea, recabar las muestras, y llevar a cabo el análisis de los parámetros de la calidad del agua subterránea en todos los pozos de monitoreo de un acuífero o de alguna zona contaminada. Es por ello que se requiere de un diseño óptimo de redes de monitoreo con las cuales se optimice la selección de los pozos y los tiempos de muestreo (en dónde, cuándo y cuántos pozos muestrear), los cuales proporcionen la mayor información para lograr estimaciones confiables de los parámetros de calidad del agua subterránea y de los niveles del agua subterránea.

Como respuesta a estas necesidades, Herrera (1998) propuso un método para el diseño óptimo de redes de monitoreo, que involucra en el diseño el espacio y el tiempo de forma combinada. Este método tiene varias ventajas: por un lado permite seleccionar tanto posiciones como tiempos de monitoreo en forma óptima; con la información obtenida de la red de monitoreo permite hacer estimaciones en espacio y el tiempo; y utiliza un modelo de flujo y transporte, el cual sintetiza el conocimiento de las características del acuífero, de la dinámica del flujo subterráneo y de la problemática de calidad del agua subterránea que se tenga, a la vez que incluye el conocimiento de las leyes físicas. Este método se probó en un caso de estudio que se reporta en Herrera et al. (2001), y otros ejemplos de aplicación se presentan en Herrera y Pinder (2005). Nos referiremos a este método como el método de Herrera y Pinder.

El método de Herrera y Pinder se propuso originalmente para diseñar redes de monitoreo de calidad del agua y se divide en dos procedimientos: 1) estima la varianza del error de la estimación que se obtendría si se tomaran muestras de agua en diferentes pozos y tiempos de monitoreo, y 2) escoge las posiciones y tiempos de monitoreo que minimizan una función de la varianza del error de la estimación, con los cuales se definen la red de monitoreo y su programa de muestreo. En este trabajo se prueba una modificación de ese método para diseñar redes de monitoreo de los niveles del agua subterránea. Dicha modificación consiste en emplear la carga hidráulica en vez de la concentración, y por tanto un modelo de flujo en vez de un modelo de transporte, como en los dos procedimientos descritos anteriormente. La formulación del método resultante se presenta para redes de monitoreo espacio-temporales y la aplicación se realiza en un caso de estudio sintético en el que la red de monitoreo se diseña para un solo tiempo.

A continuación se presenta un resumen de trabajos que proponen métodos para el diseño de redes de monitoreo del agua subterránea. Se incluyen también algunos que presentan métodos para el diseño de redes de monitoreo de la calidad del agua subterránea, ya que con algunas modificaciones pueden servir de base en el diseño de redes que permitan una adecuada caracterización de la carga hidráulica en acuíferos.

Tres puntos de vista han influenciado mucho el diseño de redes de monitoreo del agua subterránea (ASCE, 2003). En el primero de ellos, que llamamos marco hidrológico, siguiendo a Loaiciga et al. (1992), la red y su programa de muestreo se definen por consideraciones que toman en cuenta solamente las condiciones hidrológicas del sitio, sin recurrir a técnicas estadísticas o probabilísticas avanzadas. El segundo propone el análisis de datos dentro de un marco estadístico y define la red de muestreo basado en inferencias obtenidas de los datos, éste se llama marco estadístico. Al último punto de vista le llamaremos aquí marco de modelación, en donde se utilizan modelos matemáticos del agua subterránea para determinar posiciones y frecuencias de muestreo.

A continuación, primero se presentan los trabajos basados en criterios geoestadísticos, que pueden ser incluidos en el marco estadístico, y posteriormente trabajos que pueden incluirse en el marco de modelación. En cada una de estas categorías primero se describen los trabajos en los que únicamente se seleccionan de forma óptima las posiciones de muestreo y posteriormente los que incluyen también la selección óptima de frecuencias o tiempos de muestreo.

Marco estadístico

Rouhani (1985) propone el método de reducción de la varianza para aumentar una red de monitoreo existente. La red de monitoreo se aumenta agregando muestras de pozos en lugares que minimizan la varianza del error de la estimación de una variable, calculada por medio de kriging. Este trabajo no involucra la frecuencia temporal de muestreo en el diseño. La aplicación de la red de monitoreo se hace para la concentración de un contaminante, sin embargo, el mismo método se puede utilizar en el diseño de redes de monitoreo piezométrico.

Samper y Carrera (1990) presentan una comparación de diferentes métodos para el diseño de redes de observación basados en la implementación del kriging. Dividen el diseño de las redes de monitoreo en dos grandes categorías: cuando se añade un punto de medición a una red ya existente y cuando se añaden varios puntos. Incluimos aquí únicamente lo referente a la segunda categoría, que es la más relevante para este trabajo. En este caso se consideran varios métodos para añadir posiciones de puntos de medición:

• Método de inclusiones sucesivas: se van probando uno por uno puntos alternativos que se unen a un conjunto dado de pozos y se escoge el que minimiza la varianza del error de la estimación.

]]> • Intercambios sucesivos óptimos: se intercambian los elementos de un conjunto dado de puntos existentes con los elementos de un conjunto de puntos alternativos para encontrar la mínima varianza.

• Enumeración total: se analizan subconjuntos del total de pozos a monitorear; se empieza con todos los subconjuntos de un elemento, luego los que se formen con dos elementos y así sucesivamente hasta llegar al subconjunto que proporcione la incertidumbre deseada.

• "Rama y límite": inverso al anterior, el conjunto inicial corresponde al conjunto total de puntos disponibles y se van suprimiendo puntos.

Dentro de esta categoría también se considera la definición de redes de monitoreo mediante el uso de mallas regulares que pueden tener elementos con diversas geometrías (se ubican pozos de una red de monitoreo en cada uno de los elementos o nodos de una malla). En este caso se trabaja únicamente con los parámetros que definen la malla regular y se establece una relación entre la varianza y el costo asociado con estos parámetros.

Del análisis realizado para el diseño de una red se concluye que el método de inclusiones sucesivas, aunque es subóptimo, da buenos resultados en la elección de puntos de estimación, además de que resulta ser una metodología muy económica. Es por lo anterior que esta metodología se ha utilizado frecuentemente, obteniendo buenos resultados.

Cameron y Hunter (2000) propusieron un método para reducir la redundancia espacial y temporal por medio de dos algoritmos: uno temporal y otro espacial. El algoritmo temporal combina series de tiempo de datos de muchos pozos para construir un semivariograma temporal compuesto y usa el semivariograma para definir frecuencias de muestreo que proporcionen datos temporales sin correlación. En el algoritmo espacial se obtiene una estimación inicial de la pluma del contaminante usando kriging con todos los pozos existentes. Se asignan pesos numéricos, llamados pesos de kriging global, a la posición de los pozos, para evaluar su contribución relativa a la estimación de la pluma. Subconjuntos de pozos con los pesos de kriging global más bajos se remueven temporalmente y se hace una nueva estimación de la pluma. La variancia del kriging de la nueva estimación se compara con la de la estimación inicial para ver si la incertidumbre ha aumentado substancialmente. Si no es así, se quita el subconjunto de los pozos; en caso contrario, el subconjunto de pozos no se remueve. En este trabajo no se considera en el análisis la correlación cruzada espacio-tiempo de la concentración de contaminante.

Nunes et al. (2004) buscan optimizar redes de monitoreo del agua subterránea, considerando una reducción en la redundancia espacial y/o temporal. Se proponen tres modelos de optimización para seleccionar el mejor subconjunto de estaciones de una red de monitoreo del agua subterránea: 1) uno que maximiza la certidumbre espacial; 2) uno que minimiza la redundancia temporal, y 3) un modelo que maximiza la certidumbre espacial y minimiza la redundancia temporal. Los modelos de optimización propuestos se resuelven con recocido simulado, junto con un algoritmo de parametrización, usando entropía estadística. Los tres modelos se derivan de una ecuación que considera dos términos: uno espacial y otro temporal; los modelos empleados en el artículo resultan de simplificaciones a la ecuación objetivo. La ecuación general de la que se derivan los modelos contiene un término de varianza y un término que considera series de tiempo; las series de tiempo son representadas por funciones matemáticas comunes basadas en juicio empírico y experiencia.

Kumar et al. (2005) diseñaron una red de monitoreo de la carga hidráulica basándose en la teoría de las variables regionalizadas. Se le quita la tendencia a los datos al ajustarle un plano. Se ajusta un modelo lineal al semivariograma de los residuos obtenidos. Se estiman valores mediante kriging universal y se determinan los errores en la estimación mediante un procedimiento de validación cruzada (jacknife).

Faisal et al. (2007) optimizaron una red de monitoreo de la carga hidráulica para la estimación de un balance de aguas subterráneas en una cuenca de la India. Para la determinación del balance se utilizó el método de la fluctuación del doble espejo de agua (llamado así porque considera época de secas y de lluvias) propuesto por Maréchal et al. (2006). Se llevó a cabo un análisis geoestadístico para los datos recabados de carga hidráulica. El criterio de optimización de la red se basó en la validación cruzada del modelo de semivariograma propuesto. Se le asigna un orden de importancia a cada pozo de acuerdo con el tamaño del error promedio de tres años, resultado de la validación cruzada en cada sitio. Se realizó el cálculo del balance de aguas subterráneas para distintos subconjuntos de pozos (siguiendo el orden de importancia obtenido de la validación cruzada). Se seleccionó el número de pozos suficiente para no tener pérdida de información en la determinación del balance del agua subterránea.

Para el rediseño de una red de monitoreo de los niveles del agua en el acuífero de la Comarca Lagunera, México, Esquivel-Victoria (2008) utilizó un método geoestadístico propuesto por Ahmed (2004) para el diseño óptimo de redes de monitoreo de la temperatura del aire. En este método, el área de estudio se divide en una malla regular, y se obtienen la estimación y desviación estándar del error de la estimación para cada bloque utilizando kriging ordinario para los bloques con el total de datos existentes. La optimización se basa en la comparación de la desviación estándar del error de la estimación sobre los bloques, con un valor máximo preestablecido por el usuario, conforme se van quitando o añadiendo datos en la estimación.

Marco de modelación

Al parecer, Loaiciga (1989) fue el primer autor que propuso un método para el diseño espacio-temporal de redes de monitoreo del agua subterránea en el que la redundancia espacial y temporal de una red de monitoreo se analizaron en conjunto. El método combina kriging con una ecuación de transporte estocástica para predecir la incertidumbre del error de las concentraciones de un contaminante estimadas. Él usó el método de kriging generalizado para incluir el dominio en espacio y tiempo, y de esa manera evaluar cómo una muestra tomada de un pozo dado en un tiempo dado reduce la incertidumbre de todas las posiciones y tiempos en los que se obtendrían estimaciones. Para poder resolver el sistema de kriging es necesario contar con la matriz de covariancia espacio-temporal del error de la concentración en las posibles posiciones y tiempos de muestreo. Loaiciga deriva los elementos de la matriz de covariancia de la ecuación de advección-dispersión que gobierna el transporte de masa, discretizada numéricamente. Usó un término aditivo de error aleatorio en la ecuación de transporte que considera la incertidumbre del modelo; este término cambia de un paso de tiempo a otro en el modelo, y en cada nuevo paso es independiente de todos los errores previos. Para la optimización de la red, Loaiciga busca determinar dónde y cuándo muestrear, para minimizar la variancia del error de la concentración estimada (obtenida del sistema de kriging espacio-temporal) en las posiciones de interés, sujeto a constricciones de costo e insesgo. La determinación de un plan de muestreo óptimo se plantea como un problema de programación entera mixta.

Yangxiao et al. (1991) combinaron el procedimiento de estimación de parámetros propuesto en un trabajo previo (Van Geer et al., 1991) con un método para el diseño de redes de monitoreo de los niveles del agua subterránea. El objetivo del diseño de la red de monitoreo es observar los cambios espacio-temporales de los niveles del agua causados por la extracción de agua subterránea. Se utiliza un modelo de flujo determinista, al que se le suma un error incierto con correlación espacial. Los autores proponen utilizar un Filtro de Kalman para estimar los parámetros de la ecuación de flujo (llamados parámetros deterministas) y los del modelo de la matriz de covarianza del error (llamados parámetros estocásticos). La calibración se realiza para un periodo en el que las condiciones estacionarias se cumplen y en la que todas las matrices necesarias en el algoritmo del Filtro de Kalman no cambian con el tiempo. El método para la estimación de los parámetros se aplica de forma secuencial: primero se propone una estimación del modelo de la matriz de covarianza del error y la primera ronda de calibración de los parámetros deterministas se lleva a cabo; estos parámetros se mantienen constantes y se obtiene la primera ronda de calibración de los parámetros estocásticos. Estos dos pasos se repiten varias veces hasta que un criterio de error preestablecido para ambos conjuntos de parámetros se cumple. Se analizaron dos ejemplos de aplicación para el diseño de redes de monitoreo; en ambos, las frecuencias de muestreo se mantuvieron fijas y la densidad de la red se minimizó con la restricción de un valor umbral para la desviación estándar del error de la estimación. Se analizaron varias alternativas y la mejor fue elegida por inspección. Los autores discuten la importancia relativa de la densidad de la red espacial y la frecuencia de muestreo en relación con el tiempo de respuesta del sistema. Encontraron que si el sistema reacciona rápido, la optimización espacial de la red es importante. Si el sistema reacciona lentamente, tanto la optimización temporal como la espacial son importantes.

Herrera (1998) propuso una metodología para el diseño óptimo espacio-temporal de redes de monitoreo de la calidad del agua subterránea. El método combina un Filtro de Kalman estático con un modelo de transporte estocástico para predecir la varianza del error de las concentraciones estimadas. Una de las entradas para el Filtro de Kalman es la matriz de covarianza espacio-temporal de la concentración del contaminante en todas las posiciones y tiempo de muestreo posibles. Los elementos de esta matriz de covarianza se derivan de la ecuación de un modelo numérico de transporte de masa con advección y dispersión que gobierna el transporte de masa numéricamente discretizada. En este modelo, las variables aleatorias que se consideran son la conductividad hidráulica y la fuente del contaminante. Para optimizar la red se busca determinar dónde y cuándo muestrear para minimizar una función de la varianza del error de la concentración estimada (obtenida del Filtro de Kalman) en las posiciones y tiempos de interés. La determinación de un plan de muestreo óptimo se obtiene a través de un procedimiento secuencial, que selecciona los puntos de muestreo en espacio-tiempo que minimizan la función en cada paso, y se detiene cuando la función alcanza un valor predeterminado. Este método de optimización se puede considerar como uno de inclusiones sucesivas, que incluye tanto la selección de las posiciones de los pozos como los tiempos de monitoreo. En este trabajo se presentaron algunos ejemplos sintéticos de aplicación.

Posteriormente, Herrera et al. (2001) aplicaron la metodología propuesta por Herrera (1998) en el diseño de una red de monitoreo de la calidad del agua para el acuífero Toms River en New Jersey, Estados Unidos. Y en Herrera y Pinder (2005) incluyeron, además de la varianza, al coeficiente de variación como una medida de la incertidumbre. Como se mencionó antes, la metodología utilizada en el presente trabajo es una modificación de la propuesta por Herrera y Pinder.

Zhang et al. (2005) diseñaron una red de monitoreo óptima de la calidad del agua subterránea, combinando un Filtro de Kalman y un algoritmo genético para reducir al máximo el coeficiente de variación en posiciones y tiempos establecidos. La metodología propuesta es una modificación de la de Herrera y Pinder. En este caso, para diseñar redes de monitoreo de la calidad del agua: a) se generan realizaciones de la conductividad hidráulica usando la técnica Latyn Hypercube Sampling (LHS), suponiendo a la conductividad como un proceso aleatorio; b) se simula el flujo y transporte del agua subterránea para generar realizaciones de la concentración del contaminante, usando las realizaciones de la conductividad obtenidas en el paso anterior; c) se calcula una matriz de covarianza espacio-temporal de las concentraciones del contaminante; d) se usa un algoritmo genético para seleccionar posiciones y tiempos de muestreo que reduzcan el coeficiente de variación al máximo, en posiciones y tiempos de estimación. Se usa un Filtro de Kalman para actualizar la matriz de covarianza espacio-temporal después de que se ha tomado una muestra; e) se sigue con el proceso de selección hasta que se logran los valores objetivo del coeficiente de variación y se alcanza un escenario óptimo. Se realizó un problema ejemplo en Tomas River, New Jersey, Estados Unidos.

Hasta donde sabemos, no existen trabajos publicados en los que, al igual que en este trabajo, se utilicen modelos de flujo subterráneo que incluyan la incertidumbre de parámetros con correlación espacial para diseñar redes de monitoreo óptimas de los niveles del agua subterránea, ya que Herrera (1998) mostró que el modelo que utilizan Yangxiao et al. (1991) no es adecuado para representar los errores obtenidos como efecto de estos parámetros con correlación espacial.

Metodología de diseño óptimo espacio-temporal de una red de monitoreo piezométrica

1. Estimación de la carga hidráulica y de la incertidumbre de su error. Esta parte del método emplea un Filtro de Kalman y un modelo estocástico de flujo para obtener una estimación de h y la incertidumbre que el error de ésta tendría si datos de h de muestras tomadas en un conjunto de posiciones en diferentes tiempos se usaran para estimar la carga hidráulica. Primero se emplea el modelo estocástico de flujo para calcular una estimación inicial (o a priori) de h en diferentes posiciones y tiempos, así como su matriz de covarianza espacio-temporal. El modelo estocástico tiene la ventaja de que sintetiza el conocimiento que se tiene del acuífero y a la vez permite reconocer áreas con diferentes niveles de incertidumbre (a través de la varianza). Posteriormente, con el Filtro de Kalman se actualizan tanto la estimación inicial como la matriz de covarianza, al incluir los datos de piezometría disponibles y las posiciones de los pozos correspondientes. Es importante recalcar que para obtener la incertidumbre del error de la estimación no es necesario conocer el valor de los datos, sino únicamente las posiciones de las que se obtendrían los mismos; en este sentido, es posible predecir esta incertidumbre sin haber tomado las muestras.

2. Selección óptima de las posiciones de los pozos y los tiempos de muestreo. Una función de la varianza del error de la estimación, calculada con el procedimiento 1, se utiliza como criterio para escoger los pozos de muestreo y los tiempos de muestreo en cada pozo. La función que se utiliza depende de los objetivos del diseño. Para minimizar la función se pueden utilizar diferentes métodos de optimización; hasta ahora se ha usado un método secuencial que escoge los puntos espacio-temporales que la minimizan en cada paso y se detiene cuando el usuario lo decide.

Modelo estocástico

La ecuación de flujo se emplea para describir la evolución de h:

en donde h es la carga hidráulica [L]; K es la conductividad hidráulica [L/T]; Q representa fuentes o sumideros [L/T], y Ss es el coeficiente de almacenamiento específico [1/L]. Esta ecuación se discretiza usando un esquema numérico y el sistema se resuelve para el vector h de las incógnitas de la carga hidráulica en un subconjunto de puntos asociados con la malla numérica. A esta malla la llamaremos malla del modelo estocástico. Cualquiera de los parámetros, así como las condiciones iniciales o las de frontera, pueden representarse como variables aleatorias. En general se reconoce que uno de los parámetros con mayor incertidumbre en esta ecuación es la conductividad hidráulica. Por tal razón, en el ejemplo presentado en este trabajo, la conductividad hidráulica se representa como un campo aleatorio espacialmente correlacionado. Para caracterizar estadísticamente a K se calcula su media y se realiza un análisis geoestadístico, a fin de obtener un modelo de semivariograma que representa su estructura espacial.

El Filtro de Kalman es un conjunto de ecuaciones matemáticas que obtienen recursivamente estimaciones lineales no sesgadas y con varianza mínima para el estado de un sistema utilizando datos con ruido (Jazwinski, 1970). Asimismo, se puede emplear para actualizar la estimación en el momento en que se cuente con datos nuevos de la variable, sin necesidad de utilizar datos anteriores. El término recursivo significa que el filtro recalcula la solución cada vez que una nueva observación o medida es incorporada en el sistema.

En su forma general, el filtro se basa en dos ecuaciones: una ecuación dinámica y una ecuación de medición. La ecuación dinámica contiene un error aleatorio que no tiene correlación en el tiempo. Sin embargo, Herrera (1998) mostró que cuando se usa un campo aleatorio con correlación espacial para representar a la velocidad y se deriva con base en éste el modelo de transporte como una ecuación con coeficientes deterministas más un error aleatorio, el error puede tener una fuerte correlación en el tiempo. Por esta razón, el método que propuso Herrera (1998) emplea lo que llamamos el Filtro de Kalman estático, que únicamente emplea la ecuación de medición, incorporando el tiempo a través de la utilización de vectores espacio-temporales.

La ecuación lineal de mediciones del Filtro de Kalman discreto, que relaciona al vector h de la variable en las posiciones y tiempos en que se quiere estimar con las muestras z es:

en donde {z_i, j = 1, 2,...} es una secuencia de mediciones de la carga hidráulica. La matriz de H_j, es una matriz de 1 * N, que no es cero únicamente en la posición correspondiente a la entrada de h, en donde se toma el muestreo j, y N es la dimensión del vector h. H_j es la matriz del muestreo j; h = {h_ip} es el vector espacio-tiempo con los valores estimados de la carga hidráulica en las posiciones y tiempos de interés (h_ip es la carga hidráulica en la posición x_i en el tiempo t ); {v_j = 1, 2,...} es un escalar que representa el error de medición. Esta es una secuencia blanca Gaussiana, con media cero y varianza r_j. La secuencia del error de medición {v_j} y el vector h son independientes.

La matriz de covarianza del error de la estimación es:

en donde es el valor esperado de h, dadas las mediciones z₁, z₂,...,z_n y T denota transpuesta. En esta notación, el superíndice identifica el número " de mediciones que se usan para obtener la estimación.

En la aplicación del Filtro de Kalman, se estima h en una submalla de la malla del modelo estocástico, a la que llamamos malla del Filtro de Kalman. Esta malla tiene dos tipos de nodos, definidos como las posiciones de los posibles pozos de muestreo (S) y de los puntos de estimación (E), que se definirán cuando se explique el método de optimización de la red de monitoreo. Para la implementación del filtro se requiere proponer una estimación a priori de h en espacio y tiempo (), así como de la matriz de covarianza del error de estimación (P⁰). Dadas estas estimaciones a priori, la estimación lineal de mínima varianza para h se puede obtener secuencialmente a través de las siguientes fórmulas:

Estimación de los momentos a priori

El procedimiento para calcular la estimación a priori de h en espacio y tiempo (), así como la matriz de covarianza del error de esta estimación (P⁰) es la simulación estocástica o Monte Carlo (Zhang, 2002), y se explica a continuación.

1) Dadas la media y la de covarianza (o semivariograma) de la conductividad hidráulica se obtienen N realizaciones de esta variable que sean acordes con las mismas.

2) Se resuelve numéricamente la ecuación de flujo empleando cada realización, y se obtienen realizaciones de h para un periodo de tiempo.

3) Utilizando promedios de las realizaciones de h para cada nodo de la malla del Filtro de Kalman, obtenemos la media en espacio-tiempo, la que se propone como la estimación a priori; esta estimación es un vector al que denotaremos y de la misma manera se calcula su matriz de covarianza. La estimación a priori de h en espacio y tiempo (), así como la matriz de covarianza del error de la estimación (P⁰), se obtienen con las siguientes ecuaciones:

]]> en donde h^k_ip denota la k-ésima realización de h en la posición x_i en el tiempo t_p; es el vector que contiene la media para todas las realizaciones h^k_ip; Cov(e_ip, e_jq) es el elemento (i, p; j, q) de la matriz de covarianza del error de la estimación P⁰, donde el error e^k_ip = h^k_ip -; como es un vector en espacio-tiempo, éste incluye una estimación de h para todas las posiciones y tiempos de interés; por tanto, P⁰ es también una matriz de covarianza en espacio-tiempo.

Para el ejemplo que se presentará, el número total de entradas que integran la matriz de covarianza del error de la estimación de la carga hidráulica es Nt = (m*p)², donde p es el número total de tiempos de estimación y posibles tiempos de muestreo, y m es el número de nodos de la malla del Filtro de Kalman; m = NE + NS, donde NE es el número de puntos de estimación y NS es el número de pozos de muestreo, ya que no existe intersección entre los puntos de estimación y los de muestreo.

Función para minimizar la varianza del error de la estimación

Como se explicó antes, para la selección óptima de las posiciones de los pozos y los tiempos de muestreo se utiliza un método de optimización, en el que se minimiza una función de la varianza del error de estimación. En este caso, esta función es la suma de la varianza del error sobre todas las posiciones y tiempos de estimación. A esta función la llamamos la varianza total de la estimación de la carga hidráulica y se denota por σ²_T(n). Ésta, cuando se calcula con el plan de monitoreo que incluye n muestras, es:

en donde σ²_ip(n) es la varianza del error de la estimación en el i-ésimo punto de estimación y en el tiempo de estimación p, y E es el conjunto de los índices de los puntos de estimación. Las varianzas σ²_ip(n) se obtienen del Filtro de Kalman después de tomar n muestras, y son los elementos de la diagonal de la matriz de covarianza Pⁿ, que corresponden a los puntos y tiempos de estimación.

El algoritmo iterativo de optimización para la selección de las posiciones y tiempos de muestreo se explica diagramáticamente en la figura 1. La selección de las posiciones y tiempos de muestreo se hace en forma secuencial. El método utiliza un algoritmo iterativo de optimización que escoge una a la vez las posiciones de los pozos y las fechas de muestreo que minimizan la varianza total de la estimación, calculada como la suma de la varianza sobre todos los nodos y tiempos en que se va a estimar la h. En otras palabras, dadas las posiciones x₁,x₂,...,x_S de los posibles pozos de muestreo, y la matriz de covarianza a priori P⁰, en el primer paso del procedimiento, para x_i, i = 1,...s, con el Filtro de Kalman se calcula la varianza del error de la estimación que se obtendría al agregar el dato de cada posición de muestreo y se escoge la posición y_j^¿, que da la varianza total menor. A la matriz de covarianza resultante de aplicar el filtro con P⁰ y y_j^¿ la denotamos P_j^¿. Una vez que ya se seleccionaron las posiciones de los pozos, y_j^¿, j = 1,...,k para formar parte de la red óptima, se toman las posiciones que no han sido escogidas, y con la matriz de covarianza actualizada P_k^¿, que incluye los cambios debidos al muestreo de las y_j^¿, i = 1,...,k posiciones anteriores, se aplica el Filtro de Kalman y nuevamente se escoge la posición que da la varianza total menor. Las fórmulas empleadas para minimizar la varianza total se presentan en Herrera (1998), y Herrera y Pinder (2005).

Debido a que el proceso es secuencial, se realiza un análisis de la contribución en la reducción de la varianza de cada pozo para determinar cuándo parar el proceso y así quede definido el número total de pozos de la red de monitoreo óptima.

La metodología descrita se probó en un caso sintético basado en una representación simplificada del acuífero del Valle de Querétaro (AVQ). Como una primera prueba del método, se rediseñó una red de monitoreo de los niveles del agua subterránea, suponiendo que los cambios en los mismos son insignificantes. Para esto se utilizó un modelo numérico de flujo en estado estacionario.

El AVQ está localizado en la porción sur oriental del estado de Querétaro, México. Se ubica en el subsuelo de la ciudad de Querétaro y su zona conurbada (ZMCQ) (Mendoza y Herrera, 2007). Su población es la mayor de todo el estado, estimada en 962 240 habitantes en 2007 (Urquiza, 2008). La mayor explotación de agua subterránea en el AVQ es para el uso público urbano, con la cual se abastece la ZMCQ.

Objetivo de la red de monitoreo

El objetivo del diseño de la red de monitoreo de los niveles del agua subterránea consistió en seleccionar, de entre 33 pozos preexistentes en el AVQ, aquellos que no dieran información redundante para la estimación de h en todo el acuífero. De acuerdo con este objetivo se estimará h en todo el acuífero, dándole el mismo peso a todas las zonas. A continuación se explica cómo se aplicó el método.

Simulación secuencial Gaussiana

Con base en un análisis geoestadístico de 46 datos de K existentes en la zona de estudio, se ajustó un modelo exponencial para el semivariograma muestral de Ln K:

en donde γ LnK es el semivariograma de LnK; c₀, la pepita o nugget; c₀ + c₁ el sill o la meseta del semivariograma; d, la distancia, y a es la escala de correlación. Los valores ajustados para los parámetros fueron c₀ = 0.43, c_l = 0.57, a = 8 900 m y la distancia máxima considerada es de 1 262.05 metros.

Modelo determinista y estocástico de flujo de agua subterránea

La simulación de flujo se realizó con un simulador determinista convencional de elemento finito llamado Princeton Transport Code (Babu et al., 1993). Éste se adaptó utilizando el modelo calibrado por Simuta (2005). Las características del modelo del AVQ se pueden ver en la figura 2. La zona de modelación corresponde a la del acuífero del Valle de Querétaro e incluye también al Valle de Obrajuelo, en el estado de Guanajuato, con un total de 426.87 km². Este último acuífero se incluye debido a que aunque por razones administrativas la Comisión Nacional del Agua incluye al Valle de Obrajuelo como parte del acuífero de Celaya, Simuta (2005) concluye que existe conexión hidráulica entre ellos y es necesario modelarlos en forma conjunta. La malla del modelo es triangular y cuenta con 6 676 nodos, 13 066 elementos y una capa. Los valores de conductividad hidráulica en cada realización se asignaron a los nodos de la malla numérica del modelo determinista, encontrando el nodo más cercano al centro de la malla de elementos cuadrados utilizada para las realizaciones de K. Los parámetros hidráulicos del modelo son: el almacenamiento específico es uniforme en toda la zona de modelación (Ss = 0.001). La recarga se considera variable en la zona con rangos de 0.011 a 0.067 metros/día para la recarga vertical originada por la precipitación y de 2.45 e^-05 a 0.0003229 metros/día en la recarga originada por el riego. En el caso de estudio no se considera extracción por bombeo. Se tienen cinco zonas con condiciones de frontera de carga asignada, en promedio cien metros por debajo de la topografía.

El modelo estocástico de flujo tiene la misma malla que el modelo determinista de flujo. El diseño de la malla cumple el criterio sugerido por Ababou (1998) para que las simulaciones del LnK no tengan una distorsión significativa debido a la discretización. Éste es, donde Δx es el tamaño de los elementos de la malla y λ, es la escala de correlación que para el modelo exponencial es a/3. Las condiciones de frontera, así como todos los parámetros del modelo se tomaron como deterministas, con excepción de K, que es considerada como un campo aleatorio con media de 1.63 metros/día y correlación espacial dada por el semivariograma exponencial descrito en la sección anterior. Para obtener la solución del modelo estocástico se generaron cuatro mil realizaciones de K y sus correspondientes realizaciones de h.

Las posiciones en las que se estimará h con los datos que se obtendrán de la red de monitoreo se asocian con los nodos de lo que llamamos malla de estimación. De acuerdo con lo explicado anteriormente, en estos puntos de estimación es en donde se minimiza la varianza del error de estimación de h en el proceso de optimizar la red de monitoreo. De acuerdo con el objetivo de la red de monitoreo, para darle el mismo peso a toda el área del acuífero en la optimización se propuso como malla de estimación una malla compuesta por 66 puntos de estimación aproximadamente equi-espaciados sobre el AVQ (figura 3). Las estimaciones de h se calcularon en los nodos de la malla de estimación, y por lo tanto la matriz de covarianza incluyó únicamente estos nodos y las posiciones de los 33 posibles pozos de monitoreo, que también se pueden ver en la figura 3.

Para verificar que el número de las realizaciones del modelo estocástico fueran suficientes, se estableció el criterio de convergencia para la matriz de covarianza basado en las siguientes ecuaciones:

en donde P_m es la matriz de covarianza calculada con m realizaciones, Δ es el incremento en el número de las realizaciones y N es el número de entradas de la matriz de covarianza. En la prueba, el incremento A utilizado fue de 20. Con base en estos criterios se considera que cuatro mil realizaciones son suficientes (ver figura 4) y se toma como matriz de covarianza a priori para continuar con el diseño de la red óptima de monitoreo piezométrico a la realización número cuatro mil.

Resultados

La figura 5 muestra la varianza del error de la estimación inicial, ésta es la varianza de h obtenida de la simulación estocástica. Como se puede observar, tiene valores entre 0 y 9 m², con valores mayores hacia el límite oeste de la zona en la que se obtendrán las estimaciones y en la región noreste también se tiene un máximo local. Esto indica que la carga hidráulica es más sensible a los cambios en la conductividad hidráulica en esas regiones.

Como se mencionó anteriormente, el método de optimización propuesto es secuencial y se añade un pozo de monitoreo en cada paso. Esto produce los resultados que se muestran en forma de gráfica en la figura 6, en donde se representa la varianza total en los puntos de estimación contra el número de pozos incluidos en la red. Se observa que los primeros pozos elegidos proporcionan mucha información y la varianza total se reduce significativamente; conforme se van añadiendo pozos a la red de monitoreo, la información obtenida por cada pozo va disminuyendo. La varianza total inicial (m²) para los 66 puntos de estimación es de cerca de 120, y conforme se van añadiendo pozos a la red de monitoreo, esta varianza disminuye hasta un valor cercano a 40. Esto indica que al utilizar todos los pozos de la red de monitoreo, la varianza total no baja de este valor, que al promediar sobre los 66 puntos de estimación da una varianza cercana a 0.6. Sin embargo, esta figura también hace ver que los últimos pozos contribuyen muy poco a reducir la varianza, por lo que es natural analizar la pertinencia de dejarlos en la red de monitoreo.

Por la forma en la que se construye la red de monitoreo, el orden en el que se eligen los pozos indica el grado de importancia de cada pozo en la disminución de la varianza; por ese motivo, a éste le llamamos orden de importancia de los pozos. El orden de importancia de los pozos se puede ver en la figura 7. Como se puede observar, el pozo con orden de importancia 1 se encuentra en la parte oeste de la zona establecida para optimizar la red de monitoreo; ello es lógico, ya que ésta es la región con mayor varianza. De igual forma, el pozo con orden de importancia 2 se ubica al noroeste de la zona de estudio, en donde se señaló anteriormente que existe un máximo local. Es importante recalcar que al seleccionar estos pozos se toma en cuenta la varianza del error de la estimación en toda la zona de influencia del pozo y no únicamente en el punto que se escoge. De esta manera, se escogen también los pozos con orden de importancia 3, 4 y 5, en posiciones en las que se tienen varianzas altas, y en las que no se han escogido pozos anteriormente. A partir del pozo con orden de importancia 6, van cobrando más importancia las selecciones realizadas anteriormente y las zonas de influencia de los pozos (zonas con correlaciones altas obtenidas de las matrices de covarianza resultantes de aplicar el Filtro de Kalman conforme se van añadiendo los pozos), y se inicia la selección de pozos cercanos a otros ya elegidos o en la zona central, que inicialmente tiene varianzas pequeñas.

Análisis para determinar el número total de pozos en la red de monitoreo

Como se comentó antes, en la optimización de esta red de monitoreo el objetivo es determinar los pozos que no sean redundantes. Para hacer esto se lleva a cabo un análisis en el que se produce una estimación inicial y se actualiza con los datos asociados con una red de monitoreo con un determinado número de pozos, mediante la aplicación del Filtro de Kalman. En el caso de estudio sintético, este procedimiento se realizó con el objeto de mostrar cómo se haría este análisis en un caso real, y para ello se tomó al azar una de las realizaciones de h del modelo estocástico, como si fuera la carga hidráulica en campo que se quiere estimar.

Como un primer paso, se muestran los valores de h en los 6 676 nodos de la malla del modelo (figura 8a). También se muestran las estimaciones de producidas con el Filtro de Kalman en los 66 nodos de la malla de estimación utilizando 12 (figura 8b), 22 (figura 8c) y 33 (figura 8d) datos de h que corresponden a los pozos con el orden de importancia del 1 hasta el número indicado. La estimación se realizó únicamente para la zona del estado de Querétaro (el Valle de Obrajuelo no se incluyó en la estimación). La comparación entre estas piezometrías indica, por un lado, que la estimación generada con el filtro, empleando 12, 22 y 33 datos de los pozos de muestreo de la red propuesta, generan piezometrías acordes con la realización escogida de h; por otro lado, también las piezometrías generadas con 12, 22 y 33 datos son muy similares.

En segundo lugar, se realizó una comparación entre la varianza inicial del modelo (figura 5a) y las varianzas finales obtenidas con el Filtro de Kalman después de incorporar los datos de 12, 22 y 33 pozos de muestreo de la red de monitoreo propuesta. Las figuras 5b, 5c y 5d muestran las varianzas que se obtienen del Filtro de Kalman después de muestrear 12, 22 y 33 pozos. La varianza total después de muestrear 33 pozos, disminuyó considerablemente y tiene valores de entre 0 y 3. También se puede ver que las varianzas del error de la estimación usando 22 pozos y usando 33 pozos son muy similares, lo que confirma que existe una redundancia importante entre los 11 pozos que se añaden.

Adicionalmente, se calcularon los errores en los 66 puntos de estimación, restando a la realización escogida la estimación producida con el filtro. Los errores más grandes detectados se encuentran ubicados en las zonas NE y SW y varían de -1.29 a 1.93 metros, y se ubican en las zonas donde se encuentran las varianzas más grandes; ello se debe a que en estas zonas no se cuenta con ningún pozo de muestreo que aporte datos para reducirlas. Esto indica que se requerirían pozos adicionales en tales zonas para obtener una mejor estimación.

Para comparar las alternativas entre la opción de una red de monitoreo con 12, 22 o 33 pozos, se realizó el cálculo de los errores. El error e. se calculó mediante la diferencia de los datos de la realización de h y de la h estimada mediante el Filtro de Kalman en los puntos de la malla de estimación con cada una de las opciones de números de pozos. Los errores calculados son el error medio (EM), el error cuadrático medio (ECM) y el error cuadrático medio estándar. El EM es un indicador de la sobrestimación o subestimación de los valores de la realización de h. El error cuadrático medio (ECM) es una medida de la magnitud del error de estimación. La teoría indica que mientras más cercanos estén el EM y el ECM a cero, más adecuado será el ajuste entre la realización de h y la estimación de con 12, 22 o 33 pozos (Armando, 1994). El error cuadrático medio estándar (ECME), mientras más cercano esté a 1 indica el nivel de confianza entre la realización de h y los valores de la estimación h. Los resultados de estos errores para cada uno de los casos de estudio se muestran en cuadro 1. Se puede observar cómo el ECM y ECME tienen valores cercanos a lo que indica la teoría, sin embargo, es importante notar que cuando el número de pozos que conforman la red disminuye, también decrece el ECME.

Red de monitoreo final

Del análisis anterior se concluye que al utilizar una red de 22 pozos se obtiene una estimación con errores similares a los obtenidos con los 33 pozos, esto indica que los once pozos con orden de selección de 23 a 33 se pueden considerar como redundantes. Además, los resultados de la piezometría indican que con 22 pozos obtenemos una representación piezométrica muy parecida a la real, con errores de estimación del orden de 1.92 a -1.33 metros, y con un error cuadrático medio estándar de 0.97, muy cercano a 1. Es importante tener en cuenta que el número de pozos que conformarán la red de monitoreo final estará definida por la autoridad federal, estatal o municipal de aguas, o en su caso por los contaminadores potenciales del agua subterránea, como empresas que requieran de algún monitoreo de calidad del agua en sus instalaciones. La red está en función del tipo y volumen de datos requeridos, y varía considerablemente en cuanto al aspecto de gestión de que se trate, pero inevitablemente también dependerá de los recursos financieros disponibles y del nivel de riesgo que se esté dispuesto a correr.

Conclusiones

Del análisis anterior se concluye que el diseño de la red de monitoreo piezométrica para este caso de estudio mediante el método de Herrera y Pinder es satisfactorio, y que se tienen buenas perspectivas para aplicar exitosamente esta metodología al diseño de redes de monitoreo de los niveles del agua subterránea en casos reales.

Como se mostró en los resultados del análisis de la varianza en la figura 6, el modelo estocástico de flujo es una valiosa herramienta para evaluar la incertidumbre en las estimaciones del modelo, haciendo la metodología muy efectiva en la selección de las posiciones y los tiempos de muestreo que minimizan la incertidumbre de la estimación. Asimismo, la metodología es capaz de sugerir en qué zonas se requiere la construcción y/o habilitación de pozos de monitoreo que aporten datos de carga hidráulica para lograr una mejor estimación en la zona de estudio.

Es importante recalcar que la formulación del método se hizo para seleccionar en forma óptima tanto posiciones como tiempos de monitoreo, y la aplicación se hizo para la selección de posiciones de monitoreo únicamente. En un trabajo futuro se planea realizar una prueba similar para el diseño espacio-temporal de una red de monitoreo piezométrica.

Finalmente es importante hacer notar que al aplicar este tipo de métodos a casos reales, además de utilizar criterios de optimización, es indispensable usar criterios geohidrológicos, que dependen del conocimiento del área de estudio.

Agradecimientos

Jessica Briseño y Hugo Júnez agradecen al Consejo Nacional de Ciencia y Tecnología (Conacyt) por el otorgamiento de sus becas para la realización de estudios de doctorado, lo que permitió el desarrollo de la investigación como parte de sus tesis doctorales.

Referencias

ABABOU, R. Three-dimensional flow in random porous media. Ph.D. diss. Cambridge: Massachusetts Institute of Technology, 1988.         [ Links ]

AHMED, S. Geostatistical estimation variance approach to optimizing an air temperature monitoring network. Water, Air and Soil Pollution. Vol. 158, 2004, pp. 387-399.         [ Links ]

ARMANDO, A.J. Introducción al tratamiento de series temporales: aplicación a las ciencias de la salud. Ediciones Díaz de Santos, 1994, 585 pp.         [ Links ]

ASCE. Long-Term Groundwater Monitoring: The State of the Art. Task Committee on the State of the Art in Long-term Groundwater Monitoring Design, American Society of Civil Engineers, USA, 2003.         [ Links ]

BABU, D., PINDER, G.F., NIEMI, A., AHLFELD, D.P., and STOTHOOFF, A. Chemical transport by three dimensional groundwater flows (PTC-Princeton Transport Code). Rep. 84-WR-3. Princeton: Princeton University, Department of Civil Engineering, 1993.         [ Links ]

CAMERON, K. and HUNTER, P. Optimization of LTM Networks Using GTS: Statistical Approaches to Spatial and Temporal Redundancy. Brooks, USA: Tech. Rep., Air Force Center for Environmental Excellence, 2000.         [ Links ]

DEUTSCH, C.V. and JOURNEL, A.G. GSLIB Geostatistical Software Library and User's Guide. 2^nd. ed. New York: Oxford University Press, 1998, 368 pp.         [ Links ]

ESQUIVEL-VICTORIA, E. Aplicación geoestadística al rediseño de la red de monitoreo en el acuífero de la Comarca Lagunera, Coahuila, Durango, México. Tesis para obtener el grado de maestría. México, D.F.: Universidad Nacional Autónoma de México, 2008, 53 pp.         [ Links ]

FAISAL, K., ZAIDI, S., and DEWANDEL, B. Optimizing a piezometric network in the estimation of the groundwater budget: a case study from a crystalline-rock watershed in southern India. Hydrogeology Journal. Vol. 15, No. 6, 2007, pp. 1131-1145.         [ Links ]

HERRERA, G.S. Cost Effective Groundwater Quality Sampling Network Design. Ph. D. thesis. Vermont: University of Vermont, Burlington, 1998, 172 pp.         [ Links ]

HERRERA, G.S., GUARNACCIA, J., PINDER, G. y SIMUTA, R. Diseño de redes de monitoreo de la calidad del agua subterránea eficientes. Proceedings of the 2001 International Symposium on Environmental Hydraulics, Arizona, USA, 2001.         [ Links ]

HERRERA, G.S. and PINDER, G.F. Space-time optimization of groundwater quality sampling networks Water Resources Research. Vol. 41, No. 12, W12407, 10.1029/2004WR003626. 2005.         [ Links ]

JAZWINSKI, A.H. Stochastic Processes and Filtering Theory. San Diego: Academic Press, 1970.         [ Links ]

KUMAR, S., SONDHI, S.K., and PHOGAT, V. Network design for groundwater level monitoring in Upper Bari Doab canal tract, Punjab, India. Irrigation and Drainage. Vol. 54, 2005, pp. 431-442.         [ Links ]

LOAICIGA, H.A. An optimization approach for groundwater quality monitoring network design. Water Resources Research. Vol. 25, No. 8, 1989, pp. 1771-1782.         [ Links ]

LOAICIGA, H.A., CHARBENNEAU, R.J., EVERETT, L.G., FOGG, G.E., and HOBBS, B.F. Review of groundwater quality monitoring network design. Journal Hydraulic Engineering. ASCE. Vol. 118, No. 1, 1992, pp. 11-32.         [ Links ]

MANCHUK, J., LEUANGTHONG, O., and DEUTSCH C. A New Approach to Direct Sequential Simulation that Accounts for the Proportional Effect: Direct Lognormal Simulation. Sixth Annual Report of the Centre for Computational Geostatistics. University of Alberta. 2004. [citado el 13 de abril del 2010] Disponible para World Wide Web: http://www.uofaweb.ualberta.ca/ccg/pdfs/2004%20104-Direct%20Lognormal%20Simulation.pdf        [ Links ]

MARÉCHAL, J.C.,DEWANDEL,B.,AHMED,S., GALEAZZI, L., and ZAIDI, F.K. Combining the groundwater budget and water table fluctuation methods to estimate specific yield and natural recharge. Journal of Hydrology. Vol. 329, No. 1-2, 2006, pp. 281-293.         [ Links ]

MENDOZA, E.Y. y HERRERA, G. Estimación multivariada espacio-tiempo de la carga hidráulica en el Valle de Querétaro-Obrajuelo. Ingeniería hidráulica en México. Vol. XXII, núm. 1, enero-marzo de 2007, pp. 63-80.         [ Links ]

NUNES, L.M., CUNHA, M.C., and RIBEIRO, L. Groundwater monitoring network optimization with redundancy reduction. J. Water Resour. Plng. and Mgmt. Vol. 130, No. 1, 2004, pp. 33-43.         [ Links ]

ROUHANI, S. Variance reduction analysis. Water Resources Research. Vol. 21, No. 6, 1985, pp. 837-846.         [ Links ]

SAMPER, F.J. y CARRERA J. Geoestadística, aplicaciones a la hidrogeología subterránea. Barcelona: Centro Internacional de Métodos Numéricos en Ingeniería, 1990.         [ Links ]

SIMUTA, R. Modelo en elemento finito para el flujo del acuífero del valle de Querétaro. Tesis de maestría. México, D.F.: Universidad Nacional Autónoma de México, 2005.         [ Links ]

URQUIZA, M. Uso sostenible del agua en la ciudad de Querétaro [en línea]. Expo Zaragoza 2008. Semana temática: "Agua y servicios de abastecimiento y saneamiento". Eje temático: Foro de los buenos ejemplos. Zaragoza, España, de junio a septiembre del 2008 [citado el 9 de marzo del 2010]. Disponible para World Wide Web: http://www.zaragoza.es/contenidos/medioambiente/cajaAzul/27S6-P5-Manuel%20UrquizaACC.pdf.         [ Links ]

VAN GEER, F., TE STROET, C., and YANGXIAO, Z. Using Kalman Filtering to Improve and Quantify the Uncertainty of Numerical Groundwater Simulations 1. The Role of System Noise and Its Calibration. Water Resources Research. Vol. 27, No. 8, 1991, doi:10.1029/91WR00509.         [ Links ]

YANGXIAO, Z., TE STROET, C., and VAN GEER, F. Using KalmanFilteringtoImproveand Quantify theUncertainty of Numerical Groundwater Simulations 2. Application to Monitoring Network Design. Water Resources Research. Vol. 27, No. 8, 1991, doi:10.1029/91WR00510, 1991.         [ Links ]

ZHANG, D. Stochastic Methods for Flow in Porous Media. Coping with Uncertainties. San Diego: Academic Press, 2002.         [ Links ]

ZHANG, Y., PINDER, G.F., and HERRERA, G.S. Least cost design of groundwater quality monitoring networks. Water Resources Research. Vol. 41, W08412, doi:10.1029/2005WR003936, 2005.         [ Links ]