GERARDO ADRIAN RICCARDI ERIK DANIEL ZIMMERMANN Investigador Consejo Investigaciones UNR Investigador Asistente CONICET e_mail: riccardi@fceia.unr.edu.ar e_mail: erikz@fceia.unr.edu.ar Centro Universitario Rosario de Investigaciones Hidroambientales Facultad de Cs. Exactas, Ingeniería y Agrimensura. Universidad Nacional de Rosario Riobamba 245 bis. 2000 Rosario. ARGENTINA telefax 54 (0)341 480 8541 e_mail: curiham@fceia.unr.edu.ar

The development and application of two-dimensional mathematical model for groundwater flow in saturated porous media is described. Model structure is based on cell schemes, frequently used in free surface flow simulations. This work is part of a project to develop an integral system of soil-vegetation- atmosphere modelling, which will be used for dynamic simulations of the dominant hydrological processes in the region. Conceptual scheme, mathematical formulation, numerical solution and applications of the model are described. The results were compared with another, which were computed by means of methodologies and models of grateful quality. The problems analysed were: an aquifer between rivers, a two-dimensional flow example and an aquifer with extraction of water by means of pumping in deep wells. A study of quantitative affectation of aquifer in a rural area of the region subject to an irrigation project is also described. The achieved results indicate that the approach by means of cell schemes is satisfactory for the simulation of the movement of groundwater flow in porous media.

Se describe el desarrollo y aplicación de un modelo matemático bidimensional para escurrimiento subterráneo en medio poroso saturado. La estructura del modelo se basa en esquemas de celdas también utilizados en la simulación de escurrimiento a superficie libre. Este trabajo forma parte del proyecto de desarrollo de un sistema de modelación integral dentro del complejo suelo vegetación-atmósfera, apto para la simulación temporal, continua de los procesos hidrológicos dominantes en Ia región. Se plantea el modelo conceptual, la formulación y solución numérica, y una síntesis de testeos donde se compararon sus resultados con los que se computaran mediante metodologías y modelos de calidad reconocida. Los testeos comparativos aquí documentados comprenden: un macizo interfluvio, una gran área de flujo bidimensional y un acuífero con extracción de agua mediante bombeo en pozo profundo. Además se describe un estudio deafectación cuantitativa de acuífero en una zona rural de la región donde se proyecta establecer un dispositivo de riego intensivo. Los resultados logrados indican que la aproximación mediante esquemas de celdas es satisfactoria para la simulación del movimiento de flujo subterráneo en medio poroso saturado.

El movimiento de flujo subterráneo tiene características claramente bidimensionales y registra variaciones pequeñas y lentas de alturas a través del tiempo en tanto que los mecanismos friccionales y gravitacionales son preponderantes en el movimiento. En función de las fuerzas dominantes en el movimiento de flujo, puede afirmarse que resulta adecuado para su simulación el planteo de las hipótesis y ecuaciones gobernantes en que se basan los esquemas de celdas. Estos modelos simulan el movimiento bidimensional mediante el intercambio de flujo entre celdas con cualquier dirección contenida en el plano, pero con leyes de intercambio de tipo unidimensional.

Este trabajo específico se enmarca en un proyecto general de investigación que tiene como objetivo el desarrollo de un sistema de modelación dentro del complejo suelo-vegetación-atmósfera, apto para la simulación continua en el tiempo (Zimmermann y Riccardi, 1995), comprendiendo los procesos de: precipitación, evapotranspiración, escurrimiento superficial, infiltración, percolación y escurrimiento subterráneo. A la fecha, los autores ya han desarrollado y corroborado modelos que permiten la simulación de escurrimiento superficial (Riccardi, 1994/1997 y 1997), infiltración y evapotranspiración (Zimmermann y Riccardi, 1995). Estos modelos fueron estructurados de acuerdo a esquemas de unidades espaciales (celdas) interconectadas, por lo que se buscó generar para el caso de escurrimiento subterráneo un algoritmo de simulación de igual estructura que los correspondientes a los otros procesos a los efectos de permitir el ensamble de modelos dentro de un marco de total compatibilidad y escala de simulación.

Las ecuaciones que gobiernan el modelo son la de continuidad y la de descarga entre celdas vinculadas como ecuación de momento. La ecuación de continuidad se plantea para cada celda (Fig. 1) y se deriva a partir de la definición del incremento del volumen de agua almacenada desde consideraciones geométricas y condiciones de descarga (Cunge, 1975):

