1–D numerical modelling of the infiltration dynamics in the unsaturated zone, Toluca valley aquifer

Javier Salas–García^1,*, Jaime Gárfias¹, Hilario Llanos² y Richard Martel³

¹ Universidad Autónoma del Estado de México, Facultad de Ingeniería, Cerro de Coatepec, Ciudad Universitaria, 50130 Toluca, Estado de México, México.

² Universidad del País Vasco, Departamento de Geodinámica, Vitoria–Gasteiz, País Vasco, España.

³ Institut National de la Recherche Scientifique, Centre Eau, Terre et Environnement, 490 de la Couronne, Québec (Qc), G1K 9A9 Canadá. * proyectos@javiersalasg.com

RESUMEN

En el valle de Toluca se han producido importantes desplazamientos del terreno relacionados con la activación de fracturas de reciente generación que localmente condicionan la infiltración de aguas superficiales contaminadas. Tal es el caso de la zona de fracturas existente en el municipio de Santiago Tianguistenco (Estado de México) a la que afluye un canal de aguas residuales, principalmente domésticas, procedentes de unos 14,000 habitantes. El objetivo de la presente investigación consiste en cuantificar el proceso de infiltración en la zona no saturada de la mencionada área de fracturas. Para ello, se siguió un esquema de trabajo que combina tanto el desarrollo tecnológico como la modelización numérica, ambos sustentados en un intenso trabajo de campo. Así, los resultados de los sondeos eléctricos verticales efectuados contrastados con la información aportada por la perforación de un pozo de monitoreo realizado durante el propio proyecto, han puesto de manifiesto las características litológicas de los materiales involucrados. Éstos consisten en basalto fracturado, recubierto por un horizonte limo–arcilloso de espesor variable. Además, la precipitación, la evaporación, el contenido de humedad en el suelo y la variación del nivel freático se cuantificaron mediante un conjunto de dispositivos automáticos. Los valores obtenidos se emplearon para calibrar y establecer las condiciones iniciales y de frontera de un modelo numérico de flujo y transporte de contaminantes. En dicho modelo se conceptualiza la zona de estudio como una serie de capas equivalentes a medios porosos continuos. Los resultados obtenidos a través de las diferentes simulaciones evidencian que la naturaleza de la primera capa es determinante en la tasa de infiltración. En los últimos dos meses de observación hubo un sellado parcial de las fracturas debido a los sólidos en suspensión presentes en las aguas residuales. En principio esto contribuyó a retardar el proceso de infiltración. Las simulaciones para condiciones extremas muestran que un contaminante conservativo (sin reacciones) podría atravesar la zona vadosa en un rango de dos a cuatro años para el peor y mejor de los casos considerados, respectivamente.

Palabras clave: infiltración, modelización, aguas residuales, zona vadosa, fracturas, valle de Toluca, México.

ABSTRACT

In the Toluca valley, important surface displacements have occurred in response to the activation of recently generated fractures, which locally condition the infiltration of contaminated surface water. This is the case of a fractured zone in the community Santiago Tianguistenco (Mexico State), to which wastewater, mainly domestic, from about 14,000 habitants is discharged through a canal. The objective of this paper is to quantify the infiltration process in the unsaturated zone of the mentioned fractured zone by combining technological development and numerical modeling, both supported by intense field work. In the course of the study, vertical electrical soundings were carried out and a monitoring well was drilled, which provided information to define the lithological features of the materials involved. These consist of fractured basalt, covered by a silt–clay horizon of variable thickness. Also, precipitation, evaporation, soil water content and water table variations were quantified by automatic devices, and the resulting values were used to calibrate and set the initial and boundary conditions in a numerical model of flow and transport. In this model, the studied zone is conceptualized as consisting of several layers, similar to a continuum porous media. The results obtained in the different simulations show that the nature of the first layer is determinant for the infiltration rate. In the last two months of observation, the fractures were partially sealed as a result of settling of solids suspended in the wastewater, a process that contributed to retard the infiltration process. The simulations for extreme conditions show that a conservative contaminant (without reactions) would go through the vadose zone in a time range of two to four years for the worst and best case, respectively.

INTRODUCCIÓN

El agua tiene un papel determinante en el sustento y desarrollo de los seres vivos, por lo que desde sus inicios los asentamientos humanos han dependido en gran medida de la proximidad a fuentes de abastecimiento. Esta circunstancia ha provocado que la demanda de agua de buena calidad se haya incrementado con el paso del tiempo paralelamente al incremento demográfico. Por eso, sin el agua subterránea actualmente no sería posible garantizar el abastecimiento doméstico de agua potable, ni el destinado a la agricultura y la industria. De hecho, en muchas regiones la mayor parte del agua tiene origen subterráneo, hasta un 80% en Europa y Rusia e incluso más en el Medio Oriente y en el norte de África (IYPE, 2005).

En México existen graves problemas de abastecimiento de agua en cantidad y calidad en el sector central de la república, en parte debido a la alta densidad de población y al acelerado crecimiento (UNESCO, 2006). Lo anterior ha provocado una extracción intensiva de los acuíferos, lo cual a su vez ocasiona problemas colaterales importantes. Entre ellos, se puede mencionar la subsidencia asociada a la extracción del agua subterránea, de la que existen ejemplos bien documentados en algunas cuencas de México (Ortega–Guerrero et al., 1993; Ortega–Guerrero et al., 1999; Rudolph y Frind, 1991; Rivera et al., 1991). Sin embargo, hay muy pocos estudios referentes al valle de Toluca (Figura 1; Rudolph et al., 2006), en el que se ubica la zona de estudio del presente trabajo. Precisamente es allí donde se han observado en los últimos años importantes desplazamientos del terreno asociados a la aparición de fracturas con desigual comportamiento dinámico, constituyendo localmente vías preferenciales que favorecen la infiltración de aguas superficiales contaminadas (Figueroa–Vega, 2004).

