Estimación espacio-temporal de la carga hidráulica utilizando el concepto de función aleatoria espacio-tiempo

Spatiotemporal estimation of hydraulic head using a single spatiotemporal random function model

Edgar Yuri Mendoza-Cázares¹, Graciela del Socorro Herrera-Zamarrón²

¹ Instituto Mexicano de Tecnología del Agua.

² Universidad Nacional Autónoma de México.

Dirección institucional de los autores

Dr. Edgar Yuri Mendoza-Cázares

Instituto Mexicano de Tecnología del Agua
Coordinación de Hidrología
Subcoordinación de Aguas Subterráneas
Paseo Cuauhnáhuac 8532
Colonia Progreso, 62550 Jiutepec, Morelos, México
teléfono: + (52) (777) 329 3600, extensión 339
fax: + (52) (777) 329 3682 ]]> edgar_mendoza@tlaloc.imta.mx

Dra. Graciela del Socorro Herrera-Zamarrón

Universidad Nacional Autónoma de México
Ciudad Universitaria
Delegación Coyoacán
04510 México, D.F.
teléfonos: + (52) (55) 5622 4133 y 56 22 41 35
fax: + (52) (55) 5622 4198
ghz@geofisica.unam.mx

Resumen

En este trabajo se presenta un método geoestadístico para realizar estimaciones espacio-temporales y se aplica a la estimación de la carga hidráulica del acuífero de Querétaro-Obrajuelo en el periodo 1981-2004. Las estimaciones para los años 1993, 1995 y 1999 se comparan con las de los métodos de krigeado ordinario y krigeado multivariado, utilizando el método de validación cruzada, reportándose los valores del error medio (EM), error cuadrático medio (ECM) y error cuadrático medio estándar (ECME). El promedio del error medio para los tres años elegidos, para el caso de krigeado ordinario, resulta ser el menor (-0.23); el valor promedio del error cuadrático medio (ECM) es más bajo en el caso del krigeado espacio-temporal (224.29), y el valor promedio del error cuadrático medio estándar (ECME) es mejor para el método de cokrigeado, ya que sus valores de ECME son cercanos a la unidad (0.95). El ECME para el método espacio-tiempo es de 0.8 al considerar todos los tiempos; pero en lo particular para los años elegidos, este índice disminuye. Los resultados de las varianzas del valor estimado son siempre más pequeños utilizando el método espacio-temporal, ya que utiliza más información para la estimación; asimismo, es posible hacer las estimaciones en todo el espacio para todos los tiempos. Se concluyó que la herramienta es poderosa, porque considera toda la información disponible al efectuar las estimaciones.

Palabras clave: modelo espacio-temporal, geoestadística, estimación carga hidráulica.

Abstract

This paper presents a geostatistical method for performing estimates using space-time kriging which is applied for the estimation of the water level of the Queretaro-Obrajuelo aquifer in the 1981-2004 period. The estimates obtained by this method for the years 1993, 1995 and 1999 are compared with ordinary kriging and cokriging using cross validation. The average mean error (ME) for the three chosen years was lowest for ordinary kriging (-0.23), the average mean squared error (MSE) is the lowest in the case of the space-time method (224.29) and the value of the average squared mean standard error (SMSE) is better for the method of cokriging (0.95), because the SMSE values are close to unity. The SMSE for the space-time method is 0.8 when considering all the times, but it decreases in particular for the selected years. The results of the estimated variances are always smaller using the space-time method, because it uses more information to get the estimate. Also, it is possible to make estimates in the whole space for all times. Therefore, it was concluded that the tool is powerful, since it considers all available information to produce the estimates.

Keywords: single spatio-temporal model, geostatistics, hydraulic head estimation.

Introducción

Los métodos geoestadísticos de krigeado han sido ampliamente utilizados en la hidrogeología subterránea. Numerosos investigadores usan los métodos geoestadísticos para estimar campos de transmisividad (Ahmed y De Marsily, 1987) la carga hidrostática de un acuífero (Delhomme, 1978; Gambolati y Volpi, 1979a y 1979b; Aboufirassi y Mariño, 1983; Rouhani, 1986; Gambolati y Galeati, 1987) o la precipitación (Hevesi et al., 1992a y 1992b).

La constante en estos trabajos es que la variable se estima en el dominio del espacio; sin embargo, se sabe que las variables son dinámicas, lo cual indica que tienen una componente temporal de importancia. En algunos casos, los cambios son lentos (e.g.: transmisividad) y en otros son rápidos (e.g.: carga hidráulica).

De forma tradicional, se ha realizado el análisis bajo dos concepciones: el punto de vista espacial (ASCE, 1992a y 1992b) y el punto de vista temporal (Asli y Marcotte, 1995; Larocque et al., 1998; Angelini, 1997; Lee y Lee, 2000). No obstante, la conceptualización conjunta a la que nos referiremos como espacio-temporal ha sido ámbito de estudio en los últimos 15 años y poco estudiada en el caso particular de la hidrogeología subterránea (Christakos, 2000).

Rouhani y Hall (1989) son los primeros investigadores que plantean determinar la carga hidráulica utilizando información conjunta en el dominio del espacio y el tiempo. Los autores muestran cómo calcular dicha estimación usando una modelo del tipo lineal:

donde C_st(h_sh_t) es la función de covarianza espacio-temporal; C_s(h_s), C_t(h_t) son las funciones de covarianza espacial y covarianza temporal, respectivamente. El modelo lineal presenta problemas, ya que la matriz de covarianza, para ciertos casos, es singular (Myers y Journel, 1990; Rouhani y Myers, 1990; Dimitrakopoulos y Lou, 1994).

El objetivo de este trabajo consiste en presentar una metodología para estimar la carga hidráulica desde el punto de vista espacio-temporal, utilizando un modelo espacio-temporal tipo generalizado producto suma (De Iaco et al., 2001; De Cesare et al., 2001). En la metodología se enfatiza la determinación de la deriva; el cálculo del semivariograma muestral; el modelo espacio-temporal y la forma de ajustarlo, mediante un proceso visual, y por último, en la validación del modelo a través del método cruzado (validación cruzada). La metodología fue aplicada al acuífero de Querétaro-Obajuelo y los resultados se comparan, en términos de la validación cruzada y los valores de error, con los obtenidos al utilizar sólo el método de krigeado y multivariado (Mendoza y Herrera, 2007).