ASei es el área superficial efectiva de la celda; zi es la cota de agua en la celda i respecto a un plano de referencia; Pi(t) es el intercambio externo de caudal en la celda i; Qi,k es el caudal entre celdas i y k.

El modelo evalúa el intercambio de caudal entre celdas (Fig. 2) de acuerdo con la formulación de Darcy para flujo uniforme en medio poroso saturado:

hi,k = zm i,k - zmf i,k: es el nivel medio de agua entre celdas vinculadas, medida desde el fondo del acuífero; zf es la cota de fondo de acuífero; Ki,k es la permeabilidad media entre celdas del estrato; bi,k es el ancho total de la unión entre celdas; .zi,k = zi - zk es la diferencia de cotas de nivel de agua de los centros de celda; .xi,k es la distancia entre centros de celdas.

La función F i,k se determina en una instancia previa a la modelación en función de la geometría y las características de conductividad hidráulica del estrato y posteriormente es utilizada como parámetro de ajuste. En el caso de secciones compuestas con distintos coeficientes de permeabilidad a lo largo de su perfil vertical, el parámetro F i,k se determina para cada zona componente de la sección de similar característica sumándose luego para obtener un F i,k para la sección completa.

Formulación Numérica: Primeramente se explicita la función de descarga entre celdas para introducirla en la ec. de continuidad (1). Se asume que la descarga Qi,k(zi( t ), zk( t )) es una intermedia entre Qi,k ⁽ⁿ⁾ y Qi,k ⁽ⁿ⁺¹⁾ con t variando entre n .t <= t < =(n+1) .t. Es adecuado para la resolución el uso de un esquema implícito (Cunge, 1975):

.zi es la variación de cota de agua en celda i en un intervalo de tiempo dado; ..t es el intervalo de tiempo. El valor j es la cantidad total de celdas vinculadas a la celda i. En cada paso de tiempo se determinan los .z con los cuales se calculan los niveles de agua en cada celda:

z ⁿ⁺¹i = z ⁿi + .zi

Posteriormente se calcula en forma explícita de las descargas entre las celdas: Q ⁿ⁺¹i,k . La resolución numérica se realiza mediante un algoritmo basado en la resolución matricial por el método de Gauss-Seidel. El código computacional resultante fue denominado C1SUBT, ejecutado en lenguaje Fortran 77, versión Lahey F77L3 V4.0. No se realizaron versiones comerciales del modelo.

Condiciones de Borde e Iniciales: las condiciones de borde admitidas por el modelo son: (a) Nivel en función del tiempo: z(t); ( b). Caudal en función del tiempo: Q(t) y (c). Relación cota-caudal Q = f(z)

El modelo requiere la especificación de los niveles iniciales de agua en todas las celdas.

El problema de un interfluvio sujeto a una recarga en régimen permanente (Fig. 3) tiene solución analítica planteando la ecuación de Dupuit (Mijailov, 1985). La ley de descarga en un punto genérico x tiene la forma:

W es la intensidad de infiltración por unidad de tiempo y área; h1 y h2 son los espesores de la capa acuífera en las líneas de agua de los ríos; L1-2 es el ancho del interfluvio; k es la conductividad hidráulica saturada del macizo.

La curva de depresión puede determinarse analíticamente, sobre la base de la ecuación de continuidad para la abscisa x y la ecuación de momento deducida por Dupuit, arribándose a:

El planteo teórico previo se aplicó al perfil transversal del acuífero freático que subyace en la cuenca del arroyo Ludueña, respondiendo a un modelo conceptual del acuífero planteado previamente y validado mediante modelación matemática (Zimmermann, 1994a).

El macizo de 33 km de extensión está limitado por el arroyo Saladillo al sur y el arroyo San Lorenzo al norte, con cotas de agua estacionarias en los valores de 35 m y 55 m respectivamente. El ancho considerado fue 10 km. La permeabilidad media aproximada fue adoptada en 5 m/día. La recarga por infiltración fue establecido en 20 mm/año. Ambos valores presentan representatividad regional y fueron calibrados mediante la aplicación del modelo GW8 (N.U., 1989) (Zimmermann, 1994b). La zona en estudio fue discretizada en 23 celdas (Fig. 3). Los resultados computados por el C1SUBT y las ecuaciones teóricas de Dupuit se sintetizan en la siguiente tabla:

Tabla 1. Comparación de resultados interfluvio Ao. Saladillo - Ao. San Lorenzo

En la Fig. 4 se presenta la comparación de los perfiles hidráulicos de equilibrio alcanzados. La correlación entre los valores computados y los que arroja la metodología de Dupuit indica un coeficiente de correlación de 1.00 y un error estándar de correlación de 2.0e-12.