En este contexto, el presente trabajo tiene como objetivo la estimación del tiempo que tarda un contaminante conservativo en atravesar una zona no saturada fracturada en la que se vierten aguas residuales, hasta incorporarse al acuífero subyacente de la región de Santiago Tianguistenco. Esta zona cuenta con la particularidad de estar situada en la vecindad de un pozo de abastecimiento de agua potable. Para alcanzar este objetivo, se consideró la estratigrafía de la zona no saturada y se cuantificaron los volúmenes de entrada que afectan al sistema acuífero subyacente para su incorporación en un modelo de flujo y transporte de contaminantes para medios porosos variablemente saturados desarrollado por el United States Geological Survey (USGS) de los EEUU denominado UNSAT (Fayer, 2000; Lappala et al., 1987).

DESCRIPCIÓN DE LA ZONA DE ESTUDIO

Desde 1942, la cuenca del valle de Toluca (Figura 1) ha provisto de agua potable a la ciudad de México (Sabalcagaray, 1981; UAEM, 1993). Inicialmente esto fue posible mediante la captación de sus principales manantiales. A partir de 1951, la extracción del agua subterránea se realizó a través del sistema Lerma. En la actualidad este sistema está constituido por un conjunto de 236 pozos de bombeo perforados paralelamente a la frontera oriental de la cuenca (UAEM, 1993).

La excesiva extracción del agua subterránea en la cuenca ha ocasionado, entre otros problemas, el descenso de los niveles piezométricos y el cambio en los patrones del flujo subterráneo, fenómeno sólo observable externamente por la lenta pero progresiva disminución de los cuerpos de agua lacustres existentes en la cuenca. Esto ha traído como consecuencia una pérdida del almacenamiento en el acuífero, lo que ha dado lugar a significativos descensos del terreno en los últimos 20 años, fracturas en el terreno, así como colapsos en la superficie (Figueroa–Vega, 2004).

Estas prácticas de manejo del agua residual representan un grave peligro para la calidad de los recursos hídricos subterráneos que constituyen la única fuente de abastecimiento para los núcleos poblacionales inmediatos. Así, a tan solo 900 metros de la mencionada zona fracturada opera un pozo de explotación de agua potable (Figura 2d). La aparición de estas fracturas en superficie es relativamente reciente, tal y como se desprende de un sencillo análisis histórico. En la fotografía aérea mostrada en la Figura 2a y correspondiente al año 2000, no se advierte fractura alguna afectando a la zona en cuestión; mientras que en la Figura 2e, tomada en el año 2005 sí se pueden observar.

MÉTODOS Y MATERIALES

La metodología empleada en este trabajo consistió en tres etapas. En la primera, se abordó la recopilación de información documental de la zona de estudio. Durante la segunda etapa, se emprendió una campaña de adquisición de parámetros hidrogeológicos. Incluyó la instalación de cinco dispositivos para el monitoreo de la precipitación, la evaporación, el caudal de descargas de aguas residuales, el contenido de humedad y el nivel freático. Por último, la tercera etapa consistió en el empleo de los datos obtenidos para la generación de un modelo hidrogeológico conceptual, cuya información se incorporó en un modelo de flujo y transporte de contaminantes a través de la zona no saturada.

Recopilación de la información histórica

Con el fin de ubicar el contexto geológico del fracturamiento del área de estudio y establecer la ubicación del pozo de monitoreo, se recopiló la información concerniente a la ubicación de pozos de abastecimiento, sondeos eléctricos verticales, mapas geológicos y fotografías aéreas disponibles. En lo que respecta a los datos hidrológicos, se analizó la existencia y ubicación de las estaciones meteorológicas y de aforo, a fin de identificar si contaban con datos de precipitación, evaporación y caudal en corrientes de aguas superficiales. Puesto que no hay estaciones de esta naturaleza dentro de la zona de estudio, los datos se determinaron en el campo especialmente para este proyecto y se cotejaron con datos de precipitación y evaporación de las estaciones metereológicas automatizadas del Servicio Meteorológico Nacional (SMN) más cercanas.

Caracterización hidrogeológica

Para determinar el tipo y el espesor de los diferentes niveles que conforman la zona no saturada se realizaron sondeos eléctricos verticales (SEV), cuya localización, orientación y aperturas interelectródicas A–B se muestran en la Figura 3.

Este método suministra un conocimiento de la variación vertical de la resistividad aparente de los materiales (ρ) mediante un dispositivo tetraelectródico en el que, simultáneamente a la introducción de una corriente eléctrica de intensidad conocida, se obtiene una correspondiente diferencia de voltaje. El ajuste posterior de los datos obtenidos con ayuda del programa Qwseln (Francia, Centre National de la Recherche Scientifique, Unité Mixte de Recherche 7619), permitió establecer la resistividad real de cada uno de los niveles existentes e indirectamente su naturaleza litológica, espesor y profundidad. Estas características se contrastaron y validaron por comparación con las secuencias estratigráficas obtenidas a partir del análisis petrológico de las muestras extraídas durante la perforación de un pozo de monitoreo. Éste se perforó en diciembre de 2005 especialmente para este trabajo. Su localización se muestra en la Figura 3.

Las variaciones del contenido de humedad en el suelo se establecieron mediante un dispositivo cuyas terminales de prueba fueron ubicadas a 0.5 m, 2 m y 3.8 m de profundidad respecto al nivel del suelo. Dicho instrumento almacena los valores de la conductividad eléctrica entre dos terminales contenidas en bloques de yeso. Para su calibración se efectuaron pruebas de laboratorio que incluyeron su caracterización en función de la humedad y la temperatura siguiendo métodos estándar de calibración para bloques de yeso (Dela, 2001; Salas–García et al., 2010). Se tomaron mediciones in situ de la resistencia eléctrica del área de instalación junto a muestras del mismo lugar para, posteriormente, aplicar el método gravimétrico expuesto en Sanders (1998). Éste consiste en correlacionar el contenido de humedad en el medio con la resistencia eléctrica entre las terminales de prueba.