La propuesta difiere de la de Rouhani y Hall en cuanto al tipo de modelo utilizado: producto-suma; pero en la estimación se usa el mismo sistema de krigeado (krigeado ordinario espacio-temporal).

Área de estudio

La zona de estudio comprende una superficie de 426 km², que incluye el acuífero del valle de Querétaro, ubicado en la porción sur oriental de este estado, y el acuífero del valle de Obrajuelo, localizado en el estado de Guanajuato y colindante con el acuífero anteriormente citado (figura 1).

En el área se encuentran dos corrientes poco caudalosas: los ríos Querétaro y El Pueblito. La corriente principal es el río Querétaro, que corre desde La Cañada hasta interceptarse con el río El Pueblito. La información geológica y estructural indica que los valles de Querétaro y Obrajuelo se alojan en una secuencia de fosas tectónicas (originada por el tectonismo ocurrido en el Terciario), por lo que el acuífero está limitado por un sistema de fallas sepultadas por sedimentos aluviales, flujos de lava y depósitos volcanoclásticos lacustres potentes. La disposición natural de las unidades litológicas agrupadas en los medios poroso y fracturado conforman el sistema que define el flujo del agua subterránea en la zona de estudio, lo que constituye un acuífero del tipo libre, comportándose como confinado localmente en algunas zonas por la presencia de arcilla (Guysa, 1996).

Del modelo conceptual se infiere que las principales zonas de recarga al acuífero se presentan al este, por La Cañada; al noreste, por las montañas; al norte, en la zona de Jurica; al sur se tiene una entrada en el lineamiento del río El Pueblito.

Se consideran dos sistemas de flujo: uno local, formado por zonas con elevaciones del nivel freático de poca profundidad en la parte central del valle, y un sistema de flujo intermedio que sustenta casi el total de extracción en el valle, formado también por medio poroso y fracturado.

Por otro lado, en la zona límite entre los estados de Querétaro y Guanajuato hay presencia de abatimientos importantes, además de que es imposible definir la presencia de un parteaguas subterráneo o alguna frontera física que permita separar hidrogeológicamente ambos acuíferos (Simuta, 2005). La información piezométrica analizada en el presente artículo comprende un periodo de 24 años (1981-2004), que se detallará en la sección de resultados.

Marco teórico

En este artículo, la carga hidráulica se representará como una función aleatoria dependiente del espacio y el tiempo. Así, se considera un dominio D⊆R² en el espacio y un dominio T⊆R¹ en el tiempo. Una variable aleatoria espacio-temporal (VA) Z(x,t) es una variable que puede tomar una serie de valores (realizaciones) en cualquier localización en el espacio (x∈D, en un instante de tiempo t∈T. La VA Z(x,t) quedará totalmente caracterizada al conocer su función de distribución, que da la probabilidad de la variable Z en la posición x en el espacio y en el instante t:

Una función aleatoria espacio-temporal (FA) se define como una familia de variables aleatorias {Z(x,t),(x,t)∈ DxT}. Una realización de la FA Z(x,t) es una colección de realizaciones de las VA.

Kyriakidis y Journel (1999) han definido tres puntos de vista para la modelación espacio-temporal.

Función aleatoria espacio-tiempo

Típicamente descompuesta en una componente de deriva que modela una variabilidad suave de la variable aleatoria Z(x,t) y una componente estocástica, la cual denota altas frecuencias alrededor de la deriva, ambas, en espacio y tiempo (Bilonick, 1983; Bilonick, 1988; Bilonick y Nichols, 1983; Bogaert y Christakos, 1997; De Iaco et al., 2001).

Multivariada de la variable aleatoria

En este marco se conceptualiza la variable aleatoria Z(x,t) como una colección finita de observaciones Z(x), correlacionadas temporalmente (Solow y Gorelick, 1986; Comegna y Basile, 1994; Papritz y Flühler, 1994; Mendoza y Herrera, 2007).

Multivariada de serie de tiempo (S/T)

Modelos de series de tiempo, S/T, son conceptualizados como una colección finita Z(t) correlacionadas espacialmente (Rouhani y Wackernagel, 1990; Rouhani et al, 1992, Knotters y Bierkens., 2001; Kyriakidis y Journel, 2001a, 2001b).

La metodología aquí propuesta corresponde al primer punto de vista, donde el proceso espacio-temporal se descompone en:

en donde m(x,t) es una componente determinista conocida como "tendencia". La variable aleatoria R(x,t), a la que llamaremos "residuo", cumple con lo siguiente:

a) Su esperanza matemática existe; no depende del punto de observación y es nula.

b) Para todos los incrementos tiene una varianza finita.

Como en el caso espacial, estas suposiciones se conocen como "hipótesis intrínseca" y es suficiente saber lo dos primeros momentos de la variable para caracterizarla.

Sea Z(x,t) una VA Z(·) en la posición x=(x,y), al tiempo t.

Entonces, el semivariograma muestral para los incrementos Δx=(Δx, Δy) y Δt se calcula con base en la siguiente expresión:

en donde son el número de pares separados por los incrementos Δx y Δt; por definición γ(Δx=0, Δt=0)=0.

Metodología y resultados

En esta sección se muestran los resultados. Por ello, consideramos conveniente hacer una descripción de la variable analizada desde los puntos de vista espacial, temporal y espacio-temporal. Se presentan sus estadígrafos con intención de que, al momento de modelar el semivariograma espacio-temporal, se tenga una idea de la forma en cómo varía en cada uno.

Estadígrafos en espacio

Como primer paso, se analizó en forma separada la información espacial de la temporal. De esta manera, se describieron los estadígrafos que caracterizan la carga hidráulica para cada año, utilizando las librerías desarrolladas por Deutsch y Journel (1998).

En el diagrama de barras de la figura 2 se muestra el número de lecturas anuales de pozos de observación. Cabe destacar que se cuenta con información rica en el espacio, pero en algunos tiempos, con pocas observaciones. Para el periodo 1970-2004 existen años que van desde cero lecturas (2000) hasta un máximo de 170 (2001).