Se ensayó la simulación del flujo bidimensional en planta contrastando su respuesta con la de modelos matemáticos cuya aptitud haya sido verificada. Se seleccionó como herramienta de contraste el Groundwater Software (GW8)(N.U, 1989). El planteo matemático de resolución utilizadopor el GW parte de la ecuación diferencial general que describe la filtración inestable de una corriente plana sistemática bidimensional de aguas subterráneas en un medio poroso heterogéneo, conocida como ecuación de Boussinesq, que tiene la forma siguiente:

T es la transmisibilidad de la capa; H es el nivel de agua; t es la coordenada temporal; Q es la diferencia entre los caudales extraídos y recargados por unidad de área y S es el coeficiente de almacenamiento del medio poroso. La ecuación no tiene solución general, sin embargo puede obtenerse una solución por medio de un planteamiento en diferencias finitas, con lo cual se sustituye el continuo físico del acuífero por un conjunto equivalente de elementos discretos

El ejemplo simulado consistió en una región de 30 km de largo y 13,5 km de ancho. La recarga considerada fue de 22 mm/año. La conductividad hidráulica k= 5 m/día y la porosidad efectiva promedio del estrato m= 0,10 (10%). Se impuso una cota de agua 35 m, fija en toda la frontera. La discretización topológica consistió en 210 celdas iguales de 1500 m x 1500 m. En la Fig. 5 se presentan las curvas de nivel obtenidos mediante el C1SUBT y el GW8. El coeficiente de correlación fue 0.997583, con un error estándar de correlación de 0,00144.

La solución analítica para la determinación de niveles de depresión debido a un bombeo puntual, en régimen transitorio puede establecerse a partir del método de Jacob (Custodio y Llamas, 1983):

El radio de influencia del pozo puede estimarse a partir de:

t0 es el tiempo de bombeo; ri es el radio de influencia del pozo; r es la distancia al centro del pozo; QB es el caudal de bombeo.

Se simuló el bombeo del acuífero libre, en un medio de transmisibilidad 700 m ² /d, porosidad efectiva del 10% y con un caudal constante en el tiempo de 30 m ³ /h. El nivel estático se propuso en cota 35 m y el piso del acuífero en cota 0 m. Se ensayaron cuatro discretizaciones espaciales diferentes.

analizando el grado de precisión logrado en cada caso al variar las formas, dimensiones y cantidades de las celdas de modelado.

Discretización Ortogonal: Este caso correspondió a 206 celdas cuadradas de 50 m de lado (Fig. 6), constituyendo una superficie de 427500 m ² . La condición de borde impuesta fue cota de agua fija a 35 m. La extracción fue considerada continua de 30 m ³ /h. El pozo fue considerado de radio 0,20 m.

El menor radio para el cual se computaron descensos con tolerable error ( ?0.02 m) fue de 25 m. En esta discretización se evidenció la influencia de las condiciones de borde. La influencia del pozo supera los 500 m, por lo que en la dirección de la menor longitud el cono de depresión es forzado a alcanzar la cota 35 a los 200 m, dando lugar a un cono elipsoidal.

Discretización densificada en celda de bombeo: para mejorar la aproximación a las curvas de Jacob para radios menores 25 m se implementó una variante de discretización en la zona del bombeo. La celda central de 50 m x 50 m fue dividida en cuatro celdas y en la parte central se consideró una celda de 5m x 5m. (Fig. 7). Mediante esta alternativa pudo mejorarse la respuesta del modelo hasta el pozo mismo de bombeo. En efecto considerando una perforación de 0,40 m de diámetro, las máximas diferencias encontradas respecto a Jacob en las curvas correspondientes a r: 0,20 m y 10 m no fueron superiores a 0,02 m (Fig. 8).

Discretización en celdas de tipo sector circular: Se testeó una serie de grillas constituidas por celdas con forma de 1/4 de sector circular para lograr mayor precisión. Se presentan aquí dos alternativas, representativas de todas las analizadas.

Variante con 33 celdas: Se consideró una grilla de 33 celdas concéntricas (Fig. 9) lográndose una aproximación de ?0.02 m a las curvas de Jacob hasta un radio de 100 m. Si bien se evidencia una similitud en los perfiles para gran cantidad de horas de bombeo (mayor a 100 hs.) Respecto a los resultados de la malla ortogonal, al densificarse la discretización en los lugares de mayor variación de niveles, aparecen diferencias en los perfiles para tiempos de bombeo menores. La discretización más intensa es más sensible a los cambios de niveles puesto que tiene menor capacidad de amortiguar los descensos, asimilándose en mejor forma a las hipótesis de partida (pozo sin almacenamiento).