Para medir el nivel freático se instaló un sensor de nivel en el interior del pozo de monitoreo. Dicho freatímetro mide el nivel del agua a través de un transductor que cuantifica la presión diferencial entre la columna de agua asociada y la presión atmosférica, según se describe en Salas–García et al. (2010). La precipitación se midió con un pluviómetro ubicado a 900 m de la zona de descargas de aguas residuales. La evaporación se calculó a partir de un tanque de evaporación tipo A (Allen y Pruitt, 1991). Para determinar las aportaciones que se descargan en la zona de fracturas procedentes del canal de aguas residuales se empleó un segundo medidor de nivel. Así, el caudal se calculó mediante una curva calibrada del nivel en el canal versus caudal. Los cinco dispositivos mencionados cuentan con la capacidad de efectuar determinaciones continuas a intervalos programables y almacenarlos para su descarga posterior en cualquier computadora personal. Los detalles relativos al diseño y funcionamiento de estos equipos se puede consultar en Salas–García et al. (2010). El conjunto de la información aportada por los dispositivos mencionados ha permitido definir un modelo conceptual de la dinámica hidrológica de la zona de estudio.

A fin de determinar la tasa máxima de infiltración en la superficie, se emprendió una campaña de siete mediciones realizadas con un permeámetro Guelph fabricado por la empresa Soil Moisture Corporation. Dicha tasa está determinada por la conductividad hidráulica en saturación (Sanders, 1998).

Modelización del flujo y transporte de contaminantes en la zona vadosa

En la última etapa de este trabajo se aplicó un modelo denominado UNSAT desarrollado por la USGS que, dadas las condiciones iniciales y de frontera, simula el flujo y el transporte de contaminantes en la zona vadosa (Lappala et al., 1987). Así mismo, de las distintas opciones actualmente existentes en la modelización que permiten describir el flujo y el transporte en la zona no saturada, en la presente investigación se optó por considerar a la zona de estudio como un conjunto de múltiples horizontes equivalentes a medios porosos continuos pero desigualmente permeables. A este respecto, tanto Larocque et al. (2000) como Ford y Williams (1989) indican que se puede aplicar este modelo a medios fracturados, incluso kársticos con baja disolución. Otros autores comparten este punto de vista (Berkowitz et al., 1988; Pulido–Bosch y Padilla–Benitez, 1998; Kresic et al., 1993; Estrela y Sahuquillo, 1997; Angelini y Dragoni, 1997; Finsterle, 2000; Royer et al. 2002; Xiang, 2005).

La forma convencional de la ley de Darcy no describe adecuadamente el flujo de agua en la zona no saturada, debido, por un lado, al rápido decrecimiento del valor de la conductividad hidráulica cuando el contenido de agua disminuye y, por otro, a la disminución del área transversal disponible para el flujo de agua en medios no saturados. Puesto que en la zona de estudio resulta patente la presencia de fracturas, es necesario recurrir a un planteamiento en el que se contemple el cálculo de la conductividad hidráulica Para lograr lo anterior, se aplicó la ecuación de Richards (1931), expresando la conductividad hidráulica, K, como una función del contenido de humedad, θ :

donde D(θ) representa la difusividad del agua (Guymon, 1994). Sin embargo, no hay relaciones universales para la conductividad hidráulica en un medio no saturado versus succión o contenido de humedad. A lo largo de la historia, se han propuesto varias relaciones empíricas (Baver et al., 1972; Childs y Collins–George, 1950a, 1950b; Brooks y Corey, 1964; Mualem, 1976; Van Genuchten, 1980). No obstante, para aplicar apropiadamente la Ecuación 1 a la zona de estudio, se requiere la estimación previa de θ, razón por la cual se instaló un medidor del contenido de humedad. Al introducir dicho parámetro en las relaciones mencionadas, se puede obtener K(θ). A su vez, para obtener la solución de la aproximación seleccionada, generalmente se recurre al método de los elementos finitos y sus variantes, como es el caso del modelo UNSAT (Lappala et al., 1987).

donde: θ es el contenido volumétrico de agua, adimensional; c es la concentración del constituyente químico, [ML^–3] (masa por unidad de volumen de agua); t el tiempo, [T]; = (∂/∂x)+(∂/∂y)+(∂/∂z) es el operador, [L^–1]; _h el tensor de dispersión hidrodinámica, [L²T^–1]; el vector de velocidad de flujo, [LT^–1] y SS es un término de fuente/sumidero, [ML^–3T^–1]. Bear (1972) describe el desarrollo de la ecuación anterior. Para el modelo empleado en este trabajo, c se consideró como la concentración normalizada de un contaminante conservativo (sin reacciones) para simular el agua superficial acumulada procedente del canal de aguas residuales. Por otra parte,

incluye el tensor de dispersión mecánica (_m ) y el de difusión molecular (). Además, si se considera que para este trabajo sólo se realiza la modelización en una dimensión (z), cualesquier otras componentes en x y en y se cancelan. De este modo, se tiene que el tensor de dispersión hidrodinámica en la dirección z (D_hz) puede escribirse como:

donde D_m = D_dτ, mientras que D_d es el coeficiente de difusión molecular del soluto en agua, [L²T^–1], y τ es la tortuosidad (adimensional) del medio poroso. El parámetro αL representa la dispersividad longitudinal del medio poroso, [L]; |v|, la magnitud del vector de velocidad, [LT^–1], y vz a la componente en dirección vertical del vector velocidad.

Para determinar la dispersión mecánica (D_m) se requiere establecer los valores de la tortuosidad (τ) y del coeficiente de difusión molecular (D_d). La tortuosidad (τ) es un factor adimensional, en un rango de 0.3 a 0.7 para la mayoría de los suelos (Van Genuchten y Wierenga, 1976). A partir de una revisión de los datos procedentes de medios no consolidados obtenidos por varios investigadores, Perkins y Johnson (1963) sugieren el valor τ de aproximadamente 0.7070, al igual que Gillham y Cherry (1982). Bear (1972), sugiere que el valor de 0.67 es el mejor para materiales granulares. Dada la naturaleza de los materiales de la zona de estudio, se empleó el último valor. Por otra parte, los valores del coeficiente de difusión molecular (D_d) y de dispersividad longitudinal (α_L) se determinaron a partir de un estudio de Jørgensen et al. (1998) en el que se emplearon trazadores para la determinación experimental de parámetros de transporte en medios fracturados similares a los de este trabajo.