La base de datos en ocasiones es rica en espacio, pero en otras muy pobre, lo que afecta directamente al cálculo de la estructura del semivariograma. Por tal motivo se procede a trabajar con aquellos años que presentan treinta o más observaciones, mismas que se muestran en la figura 3a. Se observa cómo la muestra se redujo a veinte columnas representantes de los años de muestreo, que cubren desde 1981 hasta 2004, teniendo como número mínimo de lecturas treinta, para 1997, y máximo de 173, para 2001.

Los estadígrafos que caracterizan las lecturas de la carga hidráulica se resumen en el cuadro 1. Los estadígrafos indican que en todos los años los valores de la media y mediana están separados, y se reportan valores de desviaciones estándar altas con valores de simetría que señalan sesgo en la muestra y forma platicúrtica.

Por lo antes comentado, será necesario transformar los datos, con objeto de trabajar con muestras cuyos valores asemejen una distribución simétrica.

A diferencia del análisis exploratorio espacial, el temporal toma en cuenta la variación de la variable aleatoria, carga hidráulica, en una misma posición; es decir, analiza la serie de tiempo de la carga hidráulica que se presenta en un pozo de observación.

Observando la base de datos, es evidente que la información temporal no se encuentra densamente muestreada; basta considerar que de los 322 pozos existentes, para uno hay un máximo de 31 lecturas y, como mínimo, una. La figura 4 muestra el histograma del número de observaciones en el periodo 1970-2004.

Por dicha razón, se decidió trabajar con pozos que reportaran, en el tiempo, 19 o más lecturas; así se trabajó con una submuestra que contabiliza un total de 13 pozos, tal como se observa en la figura 3b. Basta comentar que la decisión anterior obedece a que al calcular el semivariograma muestral en tiempo, se requiere de un número mínimo de observaciones para evitar que el semivariograma tenga un comportamiento errático.

El cuadro 2 presenta los estadígrafos de los trece pozos. Hacemos notar que la muestra no guarda un comportamiento de una distribución simétrica.

Estimación de la deriva

Es común que los niveles piezométricos muestren una pendiente global en la dirección del flujo, por lo que suponemos que la media de la variable sí depende de su posición en el espacio y el tiempo; es decir, existe una deriva. Se obtuvieron los semivariogramas muestrales para los niveles piezométricos para los años con más datos y se observó que éstos crecen más rápidamente que h². Se sabe que una función aleatoria intrínseca es aquella cuyos incrementos tienen esperanza matemática y varianza definidas e independientes de la posición para todo vector h (Samper y Carrera, 1990), por lo que en este caso la variable es intrínseca y es lógico pensar que existe una derivada.

Para quitar el efecto de deriva se procede a tratarla como indica la ecuación (3) y expresando la deriva como:

donde los valores de los coeficientes son:

Una vez determinada la deriva, ésta se resta de la observación inicial y se trabaja con los residuos. Suponemos que los residuos son aleatorios y se comportan como funciones aleatorias intrínsecas. Los estadígrafos de los residuos, en el espacio, se muestran en el cuadro 3 y, en el tiempo, en el cuadro 4.

En resumen, desde el punto de vista espacial, se tiene que la variable carga hidráulica presenta un comportamiento ascendente en la distancia, típico en acuíferos, que refleja un efecto de deriva. Se considera que los residuos son estacionarios de segundo orden (es decir, existe la media y no depende de la posición, y la covarianza existe y depende únicamente del vector de separación entre las posiciones), de esta forma se trabajará con ellos para determinar su estructura y, posteriormente, calcular los parámetros que los caracterizan.

Análisis estructural

El objetivo de esta sección consiste en determinar la estructura de los datos en espacio-tiempo, por lo que se procederá a calcular el semivariograma muestral de los años identificados en el cuadro 3. Se calculó el semivariograma espacio-temporal y luego se efectuó la validación cruzada del modelo.

Las propuestas de Bilonik (1983), Cameron y Hunter (2006), y Gardner y Sullivan (2004) se basan en calcular un semivariograma promedio, partiendo de cada semivariograma individual, en una distancia "lag" común. Este tipo de semivariograma es denominado "compuesto". En este trabajo se propone analizar un semivariograma que de forma conjunta integre la información espacial y temporal, tal como se indica en la siguiente sección.

Semivariograma espacio-temporal

Se procedió a calcular el semivariograma muestral (ecuación (4)), utilizando el algoritmo propuesto por De Cesare et al. (2002). El semivariograma espacio-temporal se presenta en la figura 5.

De Cesare et al. (2001) proponen el modelo que combina el producto y suma de la siguiente manera:

o, de forma equivalente, el semivariograma:

donde C_s y C_t son la función de covarianza en espacio y tiempo, respectivamente; γ_s y γ_t son las correspondientes funciones de semivariograma, en espacio y tiempo. Es de remarcar que C_st(0) es la cima de γ_st; C_s(0) es la cima de γ_s, y C_t(0) es la cima de γ_t, y por definición γ_st(0,0)=γ_s(0)=γ_t(0)=0. Por otro lado, se deben cumplir las condiciones siguientes: k₂+k₁C_t(0)=1 y k₃+k₁C_s(0)=1. Para ajustar el modelo teórico al semivariograma muestral, se requiere estimar la cima global, espacio-tiempo, misma que se puede definir a través de un procedimiento de ajuste visual. En este estudio se evaluaron los modelos teóricos (esférico, exponencial, Gaussiano, potencia y agujero) presentes en las utilerías del GSLIB.

El procedimiento para el ajuste del semivariograma espacio-temporal se puede hacer valiéndose de los semivariogramas teóricos en espacio y tiempo, y la selección de la cima global por medio de un ajuste visual que reproduzca el semivariograma muestral. Este proceso depende del analista, por lo que el valor obtenido varía dependiendo de quién lo proponga. Este proceso no es ajeno a la metodología que aplican los geoestadísticos en el caso de análisis univariados.

Como primer paso, se graficó el semivariograma muestral para diferentes incrementos de distancia en el eje del espacio y en el eje para el tiempo. Efectuado lo anterior, se propuso un semivariograma, producto suma, donde se conoce la cima de la componente espacial y temporal.

Para ello se propusieron diferentes valores de los parámetros presentes en la ecuación (8), hasta ajustar el semivariograma teórico al muestral, buscando que la media de las diferencias, que denominaremos "errores", entre el valor del semivariograma teórico y el muestral, fuera mínima. Siguiendo los pasos anteriores, se llegó al semivariograma teórico mostrado en la figura 6: el valor de error es de 12.59. Este proceso es de prueba y error hasta determinar el mejor ajuste; en el artículo no se mencionan los diferentes valores de error ajustados para los diversos modelos teóricos.