Variante con 37 celdas: Para mejorar los resultados a distancias menores de 10 m del pozo se densificó la grilla colocando 4 nuevas celdas (Fig.10), lográndose que el error en cualquier punto de la curva fuera menor a ?0.02 m. Se consideró un pozo de bombeo de 0,40 m de diámetro. En la Fig. 11 se presentan las evoluciones temporales para r: 0.20 m y 5 m

Se simularon además, bombeos intermitentes de 12 hs con 12 hs de interrupción y de forma idéntica de 6 hs. En todos los casos se trató de un caudal de 720 m ³ /día, por lo que el caudal horario de bombeo intermitente fue de 60 m ³ /h. La celda considerada para el análisis de niveles fue la central, donde se efectuó la extracción de caudales. En la Fig. 12 se ilustran respectivamente la evolución de niveles para bombeo cada 12 hs, cada 6 hs y bombeo continuo.

A modo de análisis comparativo se ha realizado la contrastación de los diferentes perfiles verticales computados por el C1SUBT y por Jacob (Fig. 13 y 14).

En la discretización ortogonal sin densificación en zona de bombeo para r > 25 m se computan cotas de acuífero superiores al método de Jacob no mayores a 0,04 m, por lo que los resultados pueden considerarse aceptables para esos radios, perdiendo precisión para radios menores.

La discretización ortogonal con densificación de grilla computa un error máximo de ?0,03 m en todo el rango de radios, incluyendo el pozo propiamente dicho. Si bien el seguimiento de la curva de Jacob no es riguroso, la aproximación es satisfactoria para cualquier valor de distancia desde el pozo. En esta alternativa con una celda pequeña central de 25 m ² y con escasa capacidad de amortiguamiento de variación de niveles, los descensos calculados son más pronunciados que en el caso anterior. La simulación se asimila de mejor forma al fenómeno real puntual de las cercanías inmediatas del pozo.

La discretización correspondiente a 33 celdas tipo sector circular aproxima de forma satisfactoria la curva de Jacob para r> 10 m, observándose errores no mayores a 0,02 m.

La discretización en 37 celdas tipo sector circular aproxima a la curva de Jacob con un error máximo de ?0,02 m en todos los rangos de distancias desde el pozo.

Es necesario destacar la ventaja que ofrecen las celdas tipo sectores de círculo para este tipo de escurrimientos. Efectivamente, con 37 celdas fue posible la simulación de una zona de más 1200 m de diámetro con una adecuada precisión. La desventaja radica en la poca probabilidad, debido a la forma, de poder ensamblar estas discretizaciones con las que habitualmente se utilizan en modelación de escurrimiento superficial.

Esta aplicación consistió en la estimación de la afectación en términos cuantitativos de un acuífero interfluvio producida por un sistema de riego abastecido por bombeo desde pozos profundos. La zona en estudio se localiza entre los ríos Carcarañá y Paraná en la provincia de Santa Fe, Argentina. La distancia entre los ríos en el área de bombeo varía entre 4,5 a 5 km en una extensión mayor a 15 km (Fig. 15). El actual parcelamiento de la zona comprende lotes de 81 ha (900 x 900 m) de área promedio. La demanda máxima de caudal para riego fue estimada sobre la base de la evapotranspiración promedio y al coeficiente de cultivo, en este caso de maíz, resultando en 5,6 mm/día. La variación temporal de la demanda fue aproximada a una forma trapezoidal con un período máximo de 120 días. La demanda máxima implica en cada lote de 81 ha un caudal de bombeo de 4536 m ³ /día. La situación aquí documentada representa una situación de máxima, en la cual debe regarse todos los días del ciclo cuatrimestral. En función de datos de otros pozos de bombeo explotados en la región, se consideraron en cada lote 4 pozos con un área de riego de 20,25 ha, lo que significó en cada pozo un caudal máximo de bombeo de 47,25 m ³ /h. Debido a la preponderancia del flujo subterráneo en la dirección transversal al sentido de la corriente en los ríos, fue considerada en la simulación una lonja interfluvia de 450 m de ancho (Fig.16). Se adoptaron celdas de 50 m en las lejanías de los pozos, de 22.5 m en las zonas adyacentes y celdas de 5 m en los pozos. La discretización quedó compuesta por 112 celdas. La condición de borde fue una cota promedio anual de 56 m. La recarga fue de 20 mm/año, la conductividad hidráulica k= 5,4 m/día y la porosidad efectiva promedio del estrato m= 0,10.