El último término de la Ecuación 2 (SS), se puede dividir en dos categorías generales: (1) masa de soluto que es introducida por fuentes o removida por sumideros y (2) masa de soluto que se introduce o remueve a partir de reacciones químicas que ocurren en el agua dentro del sistema o que tienen lugar entre el agua y la fase sólida. Para la metodología de este trabajo sólo se consideró la primera categoría, ya que no se incluyó ningún tipo de reacción química en el agua del canal de aguas residuales o con el medio circundante.

Para una mejor descripción, los resultados obtenidos se han estructurado en dos categorías principales. En la primera se tratan las cuestiones concernientes a la caracterización hidrogeológica de la zona en estudio. En la segunda se describen las condiciones asumidas para el desarrollo del modelo y la integración de los datos previos en la modelización.

Caracterización hidrogeológica de la zona de estudio

Tanto la localización de los sondeos eléctricos verticales como los resultados de la campaña se muestran en las Figuras 3 y 4. En la primera se presenta la localización, orientación y amplitud interelectródica A–B máxima alcanzada en cada ensayo. En la segunda se grafican las curvas experimentales de campo con los correspondientes ajustes a curvas teóricas mediante el empleo del programa Qwseln. Además, se incluye el resultado de la interpretación que en ambos sondeos efectuados coincide con la presencia de cinco capas de desigual resistividad y espesor.

Los sondeos se realizaron en febrero de 2005, antes de la época de lluvias, para evitar posibles alteraciones debido al alto contenido de humedad en el terreno. En enero de 2006, se perforó un pozo de monitoreo en el marco de este proyecto hasta una profundidad de 80.3 m con el fin de hacer un seguimiento temporal del nivel freático local, y determinar la columna estratigráfica tipo, a partir de los testigos extraídos a cada metro de profundidad.

La interpretación obtenida con los sondeos eléctricos resultaron casi coincidente con la información proporcionada por el pozo de monitoreo, como se pone de manifiesto al correlacionar los datos de los sondeos con la estratigrafía del pozo, que en conjunto estaría representada por cinco niveles agrupados en tres tipos básicos de materiales, todos ellos integrantes de la zona no saturada.

El valor de la resistividad en la primera capa corresponde a un recubrimiento reducido (~1–2 m) integrado por un suelo parcialmente encostrado. Esto concuerda con los valores reportados por Samouëlian et al. (2005), que indican que los valores altos de resistividad corresponden a contenidos de humedad bajos. Le seguiría un segundo nivel de relleno de unos 10–12 metros de espesor compuesto por acumulación de limo y arena fina con arcilla procedente de la alteración de basalto y del transporte de sedimentos por efecto del escurrimiento superficial hasta la zona de estudio. A continuación se encuentran los materiales basálticos en desigual proceso de alteración con un espesor de unos 15–17 metros. A éstos le siguen materiales piroclásticos de una granulometría semejante a la arena.

Los materiales descritos coinciden espacialmente con los encontrados en los cortes litológicos de dos pozos de abastecimiento perforados en la localidad de Coatepec, a una distancia de 0.9 y 2 km del pozo de monitoreo de este proyecto (Figura 2d), cuyas columnas estratigráficas se muestran en la Figura 5.

En la Figura 6 se muestran los datos de precipitación medidos con el pluviómetro durante el periodo comprendido entre el 11 de junio y 9 de septiembre del 2006. Para validar dichos datos, se compararon con los procedentes de dos estaciones meteorológicas automatizadas del SMN ubicadas a 20 y 30 km, obteniéndose valores muy cercanos, lo que ratifica que los valores determinados son confiables. En lo que respecta a la evaporación, los valores oscilaron en un rango de 1 a 5.5 mm/día. En general, puede indicarse que los valores para ambos parámetros, obtenidos con los equipos instalados, estuvieron dentro de un margen menor al 6% de variación respecto a los obtenidos en las mencionadas estaciones. La evaporación potencial empleada en los cálculos posteriores es el resultado de la evaporación diaria medida en el tanque de evaporación multiplicada por un factor típico de 0.7 empleado para este tipo de instrumento (Allen y Pruitt, 1991). La razón por la que se consideró la evaporación potencial en la modelización en lugar de la evapotranspiración real o potencial obedece a que en la zona de descarga de aguas residuales no hay vegetación y su configuración es semejante a un cuerpo de agua superficial.

En la misma Figura 6 se incluye, para efectos de comparación, la variación del nivel freático. Es de interés notar que ambos parámetros tienen un comportamiento ascendente antes del 4 de agosto, y que tras esa fecha el nivel freático comienza a descender, aun cuando la precipitación continúa. Además, se puede apreciar un desfasamiento entre los picos de la precipitación y los del nivel freático, producto del tiempo de tránsito a través de la zona vadosa. Puesto que la zona de estudio es una zona de recarga muy cercana al parteaguas del acuífero subyacente, la componente atribuida a la recarga lateral puede considerarse despreciable.

Los datos presentados en la Figura 6 permiten establecer que el flujo ocurre principalmente por la matriz porosa. En varios estudios se ha empleado el tiempo de respuesta del nivel freático a eventos de precipitación intensas como criterio para definir si el flujo preferencial ocurre principalmente a través de la matriz porosa o por las fracturas (Mathias et al., 2005). Una respuesta lenta se asocia a un flujo a través de la matriz porosa; recíprocamente, una respuesta rápida se atribuye a un flujo preferencial que ocurre principalmente por las fracturas del medio poroso o bien a un mecanismo de "desplazamiento de pistón" (Mathias et al., 2005; Headworth, 1972). En este último, el agua que ingresa al acuífero ha sido desplazada por la presión procedente de la parte superior de la zona vadosa (Price et al., 1993).