Al cuantificar el error de las mediciones entre la función de semivariograma y el respectivo modelo teórico, se obtiene una forma cuantitativa para seleccionar el mejor modelo que, en este caso, da como resultado un valor de la cima, C_rf(0,0), igual a 1 840.

Es importante remarcar que este criterio no es automático, ya que en el proceso de ajuste del semivariograma espacio-tiempo se van proponiendo diferentes valores hasta llegar a un valor que, de acuerdo con los errores, sea el menor. Sin embargo, dicho proceso permite al modelador ajustar de acuerdo con su criterio de mejor manera en las interdistancias de espacio o tiempo que muestren un mayor número de pares (ya que los valores de éstas en general debieran ser más representativas del semivariograma), y desechando las que no sean congruentes; esto se ve reflejado en la figura 6. El gráfico muestra un mejor ajuste en las interdistancias (h) más cercanas. Conforme se incrementa la interdistancia, el semivariograma muestral presenta mayor dispersión y el semivariograma teórico se ajusta en la porción central. Verificando el número de pares para los puntos con interdistancias mayores, se confirma que su valor es bajo. Este tipo de criterios son importantes para que el modelador, bajo su criterio y experiencia, decida el mejor ajuste. Esto no pasa con un proceso de ajuste automático, sea por mínimos cuadrados o mínimos cuadrados ponderados.

Es importante destacar que los resultados se ven afectados al momento de decidir el valor de la cima del modelo espacio-tiempo. Isaaks y Srivastava (1989) ejemplifican la sensibilidad de los datos al seleccionar un modelo, tratando el efecto que tienen cuando se varía la escala (diferente cima), la forma (diferente modelos), pepita (diferente pepita), el rango (distintos rangos) y la anisotropía (distintas orientaciones). Cuando se varía la cima de los modelos, no existe un efecto en los valores estimados y los pesos del krigeado. Sin embargo, sí existe una afectación en la varianza, modificando el valor de la varianza en la misma proporción (factor de escalamiento) en que se varía la cima del modelo propuesto.

Los ejemplos de Issaks y Srivastava orientan de forma intuitiva al modelador acerca de la sensibilidad en la estimación de la carga hidráulica y el valor de la varianza, al variar los parámetros del modelo.

Por esta razón es necesario llevar a cabo un análisis comparativo. En este trabajo se hizo a partir del análisis de los errores de la validación cruzada, lo cual se comenta en la siguiente sección.

Validación cruzada

Para verificar qué tan bien ajusta el semivariograma teórico, se procede a utilizar el método de validación cruzada para después calcular el error medio (EM), el error cuadrático medio (ECM) y error cuadrático medio estándar (ECME) (ver la expresión matemática de cada uno de ellos en el cuadro 5).

El EM mide la sobreestimación o subestimación de la variable; si es positivo, la variable en promedio se subestima y si es negativo lo contrario; la estimación es mejor si éste es cercano a cero. El ECM da una medida promedio de la magnitud al cuadrado de los errores; también se espera que el valor sea pequeño. Con respecto al ECME, los valores se interpretan como el ajuste que presenta la varianza del modelo del semivariograma, espacial y espacio-temporal, con los errores de estimación; por esta razón, cuando el valor del ECME es cercano a la unidad, se interpreta como que la varianza del modelo teórico y el error cuadrático de la estimación en promedio son cercanos.

Los resultados de validación cruzada para 1993 muestran que los valores de error son los siguientes: EM=-1.045, ECM=195.53, ECME= 0.408; el valor promedio de desviación estándar es de 31.60. La figura 7a muestra el diagrama de dispersión entre los valores observados de la carga hidráulica para 1993 y los estimados por el método de validación cruzada.

Los resultados de validación cruzada para 1995 muestran que los valores de error son los siguientes: EM=-0.131, ECM=137.64, ECME=0.361; el valor promedio de desviación estándar es de 18.0. La figura 7b muestra el diagrama de dispersión entre los valores observados de la carga hidráulica para 1995 y los estimados por el método de validación cruzada.

Por último, los resultados de validación cruzada para 1999 muestran siguientes valores de error: EM=-0.967, ECM=339.706, ECME=0.683; el valor promedio de desviación estándar es de 29.61. La figura 7c muestra el diagrama de dispersión entre los valores observados de la carga hidráulica para 1999 y los estimados por el método de validación cruzada.

Estimación espacio-temporal

Los resultados de la estimación de la carga hidráulica en el acuífero de Querétaro-Obrajuelo se efectuaron en una malla de 41 columnas y 51 renglones, con una distancia entre renglón y columna de 1 100 m y un delta de tiempo de una unidad. Sólo se reportan los planos de la carga hidráulica de los años 1993, 1995 y 1999, dejando de lado los años intermedios y la predicción correspondiente a los años posteriores a 2003, último año que cuenta con lecturas.

Las ilustraciones 8 a 13 muestran los mapas de la carga hidráulica y la varianza para los años comentados, y de forma global se describen los siguientes patrones: en la figura 8 se muestra la estimación de la carga hidráulica para el año 1993; se observa que las líneas equipotenciales son acordes con las mediciones reportadas para cada pozo de observación. Cabe destacar que se marca un parteaguas en la zona de la laguna El Salitre y en el área del acuífero de Querétaro se tiene una mejor definición de las zonas con mayor abatimiento. Esta definición se observa en el mapa como áreas concéntricas de una influencia local que difieren del entorno. Al verificar estos puntos corroboramos que, efectivamente, corresponden a pozos donde existe una carga hidráulica menor al entorno.

En el caso de La Cañada y zonas altas de la cuenca, se observa el valor subestimado de la carga hidráulica; sin embargo, de forma general, se mantiene el mismo patrón de las líneas de flujo en dirección al centro del valle de Querétaro.

El plano de varianza (figura 9) presenta que la porción central del acuífero está mayormente muestreada; la línea equipotencial de 400 unidades se distribuyen en un amplio sector y, las márgenes del acuífero pierden esta densidad de muestreo.