Primeramente se determinó el estado de equilibrio actual del acuífero bajo situación de no bombeo.

Posteriormente fue simulada una situación de bombeo continuo a través del tiempo comprendido en un período de 15 años. De análisis de la variación de nivel de agua subterránea en la zona central de la lonja (Fig. 17) pudo determinarse que el equilibrio bajo iguales condiciones se logra a los 12 años de bombeo intensivo, produciendo un descenso máximo de 9,26 m al fin del bombeo y un descenso mínimo de 6,91 m al comienzo del bombeo. Asimismo fue simulada una situación de recuperación por recarga natural del acuífero, suponiendo la interrupción del bombeo al cabo de 12 años. El tiempo aproximado de recuperación de su estado anterior al bombeo lleva 12 años.

Esta aplicación ha demostrado la relevante utilidad de la modelación matemática de flujo subterráneo como herramienta para la gestión y control de los recursos hídricos subterráneos.

Se ha desarrollado un modelo para el flujo subterráneo basado en un esquema de celdas, frecuentemente empleado en la modelación de flujos superficiales. El esquema de modelación comprende el planteo de la ecuación de continuidad para cada celda junto a la ecuación de momentum, planteada mediante leyes de vinculación entre celdas basadas en la ecuación de Darcy.

Se ha aplicado el esquema en problemas cuyas soluciones son conocidas, analíticamente ó a través de modelos numéricos alternativos de conocida capacidad, verificándose en todos los casos una satisfactoria aproximación entre los resultados.

El modelo permite una gran flexibilidad en la discretización espacial de las celdas permitiendo optimizar la cantidad de celdas necesarias y el tiempo empleado en la simulación computacional.

Los ejemplos analizados para las contrastaciones, ponen de manifiesto que el modelo propuesto permite dar soluciones tanto a casos de flujos unidimensionales y bidimensionales a gran escala como en situaciones de escala local, con grados de detalle importantes.

El modelo puede ser como herramienta de planificación del uso de los recursos hídricos subterráneos.

Se considera que el esquema propuesto es totalmente apto para ser combinado con esquemas semejantes de otros procesos hidrológicos. La descisión final fue considerar adecuada la inclusión de la estructura de celdas para la simulación de escurrimiento subterráneo en el marco del sistema global de modelación, ya que el grado de precisión, confiabilidad e incertidumbre de los resultados fue similar al obtenido en la simulación de los otros procesos hidrológicos involucrados.

Cunge J; (1975) Two Dimensional Modeling of Flood Plains, Cap.17, Unsteady flow in open channels, Ed. Mahmood K, y Yevjevich V, Water Resources Publications, Fort Collins ,EEUU.pp. 705-762

Custodio E, Llamas M; (1983) Hidrología Subterránea; 2da. Ed., Omega S.A.; Barcelona; España, pp. 670-675

 Mijailov L; (1985) Hidrogeología; Edit. MIR, Moscú, ISBN 5-03-000661-3., pp. 162-166 Naciones Unidas, Area de Recursos Hidráulicos, Departamento de Cooperación Técnica, (1989) Manual de Usuario del Software Ground Water; UN; NYC; EEUU, 193p.

Riccardi G; (1994/1997) Modelos matemáticos incorporados a la modelación integral del sistema suelo-vegetación- atmósfera, Informes de Investigación Consejo de Investigaciones UNR, Rosario, Argentina, 450p. (inéditos).

Riccardi G. (1997) The mathematical modelling of flood propagation for the delimitation of inundation risk zones, Sustainability of Water Resources under Increasing Uncertainty (ed. D. Rosberg et al.) IAHS Pub.Nro 240, ISSN 0144-7815., Wallingford, Inglaterra, pp. 355-363.

Zimmermann E. (1994a) Modelación matemática de flujo subterráneo en medios porosos saturados, Informe Anual 1993-94, CONICET, DH, FCEIA,UNR, Rosario, Argentina, 165p. (inédito).

Zimmermann E (1994b) Evolución Temporal de Niveles Freáticos y de las Zonas de Interacción con la Hidrología Superficial en un Area de Llanura, II Congreso Latinoamericano de Hidrología Subterránea, Santiago, Chile, pp. 253-264..

Zimmermann E. y Riccardi G. (1995) Preliminary Modelling of Fluxes of Water in Flatland Areas, IAHS XX General Assembly of European Geophysical Society (BAHC project), Hamburgo, Alemania.