El caudal del canal de aguas residuales que descarga en la zona de fracturas muestra una respuesta inmediata a la precipitación como se puede observar en la Figura 7. El análisis de dicha gráfica indica que, después de una lluvia intensa, el canal drena el agua resultante en ese día y parte del siguiente para regresar posteriormente a su caudal base de alrededor de 3.5 m³/min.

Modelo conceptual

Con los datos mencionados se generó un modelo conceptual unidimensional, cuyo esquema representativo se muestra en la Figura 8. De acuerdo con lo expuesto, la estratigrafía del terreno se dividió en cinco capas: las dos primeras de porosidad intergranular restringida y las otras tres de porosidad doble debido al medio granular y el medio fracturado, en este último caso se debe a la moderada alteración del basalto. En el período de observación de dos años hubo cambios en las dimensiones de las fracturas en superficie. Dada la geología de la zona, esto se debe principalmente a la parcial desecación de los niveles superiores. Para reflejar los extremos de dichos cambios, el modelo conceptual considera las primeras dos capas como medios porosos equivalentes. Cada uno de ellos se asocia con dos límites de conductividad hidráulica que simulan los comportamientos observados.

En la Tabla 1 se han incluido las designaciones de los materiales de cada capa considerando a su vez diferentes escenarios en función del desigual grado de fracturación que se puede presentar. En la Tabla 2 se presentan los valores de los parámetros asociados a cada uno de los materiales de la Tabla 1. En la segunda columna de la Tabla 2 se presenta, para efectos de comparación, la correspondencia entre los valores empleados en las simulaciones de este proyecto y aquellos reportados por Freeze y Cherry (1979) así como en el programa informático Rosseta (Schaap, et al., 2001). En la última columna se consignan los valores de porosidad eficaz empleados en este trabajo, mismos que son similares a los reportados por Sanders (1998).

Calibración

Todas las simulaciones se calibraron con los datos obtenidos durante el periodo de observación de marzo a septiembre de 2006. En la Figura 9 se comparan los resultados del modelo de Van Genuchten (1980) respecto a los datos de saturación obtenidos en el campo a 3.8 m de profundidad, constatándose que los valores de los parámetros de frontera de segundo orden definidos para la parte superior del modelo fueron elegidos de forma apropiada (Figura 8).

Una vez calibrado el modelo con la información hidrometeorológica de campo generada durante el proyecto, ésta se sustituyó por los datos de la precipitación y la evaporación diaria promedio correspondientes a un período de 15 años de la estación climatológica de Santiago Tianguistenco, ubicada a 3.5 km de la zona estudiada. Es necesario aclarar que estos datos se emplearon para la validación del modelo y no para calibrarlo. Esto se debe a que la calibración requirió los datos obtenidos durante el año 2006 y, de hecho, la estación dejó de operar en 1990 (SMN, 2000). El resultado de este último análisis puso de manifiesto que la precipitación no es un parámetro tan sensible como lo puede ser la conductividad hidráulica del medio.

El promedio de la precipitación acumulada de la estación climatológica de Santiago Tianguistenco, para el período comprendido entre el 11 de junio y el 10 de septiembre de los años 1980 a 1990, tiene una variación del 5.7% respecto a la precipitación acumulada medida en la zona para el mismo período durante el año 2006. Por otra parte, el valor de la conductividad hidráulica para cada uno de los escenarios simulados, puede variar hasta en un 35%, dependiendo de la sección estratigráfica en cuestión.

Escenarios simulados

De los tres escenarios considerados (Tabla 1), el intermedio reúne las características necesarias para reproducir las circunstancias y condiciones durante el año 2006. Se estableció la frontera superior de tipo 2 (flujo) como lo indica la Figura 8. La lámina de infiltración en dicha frontera es el resultado de la precipitación sumada al aporte del canal de aguas residuales menos la evaporación potencial. En cuanto a la frontera inferior, de tipo 1, estuvo determinada por el nivel freático. Con las fronteras definidas del sistema fue posible generar los resultados de la modelización de flujo y transporte presentados en las Figuras 10 a 12 [11]. Aunque el paso de tiempo empleado para las simulaciones es de una resolución diaria, para efectos de comparación, en dichas figuras se presentan los resultados correspondientes a los meses de marzo, mayo, julio y septiembre. En la Figura 10 se muestra la variación del grado de saturación de los niveles involucrados en función de la profundidad para los meses mencionados (a), junto con la simulación de transporte (b) y su perfil estratigráfico asociado (c). Respecto al grado de saturación, se confirma que los valores más altos se alcanzan en los dos primeros niveles correspondientes a materiales escasamente consolidados con una conductividad hidráulica menor respecto a los inferiores. Además, el considerar sólo la componente advectiva y difusiva de un potencial frente de contaminación constituye el caso más crítico, cuando el medio ya no tiene capacidad de retrasar o atenuar los contaminantes.

En las curvas de saturación de la Figura 10 correspondientes a los meses de enero y marzo se puede apreciar la variación de este parámetro en los diversos estratos en función de la profundidad. Lo anterior confirma que la variación en la saturación depende, principalmente, del agua retenida en el medio que no drena por gravedad. De este modo, los primeros estratos presentan un alto grado de saturación respecto de los inferiores. La baja saturación de estos últimos para los meses mencionados indica que son capaces de drenar con suficiente rapidez el agua procedente de las capas superiores, de tal modo que su saturación no se incrementa sustancialmente.

La última etapa en la modelización comprende la simulación de los otros dos escenarios que representan las circunstancias extremas del sistema. La diferencia entre los escenarios designados como severo y conservador, radica en que para el primero se considera el conjunto de materiales como un medio equivalente que permite una conductividad hidráulica mayor. De este modo, el denominado intermedio, ya tratado, correspondería a las circunstancias hidrogeológicas observadas en el año 2006, mientras que el denominado severo, o de alto riesgo, coincidiría a grosso modo con la situación hidráulica prevaleciente en los dos años anteriores, 2004 y 2005, durante los cuales no se produjo inundación alguna en el área de fracturas estudiada. Por su parte, el escenario conservador, o de bajo riesgo, no corresponde a ningún año del que se tengan datos, por lo que sus características se dedujeron a partir de las diferencias existentes entre el año 2005 y 2006. Así, por ejemplo, la diferencia en el valor de los parámetros de la conductividad hidráulica entre el año 2005 y 2006 es aproximadamente la misma en magnitud que entre el año 2006 y el escenario conservador, pero de signo opuesto. Los valores empleados para los tres escenarios mencionados están dentro del rango empleado en otros trabajos de investigación, tales como Jørgensen et al. (1998), Rosenbom et al. (2009) y Schaap et al. (2001).