Respecto a la estimación proyectada en el plano para 1995 (figura 10), se mantiene que las estimaciones son acordes con las mediciones reportadas para cada pozo de observación. El parteaguas en la zona de la laguna El Salitre se mantiene, así como la definición de las zonas con mayor abatimiento. En el caso de La Cañada y zonas altas de la cuenca, se observa que el valor de la carga hidráulica es subestimado; sin embargo, de forma general se continúa el mismo patrón de las líneas de flujo en dirección al centro del valle de Querétaro.

La estimación proyectada en el plano para 1999 (figura 12) presenta una correspondencia entre lo estimado y las mediciones reportadas para cada pozo. Cabe señalar que para este año no se marca un parteaguas en la zona de la laguna El Salitre, igurando que los acuíferos Los Apaseos, Obrajuelo y Querétaro han sido afectados en su dinámica original por su inversión en el flujo, sobre todo en Los Apaseos, ahora centro del acuífero del valle de Querétaro. Esta hipótesis deberá ser corroborada con información adicional obtenida para esta zona. Este parteaguas es acorde con el modelo conceptual comentado en el área de estudio para el caso no estacionario, pero la posición del mismo no está bien definida y mucho menos confirmada.

Que el valor del error cuadrático medio estándar sea elevado; de hecho, el más elevado de los tres años analizados, indica de forma tácita que el modelo del semivariograma espacio-tiempo no es bueno para estimar este año en particular. Es de esperar que en los puntos donde se estima, dentro de la malla indicada, no sean siempre acordes con la realidad.

El plano de varianza (figura 13) presenta un aumento de la varianza en la zona oeste, producto de una reducción de los pozos de muestreo.

Estimación en tiempo

En este punto se muestran los resultados obtenidos con el método de validación cruzada en los pozos que cuentan con 19 o más observaciones de la carga hidráulica, mismos que se listan en el cuadro 4. La intención es mostrar las lecturas de estimación de la carga hidráulica, comparadas con los valores observados, cuando se ha quitado el dato que se estima. De esta forma, sólo se mostrarán los resultados para luego ser discutidos. La figura 14 muestra tres hidrógrafos de los pozos 759, 752 y 629.

La posición espacial de los pozos seleccionados se muestra en la figura 15. Es fácil identificar que los aprovechamientos no están en una zona restringida; por el contrario, se ubican en todo el acuífero. De forma particular, el pozo 752 se sitúa en la zona de recarga; el pozo 629, en la porción central del acuífero Querétaro y el pozo 759, en la parte alta del acuífero. La falta de información histórica en la porción del acuífero Obrajuelo se da porque el mínimo de observaciones, 19 lecturas desde el año 1981 hasta el 2004, no se cumplió. Esta imposición obedece a que es necesario contar con un número mínimo de observaciones para generar los semivariogramas en el tiempo.

Aclarado el punto anterior, se comentan los resultados de los hidrógrafos. El pozo 759 se ubica en la porción alta del acuífero y de manera general describe un descenso progresivo de la carga hidráulica. Para el año de 1981 se reporta una elevación de la carga hidráulica de 1 910.92 m y, para el año del 2003, se reportan 1 879.14 m; la diferencia en estos 23 años es de 31.78 metros.

El descenso de la carga hidráulica es constante, de forma gradual, existiendo algunos puntos de recuperación, como se distingue en el año de 1990, donde el nivel se recupera 3.0 m, tomando en cuenta la tendencia general, para después descender de manera gradual. Cabe resaltar que en el tiempo, entre los años 1999 y 2001, se mantiene constante, para después descender 6.51 metros.

Respecto a los valores estimados, se observa que de manera general se ajustan a la tendencia, presentando las mayores diferencias en las series donde no existen observaciones en el año inmediato anterior y posterior (ver periodo 1981-1990), pero en los periodos donde hay observaciones, las estimaciones presentan un mejor ajuste; de hecho, las diferencias entre ambas no sobrepasan los 1.5 metros.

Respecto a los valores estimados, es evidente que la estimación de la carga hidráulica no se ajusta a las mediciones en los pozos que se encuentran en las márgenes del acuífero. En estos sitios se identifica que las estimaciones están por debajo de las mediciones; es decir, se subestiman. La falta de información, en el tiempo, ocasiona que los puntos de estimación sean extrapolados. Los valores de varianza indican que las estimaciones presentan una incertidumbre alta.

Los valores estimados, que presentan una mayor diferencia entre los valores medidos, se identifican en los intervalos de tiempo donde no existen datos antes o después. Lo contrario sucede en los intervalos de tiempo donde hay información en tiempos cortos.

Journel y Rossi (1989) discuten el efecto que tienen las estimaciones, en los casos dónde se extrapola; aunque es aplicado en el caso espacial, el efecto es el mismo en el tiempo.

El pozo 629 se ubica en la porción central del acuífero Querétaro. De forma general describe un descenso progresivo de la carga hidráulica. Para 1981 se reporta una elevación de la carga hidráulica de 1 753.12 m, y para el 2003 se reportan 1 677.69 m. La diferencia en estos 23 años es de 75.43 m. El descenso de la carga hidráulica es gradual. Respecto a los valores estimados, se observa que de manera global se ajustan a la tendencia, que es lo buscado cuando no se cuenta con observaciones.

Discusión

Recapitulando, tenemos que en el aspecto espacial se seleccionaron 23 años (cuadro 1), los cuales cumplen con treinta o más observaciones (pozos) en cada muestreo. Este criterio se contempla considerando que al determinar la estructura de la variable se requiere de un número mínimo para efectuarlo. En la literatura, para calcular el semivariograma, se recomienda utilizar un número mayor de puntos de observación (cien observaciones); sin embargo, en la práctica, se trabaja con lo que se tiene, siempre manteniendo en mente que las decisiones que se toman tienen un costo. En el aspecto temporal, se seleccionaron los pozos que presentan 19 o más observaciones, dando como resultado 13 pozos.

Así, se procedió a efectuar el análisis exploratorio de los datos, el análisis estructural espacio-tiempo y la validación del modelo propuesto, por lo que los comentarios de la discusión se enfocan hacia la validación y la estimación.

Es importante destacar que la validación cruzada se llevó a cabo tomando en cuenta el modelo espacio-temporal, por lo que los resultados comprenden del año 1981 al 2004. De este periodo se seleccionaron los años 1993, 1995 y 1999 para mostrar el grado de ajuste que tienen las estimaciones con las observaciones. El lector debe considerar que la estimación con este método no sólo muestra resultados en el espacio, sino también en el tiempo.

La figura 7a muestra el diagrama de dispersión entre los valores medidos de la carga hidráulica y los estimados para 1993. Los valores de la carga hidráulica medidos versus los estimados presentan un índice de regresión de R² = 0.94, que indica que éstos presentan una buena correlación.

El cuadro 5 presenta un resumen de los errores de la validación cruzada obtenidos con los métodos kriging ordinario, cokriging y estimación espacio-tiempo. Una descripción detallada del procedimiento utilizado en los dos primeros métodos se puede encontrar en Mendoza y Herrera (2007). Para el caso espacio-tiempo de las variables aleatorias se utiliza la anotación ET para distinguirla de la técnica multivariada (Cok), también es importante explicar que cuando en el cuadro el método ET refiere a un año específico, se reportan los resultados de este método para el tiempo señalado (1993, 1995, 1999); es decir, la técnica estima en el espacio y también en el tiempo, utilizando los datos de todo el periodo, pero los estadígrafos mostrados en el cuadro se calcularon usando únicamente los datos para el año especificado. En contraste, lo que se reporta en la cuarta línea, con la anotación ET 19981–2004, son los estadígrafos de la validación cruzada manejando toda la información; o sea, tomando en cuenta los errores de la validación de la carga hidráulica para el periodo 1981–2004.

Aclarado lo anterior y retomando los resultados que muestra el cuadro, se identifica que los valores estimados de los residuos con el método de validación cruzada, para todos los casos, presentan un comportamiento simétrico.

Analizando los valores del error (EM) para el caso de KO, se encuentra que son reducidos, cercanos a cero y todos ellos negativos, lo que indica una subestimación de los valores de la carga hidráulica.

En el caso del Cok, la magnitud del EM es mayor que la del KO y, en todos los casos, es negativo, de la misma forma que para el KO.

Para el ET, sólo en el caso donde se consideran todas las observaciones, el valor del EM es bajo, pero para cada año en particular, el error se incrementa un poco.

En términos generales, el EM indica que para cada técnica, las estimaciones están siendo subestimadas.

Para el caso del ECM, pensando éste como un valor de dispersión, ya que es el cuadrado de los errores, se observa que tanto en KO y Cok los valores son similares; sin embargo, al compararlos con el ET, el ECM se reduce drásticamente para los años 1993 y 1995; en 1999 se mantiene. El promedio del error cuadrático medio para los tres años en el ET es de 224.29 m, mientras que en el KO es de 356.4 y en el Cok de 314.27. Esta reducción en el promedio del ECM para el ET se debe a que la estimación se realiza tomando en cuenta un número mayor de puntos de observación, no sólo en el espacio, sino también en el tiempo.

De esta forma, el ECME indica que el mejor ajuste entre el error cuadrático y la varianza lo ofrece el método de Cok, seguido del KO, el ET para el periodo 1981-2004 y, por último, el ET para cada año. El hecho de que el ECME del ET para todo el periodo tenga un valor de 0.804, mucho más cercano que los valores para cada año por separado, indica que la varianza del modelo teórico representa de mejor forma el error cuadrático para todo el periodo que para los años por separado.

Respecto a la estimación espacio-temporal proyectada en el plano, dejando al tiempo invariante, se observa que las estimaciones son acordes con las mediciones reportadas para cada pozo de observación. Cabe destacar que el método distingue un parteaguas en la zona de la laguna El Salitre y en el área del acuífero de Querétaro, y define las zonas con mayor abatimiento. Esta definición se observa en los planos como áreas concéntricas de una influencia local que difieren del entorno. Al verificar estos puntos corroboramos que, efectivamente, corresponden a pozos donde existe una carga hidráulica menor al entorno.

En el caso de La Cañada y zonas altas de la cuenca, se observa el valor subestimado de la carga hidráulica; sin embargo, de forma general, se mantiene el mismo patrón de las líneas de flujo en dirección al centro del valle de Querétaro.

Ahora bien, verificando los puntos que de acuerdo con los diagramas de dispersión (figura 7a) se alejan de la línea de 45 grados, correspondientes al caso donde la estimación es igual a la medición, tenemos que existen cinco valores: 1 721.94; 1 734.17; 1 766.72; 1 866.80 y 1 890.94, que se ubican en las márgenes del acuífero (figura 8), caracterizados por pocos datos al hacer la estimación. Por la posición, se identifica que corresponden a puntos extrapolados. Verificando el valor de varianza que se obtiene para cada uno de estos puntos (figura 9) es posible corroborar que la varianza es alta.

En el caso de la piezometría para 1995, se tiene que los valores de la carga hidráulica medidos versus los estimados, presentan un índice de regresión de R²=0.95 (figura 7b), considerándose como buena correlación.

Verificando los puntos que se resaltan en el diagrama de dispersión (figura 7b), tenemos que existen siete valores: 1 629.24; 1 658.79; 1 715.60; 1 734.20; 1 866.98; 1 874.85 y 1 878.47, que al colocarse en el plano se ubican en las márgenes del acuífero (figura 10) y se caracterizan por tener pocos datos cercanos. Por la posición, al realizar la estimación, se identifica que no corresponde a un punto de interpolación; por el contrario, es un punto extrapolado.

Verificando el valor de varianza que se obtiene para cada punto (figura 11) es posible corroborar que la varianza es alta.

Al verificar los puntos que resaltan en el diagrama de dispersión (figura 7c) existen siete valores: 1 623.33; 1 650.5; 1 669.5; 1 702.1; 1 720.05; 1 792.06 y 1 864.61, que al colocarse en el plano se ubican en las márgenes del acuífero (figura 12), caracterizados por tener pocos datos para la estimación. Por la posición, al hacer la estimación, se identifica que no corresponde a un punto de interpolación; por el contrario, es un punto extrapolado.

Verificando el valor de varianza que se obtiene para cada punto (figura 13) es posible corroborar que la varianza es alta.

Con respecto a la varianza, es necesario resaltar que para todos los casos los planos son menores que los obtenidos con el método de krigeado. Esto es congruente con la teoría, ya que al aumentar el número de observaciones se reduce la varianza de la estimación.