Para evitar ambigüedades, es necesario señalar la existencia de diferentes tiempos de simulación para cada escenario, a pesar de que en todos ellos se modela el comportamiento de un contaminante que comenzaría a infiltrarse en un instante común (enero de 2004). Si durante todo el tiempo que dura la simulación las entradas al sistema y el grado de conexión de las fracturas correspondieran con las del año 2005, el resultado sería el escenario severo, y si el sistema se comportara como en el 2006, el escenario sería el intermedio. Esta diferencia se debe al desigual grado de actividad y conexión de las fracturas, a la escasez de datos de conductividad hidráulica y a que el programa adaptado para la modelización UNSAT (Fayer, 2000) no permite modificar los parámetros del medio poroso en función del tiempo. Por consiguiente, dado que las dos primeras capas determinan la cantidad de agua que fluye hacia los estratos inferiores, los diferentes escenarios presentados simulan una cantidad variable de fracturas mediante distintos valores de los parámetros del medio en los primeros dos estratos, según se muestra en las Tablas 1 y 2.

Las simulaciones de transporte se recrearon a partir de enero de 2004, momento en que se desviaron las aguas del canal a la zona de fracturas. En la Figura 11 se muestran los resultados de la modelización para el escenario severo durante el cual la zona no se inundó, lo que se consigue simulando un mayor número de fracturas activas, es decir, aumentando la conductividad hidráulica y la porosidad en los dos primeros estratos. Desde un punto de vista cualitativo puede apreciarse en la Figura 11 un comportamiento comparable con el escenario intermedio (Figura 10), en especial, en lo que se refiere al papel de barrera que, incluso en estas adversas condiciones, siguen presentando los dos primeros estratos, y a pesar de que cuantitativamente muestran una considerable disminución en la saturación.

Este último aspecto no debe interpretarse en el sentido de que la cantidad de agua que ha ingresado en el medio es menor. Más bien, el agua es capaz de pasar más rápidamente por estos niveles, impidiendo que la matriz porosa alcance un mayor grado de saturación. Este incremento en la velocidad tiene implicaciones en lo que respecta al transporte de contaminantes, como queda reflejado en las Figuras 10 y 11. En efecto, el análisis comparativo de ambas permite visualizar adecuadamente el diferente comportamiento dinámico del transporte, corroborándose la existencia de una diferencia temporal en el momento de llegada del mismo frente potencialmente contaminante al acuífero subyacente de tres y dos años para los escenarios intermedio y severo respectivamente. A este respecto hay que indicar que las condiciones de simulación para ambos escenarios fueron las mismas, incluyendo los datos de precipitación, evaporación potencial y caudales de descarga de aguas residuales, y solo difieren en los valores de la conductividad hidráulica asignada a los dos primeros estratos.

En el otro extremo se tiene el escenario conservador que, en el mejor de los casos, se alcanzaría si tuviera lugar el relleno y colmatación de las fracturas por efecto de los arrastres de sólidos en suspensión procedentes del canal de aguas residuales y de las zonas que circundan la zona de fracturas. Teniendo en cuenta estos nuevos condicionantes, en la Figura 12 se muestran los resultados de la correspondiente simulación efectuada, asimilando la matriz de los dos primeros niveles a un conjunto de naturaleza limo–arcillosa y baja conductividad hidráulica, lo que da lugar a un aumento de la saturación del medio y, en consecuencia, a un retardo mayor del transporte. De este modo, un contaminante no reactivo alcanzaría el nivel freático en unos cuatro años. Por tanto, entre los tres diferentes escenarios planteados existe una diferencia temporal de dos a cuatro años para que un contaminante, sin reacciones, se incorpore al flujo subterráneo del acuífero subyacente.

CONCLUSIONES

La contaminación del acuífero subyacente a un área endorreica fracturada integrante del municipio de Santiago Tianguistenco, México, parece estar asociada a las descargas de aguas residuales que se vierten en dicha zona. Entre los principales factores condicionantes están la época del año, la apertura y conectividad de las fracturas, las características geológicas del medio fracturado, y el posible grado de evolución en el tiempo de dichas características. En el presente estudio se centró la atención en la caracterización de la zona no saturada al objeto de determinar bajo qué condiciones el acuífero subyacente pudiera mostrar un mayor riesgo de contaminación. En virtud de lo anterior, se obtuvieron las siguientes conclusiones.

El análisis histórico de fotografías aéreas de la zona de estudio reveló que las fracturas en superficie se generaron alrededor del año 2001. Sin embargo, la cantidad y apertura de éstas disminuyó en el año 2006, lo que apunta a un comportamiento dinámico y posiblemente cíclico en el futuro. Las observaciones de campo indican una apertura mayor en estiaje, mientras que la combinación de niveles de alta saturación así como la presencia de sólidos en suspensión disminuye la cantidad y conectividad de las fracturas.

En lo que concierne a la dinámica de la infiltración, las simulaciones de flujo indican que los niveles integrantes de los primeros doce metros a partir de la superficie del terreno son determinantes para la tasa de infiltración. En las simulaciones de transporte, considerando una fuente de contaminación puntual constante situada en la superficie, se ha podido establecer que el tiempo que tarda un contaminante conservativo en alcanzar el nivel freático del acuífero sin procesos de absorción o retardo oscila entre dos y cuatro años. Teniendo en cuenta que las fracturas se hicieron evidentes en superficie en el año 2001, de acuerdo con las simulaciones, actualmente las potenciales sustancias contaminantes ya se han incorporado al flujo subterráneo del acuífero subyacente.