Conclusiones

La estimación espacio-tiempo es una herramienta útil que puede ser tomada como alternativa para efectuar estimaciones de la carga hidráulica en un contexto espacial y temporal.

La aplicación de esta técnica para determinar la carga hidráulica es innovadora, ya que resuelve el problema de los modelos lineales. La metodología desarrollada ofrece una alternativa para tratar las variables que presentan una deriva espacio-tiempo, así como indicar la forma de ajustar el modelo producto suma, propuesto por De Iaco et al. (op.cit.).

Para realizar la verificación del modelo, se utilizó el método de validación cruzada que permite evaluar la certidumbre del modelo propuesto. Los valores de EM, ECM y ECME indican que el modelo producto suma se ajusta de manera adecuada a la información espacio-tiempo; sin embargo, al evaluar de forma individual los resultados de validación cruzada para los años seleccionados, se reportan algunos puntos con errores grandes.

Las estimaciones de la carga hidráulica para los años 1993, 1995 y 1999 son congruentes con los valores observados en campo, así como la definición del parteaguas entre los acuíferos.

Una ventaja del método espacio-tiempo es el incremento de observaciones para realizar la estimación; esta característica permite hacer estimaciones en zonas donde no se cuenta con información.

El método es poderoso, ya que también realiza estimaciones en diversos tiempos. Esto se puede apreciar en la figura 14, donde se observa que el método permite efectuar interpolaciones entre años donde se carece de información.

El método espacio-tiempo requiere de un análisis estructural adecuado, ya que a partir de éste se seleccionan los parámetros que lo caracterizan. En nuestro caso se aplica una metodología en donde se contempla la selección de los parámetros, tomando en cuenta diferentes secciones de tiempo y espacio, y se refuerzan con el índice de error. Esta propuesta es diferente a la de De Cesare et al. (op. cit.), ya que de manera gráfica se determina la cima del semivariograma de superficie. Así, ellos proponen un modelo producto suma que garantiza que el modelo es positivo definido y que no afecta en el sistema de krigeado. Estas dos innovaciones superan la propuesta inicial de Rouhani y Hall (1989).

Agradecimientos

A la Comisión Nacional del Agua, que otorgó el apoyo para realizar los estudios de doctorado del primer autor.

Referencias

AHMED, S. and DE MARSILY, G. Comparison of geostatistical methods for estimating transmissivity using data on transmissivity and specific capacity. Water Resources Research. Vol. 23, no. 9, 1987, pp. 1717-1737.         [ Links ]

ANGELINI, P. Correlation and spectral analysis of two hydrogeological systems in central Italy. Hydrological Sciences Journal. Vol. 42, no. 3, 1997, pp. 425-438.         [ Links ]

ASCE. Review of geostatistics in geohydrology I: Basic concepts. Journal of Hydraulic Engineering. Vol. 116, no. 5, 1992a, pp. 612-632.         [ Links ]

ASCE. Review of geostatistics in geohydrology II: Aplications. Journal of Hydraulic Engineering. Vol. 116, no. 5, 1992b, pp. 633-658.         [ Links ]

ASLI, M. and MARCOTTE, D. Comparison of approaches to spatial estimation in a bivariate context. Mathematical Geology. Vol. 27, no. 5, 1995, pp. 641-658.         [ Links ]

BILONICK, R.A. and NICHOLS, D.G. Temporal variations in acid precipitation over New York State What the 1965–1979 USGS data reveal. Atmospheric Environment. Vol. 17, no. 6, 1983, pp. 1063-1072.         [ Links ]

BILONICK, R.A. Risk qualifies maps of hydrogen ion concentration for the New York State area for 1966-1978. Atmospheric Environment. Vol. 17, no. 12, 1983, pp. 2513-2524.         [ Links ]

BILONICK, R.A. Monthly hydrogen ion depositation maps for the northeastern U.S. from July 1982 to September 1984. Atmospheric Environment. Vol. 22, no. 9, 1988, pp. 1909-1924.         [ Links ]

BOGAERT, P. and CHRISTAKOS, G., Sapatiotemporal analysis and processing of thermometric data over Belgium. Journal of Geophysical Research. Vol. 102, no. 22, 25, 1997, pp. 831-846.         [ Links ]

CAMERON, K. and HUNTER, P. Optimization of LTM networks using GTS: statistical approaches to spatial and temporal redundancy, Air Force Center for Environmental Excellence, Brooks AFB, TX, [en línea] Disponible para World Wide Web: http://www.afcee.brooks.af.mil/products/rpol/docs/GTSOptPaper.pdf, acceso libre [consulta, agosto de 2006]         [ Links ].

CHRISTAKOS, G. Modern Spatiotemporal Geostatistics. New York: Oxford University Press, 2000, 288 pp.         [ Links ]

COMEGNA, V. and BASILE, A. Temporal stability of spatial patterns of soil water storage in a cultivated Vesuvian soil. Geoderma. Vol. 62, 1994, pp. 299-310.         [ Links ]

DE CESARE, L., MYERS, D.E. and POSA, D. Estimating and modeling space-time correlation structures. Statistics & Probability Letteres. Vol. 51, 2001, pp. 9-14.         [ Links ]

DE CESARE, L., MYERS, D.E. and POSA, D. FORTRAN Programs for space-time modeling. Computers Geosciences. Vol. 28, 2002, pp. 205-212.         [ Links ]

DE IACO, S., MYERS D.E. and POSA, D. Space-time analysis using a general product-sum model. Statistics & Probability Letters. Vol. 52, 2001, pp. 21-28.         [ Links ]

DELHOMME, J.P. Kriging in the hydrosciences. Advances in Water Resources. Vol. 1, no. 5, 1978, pp. 251-266.         [ Links ]

DEUTSCH, C.V. and JOURNEL, A.G. GSLIB Geostatistical software library and User's guide. Applied geostatistics series. Second edition. New York: Oxford University Press, 1998, 369 pp.         [ Links ]

DIMITRAKOPOULOS, R. and LUO, X. Spatiotemporal modeling: covariances and ordinary kriging systems. Dimitrakopoulos, R. (editor). Geostatistics for the next century. Dordrecht: Kluwer Academic Publ., 1994, pp. 88-93.         [ Links ]