A partir de los resultados expuestos, se recomienda finalizar las obras de prolongación del canal de aguas residuales, suspendidas alrededor del 2003, al objeto de evitar las descargas en la zona de estudio. Mientras tanto, una medida temporal para disminuir la tasa de infiltración consistiría en rellenar las fracturas con el mismo material arcilloso de la zona de estudio. Por otra parte, y puesto que la topografía local favorece el drenaje hacia la zona de fracturas, ubicada en el punto más bajo de una pequeña cuenca endorreica, un proyecto alternativo de interés sería su aprovechamiento para la instalación de un sistema de recarga artificial. La principal ventaja, al margen de los lógicos gastos derivados del tratamiento del agua para reducir la concentración de contaminantes y el contenido de sólidos en suspensión, de modo que su calidad se encuentre dentro de la norma, sería el bajo costo de la infraestructura necesaria para la captación y conducción del agua.

AGRADECIMIENTOS

Los autores agradecen por el financiamiento de este proyecto al Consejo Nacional de Ciencia y Tecnología (CONACyT), a la Comisión Geológica del Canadá, Recursos Naturales Canadá, División Québec (CGC–Québec/INRS), a la Universidad Autónoma del Estado de México (proyecto UAEM–2229) y al Departamento de Geodinámica del País Vasco. También por las facilidades otorgadas por el Ministère des Relations Internationales de Québec, Canadá y a la Dirección de Política y Cooperación Internacional en Ciencia y Tecnología del CONACyT, así como a la Comisión Nacional del Agua (CONAGUA) por hacer disponible la información que tiene a su cargo. Además reconocemos los valiosos comentarios y sugerencias del Dr. Peter Birkle (Editor Científico) así como de los revisores que mejoraron el contenido del artículo.

REFERENCIAS

Allen, R.G., Pruitt, W.O., 1991, FAO–24 Reference Evapotranspiration Factors: Journal of Irrigation and Drainage Engineering, 117(5), 758–773.         [ Links ]

Angelini, P., Dragoni, W., 1997, The problem of modeling limestone springs: The case of Bagnara (North Apennines, Italy): Ground Water, 35(4), 612–618.         [ Links ]

Bear, J., 1972, Dinamics of Fluids in Porous Media. Dover Publications: New York., American Elsevier Publishing Co., 764 pp.         [ Links ]

Berkowitz, B., Bear, J., Braester, C., 1988, Continuum models for contaminant transport in fractured formations: Water Resources Research, 24(8), 1225–1236.         [ Links ]

Baver, L.D., Gardner, W.H., Gardner, W.R., 1972, Soil physics: New York. John Wiley and Sons, 498 pp.         [ Links ]

Brooks, R.H., Corey, A.T., 1964, Hydraulic properties of porous media: Fort Collins, CO., Colorado State University, Hydrology Papers, 3, 27 pp.         [ Links ]

Childs, E.C., Collis–George, N., 1950a, The permeability of porous materials: Proccedings of the Royal Society of London, A201, 392–405.         [ Links ]

Childs, E.C., Collis–George, N., 1950b, Movement of moisture in unsaturated soils: Transactions of the 4th International Congress of Soil Sciences, Ámsterdam, v.1, 1–4.         [ Links ]

Dela, B.F., 2001, Measurement of soil moisture using gypsum blocks: Statens Byggeforskiningsinstitut. Danish Building and Urban Research, Documentation 004, 30 pp.         [ Links ]

International Year of Planet Earth (IYPE), 2005, Groundwater – Reservoir for a thirsty planet?: UNESCO, International Year of Planet Earth, 16 pp.         [ Links ]

Estrela, T., Sahuquillo, A., 1997, Modeling the response of a karstic spring at Arteta Aquifer in Spain: Ground Water, 35(1), 18–24.         [ Links ]

Fayer, M.J., 2000, UNSAT–H Version 3.0: Unsaturated soil water and heat flow model. PNNL–13249, Pacific Northwest Laboratory, Richland, Washington, 184 pp.         [ Links ]

Figueroa–Vega, G.E., 2004, El agrietamiento en la ciudad de Toluca: Consejo Consultivo del Agua, 45 pp.         [ Links ]

Finsterle, S. 2000, Using the continuum aproach to model unsaturated flow in fractured rock: Water Resources Research 36(8), 2055–2066.         [ Links ]

Ford, D., Williams, P., 1989, Karst Geomorphology and Hydrology: London, Academic Division of Unwin Hyman, 601 pp.         [ Links ]

Freeze, R.A., Cherry, J.A, 1979, Groundwater: Englewood Cliffs, N.J., Prentice–Hall, 604 pp.         [ Links ]

Guymon, G.L., 1994, Unsaturated Zone Hydrology: Englewood Cliffs, N.J., Prentice Hall, 210 pp.         [ Links ]

Gillham, R.W., Cherry, J.A., 1982, Contaminant migration in saturated unconsolidated geologic deposits, en Narasimhan, N.T. (ed.), Recent Trends in Hidrogeology: Boulder, CO, The Geological Society of America. Special Paper 189, 31–62.         [ Links ]

Headworth, H.G., 1972, The analysis of natural groundwater level fluctuations in the chalk of Hampshire: Journal of the Institution of Water Engineers and Scientists, 26, 107–124.         [ Links ]

Jørgensen, P.R., Mc Kay, L.D., Spliid, N.H., 1998, Evaluation of chloride and pesticide transport in a fractured clayey till using large undisturben columns and numerical modeling: Water Resources Research, 34(4), 539–553.         [ Links ]

Kresic, N., Kukuric, N., Zlokilica, M., 1993, Numerical versus stochastic modelling of water balance and minimum discharge of a karst hydrogeological system, en Günay, G., Johnson, A.I., Back, W. (eds.), Hydrogeological Processes in Karst Terranes: International Association of Hydrological Sciences, 207, 253–259.         [ Links ]

Lappala, E.G., Healy, R.W., Weeks, E.P., 1987, Documentation of computer program VS2D to solve the equations of fluid flow in variably saturated porous media: United States Geological Survey, Water–Resources Investigations Report 83–4099, 184 pp.         [ Links ]

Larocque, M., Banton, O., Razack, M., 2000, Transient–state history matching of a karst aquifer ground water flow model: Ground Water, 38(6), 939–946.         [ Links ]

Mathias, S.A., Butler, A.P., McIntyre, N., Wheater, H.S., 2005, The significance of flow in the matrix of the Chalk unsaturated zone: Journal of Hydrology, 310(1–4), 62–77.         [ Links ]

Mualem, Y., 1976, Hysteretical models for prediction of the hydraulic conductivity of unsaturated porous media: Water Resources Research 12(6), 1248–1254.         [ Links ]

Ortega–Guerrero, A., Cherry, J. A., Rudolph, D.L., 1993, Large–scale aquitard consolidation near Mexico City: Groundwater, 31(5), 708–178.         [ Links ]

Ortega–Guerrero, A., Rudolph, D.L., Cherrry, J.A., 1999, Analysis of long–term land subsidence near Mexico City. Field investigations and predictive modeling: Water Resources Research, 35(11), 3327–3341.         [ Links ]

Perkins, T.K., Johnston, O.C., 1963, A review of diffusion and dispersion in porous media: Society of Petroleum Engineers Journal, 3, 70–84.         [ Links ]

Pulido–Bosch, A., Padilla–Benitez, A., 1988, Deux exemples de mod–élisation d'aquiferes karstiques espagnols: Hydrogéologie, 1988(4), 281–290.         [ Links ]

Price, M., Downing, R.A., Edmunds, W.M., 1993. The Chalkasan aquifer, en Downing, R.A., Price, M., Jones, G.P. (eds.), The Hydrogeology of the Chalk of North–West Europe: Oxford, Clarendon Press, 35–58.         [ Links ]

Richards, L.A., 1931, Capillary conduction of liquids through porous medium: Physics, 1(5), 318–333. doi:10.1063/1.1745010.         [ Links ]

Rivera, A., Ledoux, E., Marsily, G., 1991, Non–linear modelling of groundwater flow and total subsidence in the Mexico City aquifer aquitard system, en Johnson, A.I. (ed.), Land Subsidence, Proceedings of the Fourth International Symposium of Land Subsidence: International Association of Hydrological Sciences (IAHS) Publication 200, 45–48.         [ Links ]

Rosenbom, A.E., Terrien R., Refsgaard, J.C., Jensen K.H., Ernsten, V., Klint, K.E.S., 2009, Numerical analysis of water and solute transport in variably–saturated fractured clayey till: Journal of Contaminant Hydrology, 104 (1–4), 137–152.         [ Links ]

Royer P., Auriault, J.L., Lewandowska, J., Serres, C., 2002, Continuum modelling of contaminant transport in fractured porous media: Transport in Porous Media, 49(3), 333–359.         [ Links ]

Rudolph, R.S., Garfias, J., McLaren, R.G., 2006, Significance of enhanced infiltration due to groundwater extraction on the disappearance of a headwater lagoon system: Toluca Basin, Mexico: Hydrogeology Journal, 14(12), 115–130.         [ Links ]

Rudolph, D.L., Frind E.O., 1991, Hydraulic response of highly compressible aquitards during consolidation: Water Resources Research, 27(1), 17–30.         [ Links ]

Sabalcagaray, M.D., 1981, Érase una vez Chignahuapan: la primera de las tres lagunas de Lerma: Toluca, México, Boletín del Archivo General del Estado de México, No. 9.         [ Links ]

Salas–García, J., 2007, Caracterización y modelación numérica unidimensional de la infiltración en la zona vadosa fracturada en Santiago Tianguistenco, México: Universidad Autónoma del Estado de México, Centro Interamericano de Recursos del Agua. (CIRA), Tesis de maestría, 127 pp.         [ Links ]

Salas–García, J., Gárfias, J., Llanos, H., Martel, R., 2010, Diseño y aplicación de instrumentación para la caracterización hidrometeorológica e hidrogeológica: Boletín de la Sociedad Geológica Mexicana, 62(2), 233–247.         [ Links ]

Samouëlian, A., Cousin, I., Tabbagh, A., Bruand, A., Richard, G., 2005, Electrical resistivity survey in soil science: a review: Soil & Tillage Research 83, 173–193.         [ Links ]

Sanders, L., 1998, A Manual of Field Hydrogeology: New Jersey, U.S.A., Prentice Hall, 381 pp.         [ Links ]

Schaap, M.G., Leij, F.J., Van Genuchten, M.Th., 2001, Rosetta: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions: Journal of Hydrology, 251,163–176.         [ Links ]

Servicio Metereológico Nacional (SMN), 2000, Base de Datos/Programa ERIC II: Instituto Mexicano de Tecnología del Agua (IMTA).         [ Links ]

United Nations Educational, Scientific and Cultural Organization (UNESCO), 2006, Water: a shared responsibility: The United Nations World Water Development Report 2, 584 pp.         [ Links ]

Universidad Autónoma del Estado de México (UAEM), 1993, Problemática ambiental de los recursos hídricos en la cuenca alta del río Lerma: Toluca, Estado de Mexico, UAEM, Coordinación General de Investigación y Posgrado, Seminario Internacional Sobre el Ambiente, v.1, 170–181.         [ Links ]

Van Genuchten, M., 1980, A closed–form equation for predicting the hydraulic conductivity of unsaturated soils: Soil Science Society of America Journal, 44, 892–898.         [ Links ]

Van Genuchten, M., Wierenga, P.J., 1976, Mass transfer studies in sorbing porous media. I. Analytical Solutions: Soil Science Society of America Journal, 40(4) 473–480.         [ Links ]

Xiang, Y.Y., 2005, Generalized equivalent–continuum method for modeling variably saturated seepage flow in fractured porous media: Rock and Soil Mechanics 26(5), 750–754.         [ Links ]