GAMBOLATI, G. and GALEATI, G. Comment on analysis of non intrinsic spatial variability by residual kriging whit application to regional groundwater levels by Neuman, S.P. and Jacobson, E.A. Mathematical Geology. Vol. 19, no. 3, 1987, pp. 249-257.         [ Links ]

GAMBOLATI, G. and VOLPI, G. Groundwater contour mapping in Venice by stochastic interpolators. Water Resources Research. Vol. 15, no. 2, 1979a, pp. 281-290.         [ Links ]

GAMBOLATI, G. and VOLPI, G. A conceptual deterministic analysis of the kriging technique in hidrology. Water Resources Research. Vol. 15, no. 3, 1979b, pp. 625-629.         [ Links ]

GARDNER, B. and SULLIVAN, P. Spatial and temporal stream temperature prediction: modeling nonstationary temporal covariance structures. Water Resources Research. Vol. 40, 2004, W01102.         [ Links ]

GUYSA. Estudio de simulación hidrodinámica y diseño óptimo de la red de observación en los acuíferos de Aguascalientes y Querétaro. Tomo III. Querétaro: Comisión Nacional del Agua, Gerencia de Aguas Subterráneas, Geofísica de Exploraciones Guysa, S.A. de C.V. GAS-008-PR-96, 1996.         [ Links ]

HEVESI, J.A., ISTOK, J.D. and FLINT A.L. Precipitation estimation in mountainous terrain using multivariate geostatistics. Part I: Structural analysis. Journal of Applied Meteorology. Vol. 31, 1992a, pp. 661-676.         [ Links ]

HEVESI, J.A., ISTOK, J.D. and FLINT, A.L. Precipitation estimation in mountainous terrain using multivariate geostatistics. Part II: Isohyetal maps. Journal of Applied Meteorology. Vol. 31, 1992b, pp. 677-688.         [ Links ]

IBAÑEZ, G.M.V. Modelos estadísticos espacio-temporales en perimetría. Tesis doctoral. Castellón de la Plana, España: Departamento de Matemáticas, Escuela Superior de Tecnología y Ciencias Experimentales, Universitat Jaume, 2003, 322 pp.         [ Links ]

JOURNEL, A.G. and ROSSI, M.E. When do we need a trend model in kriging. Mathematical Geology. Vol. 21, no. 7, 1989, pp. 715-740.         [ Links ]

KNOTTERS, M. and BIERKENS, M.F.P. Predicting water table depths in space and time using a regionalized time series model. Geoderma. Vol. 103, 2001, pp. 51-77.         [ Links ]

KYRIAKIDIS, P.C. and JOURNEL, A.G. Geostatistical space-time models: A review. Mathematical Geology. Vol. 31, no. 6, 1999, pp. 651-684.         [ Links ]

KYRIAKIDIS, P.C. and JOURNEL, A.G. Stochastic modeling of atmospheric pollution: a spatial time-series framework. Part I: methodology. Atmospheric Environment. Vol. 35, 2001a, pp. 2331-2337.         [ Links ]

KYRIAKIDIS, P.C. and JOURNEL, A.G., Stochastic modeling of atmospheric pollution: a spatial time-series framework. Part II: application to monitoring monthly sulfate deposition over Europe. Atmospheric Environment. Vol. 35, 2001b, pp. 2339-2348.         [ Links ]

LAROCQUE, M., MANGIN, A., RAZACK, M. and BANTON, O. Contribution of correlation and spectral analyses to the regional study of a large karst aquifer (Charente, France). Journal of Hydrology. Vol. 205, 1998, pp. 217-231.         [ Links ]

LEE, J.Y. and LEE, K. Use of hydrologic time series data for identification of recharge mechanism in a fractured bedrock aquifer system. Journal of Hydrology. Vol. 229, 2000, pp. 190-201.         [ Links ]

MENDOZA, C.E.Y. y HERERA, G. Estimación multivariada para determinar la carga hidráulica en el acuífero Querétaro-Obrajuelo. Ingeniería hidráulica en México. Vol. XXII, núm. 1, enero-marzo de 2007, pp. 63-79.         [ Links ]

MENDOZA, C.E.Y. Análisis de alternativas para la estimación de la carga hidráulica utilizando métodos geoestadísticos en espacio y espacio-tiempo. Tesis doctoral. México, D.F.: Universidad Nacional Autónoma de México, 2008, 253 pp.         [ Links ]

MYERS, D.E. and JUORNEL, A.G. Variograms whit zonal anisotropies and noninvertible kriging sytems. Mathematical Geology. Vol. 22, no. 7, 1990, pp. 779–785.         [ Links ]

PAPRITZ, A. and FLÜHLER, H. Temporal change of spatially autocorrelated soil properties: optimal estimation by cokriging. Geoderma. Vol. 62, 1994, pp. 29–43.         [ Links ]

ROUHANI, S. Comparative study of ground water mapping techniques. Journal of Ground Water. Vol. 24, no. 2, 1986, pp. 207-216.         [ Links ]

ROUHANI, S. and HALL, T. Space-time kriging of groundwater data. Armstrong, M. (editor). Geostatistics. Vol. 2. Dordrecht: Kluwer Academic Publ., 1989, pp. 639–650.         [ Links ]

ROUHANI, S. and MYERS D.E. Problems in space-time kriging of geohydrological data. Mathematical Geology, Vol. 22, no. 5, 1990, pp. 611-623.         [ Links ]

ROUHANI, S. and WACKERNAGEL, H. Multivariate geostatistical approach to space-time data analysis. Water Resources Research. Vol. 26, no. 4, 1990, pp. 585-591.         [ Links ]

ROUHANI, S., EBRAHIMPOUR, M.R., YAQUB, I. and GIANELLA, E. Multivariate geostatistical trend detection and network evaluation of space-time acid deposition data I. Methodology. Atmospheric Environment. Vol. 14, 1992, pp. 2603-2614.         [ Links ]

SAMPER, J. y CARRERA, J. Geoestadística. Aplicaciones a la hidrología subterránea. Barcelona: CIMNE, 1990, 481 pp.         [ Links ]

SIMUTA, C.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 ]

SOLOW, A.R. and GORELICK, S.M. Estimating monthly streamflow values by cokriging. Mathematical Geology. Vol. 18, no. 8, 1986, pp. 785-809.         [ Links ]