REPRESENTACIÓN TEMPORAL DE LA RECARGA DE AGUA DE RIEGO EN LA MODELACIÓN RÍO-ACUÍFERO

Javier Pavese * , Leticia Rodríguez . y Carlos Vionnet . *

* Universidad Nacional del Comahue, . Universidad Nacional del Litoral, CONICET

The recharge of a phreatic aquifer involves the transfer of the infiltrated water through the unsaturated zone. Depending on the depth to the water table and the soil characteristics, the estimation of the effective recharge might be a daunting task. However, its magnitude becomes a critical data for those models that simulate the behavior of a phreatic aquifer. Actually, the hypothesis of instantaneous recharge at a constant monthly rate in regional scale modeling of river-aquifer systems produces satisfactory results in many practical cases. Nevertheless, there exist circumstances where the recharge exhibits great temporal variability, as it is the case of the periodical application of irrigation water. Consequently, it is advisable to estimate the real recharge rate and its arriving time to the aquifer. In this work, the preliminary results of a model of the unsaturated zone are presented. Based on a finite difference approximation of the governing equation, the model takes advantage of the similarity between the Richards´ equation and the transport equation. The obtained results is in the process of being incorporated into a model of a stream-aquifer system already implemented in a southern region of Argentina.

La recarga de un acuífero freático involucra la transferencia del agua infiltrada a través de la zona no saturada. Dependiendo de la profundidad de la napa y de las características del suelo, la determinación de la recarga efectiva constituye una tarea compleja. Sin embargo, la magnitud de la misma deviene un dato crítico para los modelos que simulan el comportamiento de un acuífero freático. En realidad, la hipótesis de recarga instantánea a tasa mensual constante en la modelación a escala regional de sistemas río-acuífero, produce resultados satisfactorios en un sinnúmero de casos prácticos. No obstante, existen circunstancias en las que la recarga posee gran variabilidad temporal, como es el caso de la aplicación periódica del riego, por lo que es conveniente estimar la tasa real de recarga y su tiempo de llegada al acuífero. En este trabajo se presentan los resultados preliminares de un modelo de flujo en la zona vadosa. El mismo, basado en una aproximación en diferencias finitas de la ecuación de gobierno, toma ventaja de la similitud entre la ecuación de Richards y la de transporte. Los resultados encontrados serán luego incorporados al modelo de un sistema río-acuífero ya implementado en una región del sur argentino.

La zona vadosa interviene directa o indirectamente en procesos hidrológicos tales como la infiltración, la evaporación, el escurrimiento superficial y la recarga subterránea. Durante las últimas décadas el estudio de la zona no saturada ha experimentado un crecimiento notable. Uno de los factores que motivaron tal interés es el uso de esta capa como depósito de residuos líquidos y sólidos, y la consecuente necesidad de entender su diseminación con el fin de prevenir la degradación del medio ambiente. En realidad, el estudio de la zona vadosa va más allá del problema de la calidad del agua, y radica en el hecho de que es el nexo entre el agua superficial y el agua subterránea. En consecuencia, la inclusión de esta componente en el análisis hidrológico de sistemas de agua superficial y de agua subterránea, en particular en zonas áridas y semiáridas, adquiere una real importancia para una adecuada la planificación, desarrollo y manejo de los recursos hídricos.

Colonia Centenario, ubicada en el valle inferior del Río Neuquén, Pcia. de Neuquén, es una zona semiárida dedicada principalmente al cultivo de frutas bajo riego. El ascenso del nivel freático producido por el efecto combinado del exceso de agua de riego y la recarga proveniente del Río Neuquén que bordea la Colonia, estaría influenciando negativamente la productividad de los frutales. Hasta el presente se han encarado dos estudios de modelación, el primero llevado a cabo por Dufilho y otros (1996), quienes analizaron sólo la componente subterránea del problema.

El segundo trabajo fue desarrollado por Pavese y Rodríguez (1998), quienes modelaron el sistema río-acuífero en su conjunto. En este último caso, y como una primera aproximación al problema, se asumió que la infiltración del agua de riego llegaba al acuífero en forma instantánea y a una tasa mensual constante. Si bien los resultados obtenidos fueron satisfactorios (Rodríguez et al. 1998), se planteó la necesidad de incorporar ciertos mecanismos no tenidos inicialmente en cuenta en el análisis con el fin de mejorar las simulaciones, recreando el esquema de riego periódico programado e incorporando el flujo en la zona no saturada. El objetivo básico perseguido es la obtención de una función de transferencia entre la infiltración del agua de riego y la recarga que responda a las prácticas discontinuas de riego, y que permita estimar más fehacientemente la tasa a la cual el agua llega al sistema subterráneo.

En este trabajo se presentan los resultados preliminares del modelo de flujo en la zona no saturada. Dicho modelo numérico se basa en una aproximación en diferencias finitas de la ecuación de gobierno, tomando ventaja de la similitud entre la ecuación de Richards y la ecuación de transporte. Los resultados que aquí se incluyen serán luego incorporados al modelo río-acuífero ya implementado de manera de contar con una representación integral del sistema que incorpore más acabadamente la física del problema.

Como se mencionara en la Sección 1, al momento se cuenta con código numérico en elementos finitos del sistema río-acuífero en Colonia Centenario. El objetivo general es analizar el efecto combinado de la recarga del exceso de riego y de los caudales del Río Neuquén sobre el ascenso de la napa. La Figura 1 ilustra las variables de entrada utilizadas en las simulaciones, donde se destaca la distribución mensual a tasa constante de la recarga. Por otro lado, la Figura

2 muestra la discretización espacial en elementos finitos triangulares sobre la cual se aplicará el modelo de transferencia (o de flujo vertical) que se describe en las Secciones subsiguientes.

Nótese que la discretización fina en la vecindad del río tiende a capturar la capa límite generada por la interacción río-acuífero (Vionnet y Rodríguez 1998).

Si bien el flujo en la zona vadosa es de carácter tridimensional, no es difícil establecer que el flujo es principalmente vertical. La ecuación universalmente conocida que gobierna el flujo de agua vertical, isotérmico y transitorio en la zona no saturada es la clásica ecuación de Richards.

Dicha ecuación se deduce a partir de la ecuación de continuidad y la ley de Darcy, despreciando los efectos de compresibilidad del suelo. Expresada en función de la presión de succión, dicha ecuación está dada por (Bear 1972)

donde c( h) =  d O (h) / d h. es la capacidad específica de agua en el suelo (L ^-1 ), . es el contenido de humedad volumétrico del suelo (L ³ L -³ ), h es la presión de succión (L), z es la coordenada vertical (positiva hacia abajo) (L), t es el tiempo (T) y K(h) es la conductividad hidraúlica (L/T). En la ecuación (1) intervienen dos parámetros que caracterizan el tipo de suelo, ellos son K h () y c h (). Toda vez que se quiera resolver (1), se debe contar con información acerca de la distribución de estos parámetros. La disponibilidad de valores confiables de K(h) no es tan sencilla ya que el parámetro presenta una gran variabilidad espacial y su medición, de por sí costosa, insume un tiempo considerable.

Durante las últimas décadas, numerosos investigadores se abocaron a la tarea de relacionar la humedad y la conductividad hidráulica no saturada con la tensión (Gardner 1958, Brooks y Corey 1964, y Van Genuchten 1980, entre otros). Es sabido que si se utiliza la fórmula de Gardner

donde s K es la conductividad hidráulica del suelo saturado (cm/día) y a (cm ^-1 ) es un parámetro que caracteriza cada tipo de suelo, la misma no proporciona un buen ajuste con datos experimentales. A pesar de ello, es la expresión adoptada por la mayoría de los modelistas matemáticos debido a su simplicidad.

En el presente ejercicio numérico se adoptó la fórmula de Gardner para la representación de K(h), mientras que para la curva característica de humedad de suelo se eligió una relación también de tipo exponencial

donde las variables ya fueron previamente definidas.

La ecuación de Richards es altamente no lineal y por lo tanto difícil de resolver analítica y numéricamente (ver, por ejemplo, Barry et al. 1993). La no linealidad, representada por la relación funcional de los parámetros K y c con la variable de estado h, puede ser parcialmente eliminada mediante un cambio de la variable dependiente a través de la transformación integral de Kirchoff.

Esta transformación se define mediante la expresión

siendo u(h) la llamada variable de Kirchoff. Introduciendo la definición anterior en la ecuación (1) es posible reescribirla como

Si bien la ecuación (5) mantiene la no linealidad, la misma es más tratable numéricamente.

Esta transformación fue previamente utilizada por Haverkamp and Vauclin (1981), Lomen y Warrick (1978), y Kutilek et. al (1991), entre otros. Introduciendo las expresiones (2) y (3) en (5), y teniendo en cuenta que c( h) =  d O (h) / d h , se obtiene

donde los coeficientes que preceden a las derivadas espaciales son constantes y función de los parámetros físicos del suelo. Además, acorde a la transformación de Kirchoff y a la relación de conductividad hidráulica adoptada, y llamando

la ecuación (6) se reduce a

que no es más que la ecuación lineal de advección-dispersión. La correspondencia de la ecuación de Richards escrita en función del contenido de humedad y la ecuación unidimensional de transporte ha sido analizada por varios autores, entre ellos Yortsos (1987) y Barnes (1989).

En cuanto a la resolución numérica de la ecuación de transporte, existen variadas metodologías cuya descripción escapa al alcance de la presente comunicación. En este trabajo, la aproximación numérica de la ecuación (8) se realizó mediante el método de diferencias finitas mientras que el avance temporal de la solución se obtuvo a través de un esquema explícito. El modelo numérico fue convenientemente validado con soluciones analíticas de la ecuación de transporte. Igualmente se investigó la bondad del modelo para reproducir resultados analíticos obtenidos a partir de soluciones cerradas de la ecuación de Richards, encontrándose una excelente comparación entre valores numéricos y analíticos (Pavese 1999).

Los parámetros adoptados para el cálculo corresponden a las características de los suelos de la zona de Colonia Centenario. Ellos son . r = 0,12; . s = 0,40; Ks = 0,058 cm/s y a = 0,005 cm ^-1 .

Se simuló el flujo en una columna de 150 cm de alto, con una discretización espacial de 25 cm y un intervalo de cálculo de 5 segundos, compatibles con las restricciones de estabilidad impuestas por el esquema explícito, y dirigidos a obtener una adecuada representación espacio-temporal del proceso. Como condición inicial se asumió h = -100 cm en todo el perfil. En la superficie de la columna de suelo se introdujo una condición de borde periódica, consistente en la aplicación de un caudal de 0,1 cm ³ /s durante 500 segundos, alternado don períodos de igual duración con caudal nulo (Figura 3). En la misma figura se graficó la distribución de la presión de succión h en función del tiempo a tres profundidades distintas de la columna de suelo. Por un lado se observa que h alcanza su valor máximo justo antes de que se interrumpa el flujo en superficie, y por otro lado se observa una atenuación de la misma con la profundidad. La declinación presente en la distribución h-t se debe a que el modelo aún no ha alcanzado el régimen de equilibrio periódico impuesto por la condición de borde.

En la Figura 4 se graficó el caudal en superficie conjuntamente con el caudal simulado a una profundidad de 150 cm, profundidad promedio a la cual se encuentra habitualmente la napa freática en Colonia Centenario durante la estación de riego. Se puede observar el desfasaje de los caudales en profundidad a los pulsos de riego en superficie. El máximo Q a 150 cm se produce luego del pico de Q a 0 cm, cuando el perfil de suelo ha drenado lo suficiente.

Si bien los resultados obtenidos hasta ahora con el modelo de advección-dispersión son mayormente de corte académico, los mismos reflejan situaciones que se asemejan a los ciclos de riego. Así, una vez que se incorporen las magnitudes y la variación estacional de los caudales de riego, será posible reemplazar la curva de recarga hasta ahora utilizada (Figura 1), por una similar a la exhibida en la Figura 4. Resultados específicos se mostrarán en la Conferencia.

Teniendo en cuenta el contexto de la modelación río-acuífero en la cual se enmarca este trabajo, se considera que esta forma de estimar la recarga constituye un salto cualitativo respecto de los cálculos tradicionales basados en un balance hídrico de paso mensual. No obstante, la discretización temporal fina de la recarga generará nuevas dificultades numéricas debido a la aparición de una nueva escala temporal superpuesta a las ya analizadas (Rodríguez et al., 1998, Vionnet y Rodríguez, 1998)

Una combinación de subsidios provenientes de diferentes fuentes hicieron posible el desarrollo del presente trabajo. En particular se agradece a la Universidad Nacional del Comahue y a la Universidad Nacional del Litoral por la ayuda recibida, así como el soporte del CONICET a través del subsidio PIP 0616/98.

8.

Barnes, C.J., 1986. Equivalent formulations for solute and water movement in soils. Water Resources Research, 22(6): 913-918.

Barry, D.A., Parlange, J.Y., Sander, G.C., and Sivapalan, M. 1993.A class of exact solutions for Richards’ equation . Journal of Hydrology, 142, 29-46.

Bear, J., 1972. Dynamics of fluids in porous media. Dover Publications Inc., New York. 764pp.

Brooks, R.H. and Corey, A.T., 1964. Hydraulic properties of porous media. Colorado State University Hydrology Paper No. 3, March, 27 p.

Dufilho, A.C., Stangaferro, S.R., Horne, F. y Olivera, D.A., 1996. Aplicación de un modelo matemático y un SIG para simular el flujo subterráneo en el valle superior del Río Negro.

Trabajo No. 96. XVI Congreso Nacional del Agua. San Martín de los Andes, Neuquén, Noviembre 1996.

Gardner, W.R., 1958. Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Science, Vol. 85, No. 4, pp. 228-232.

Haverkamp, R., and Vauclin M., 1981. A comparative study of three forms of the Richards´ equation used for predicting one-dimensional infiltration in unsaturated soil. Soil Sci. Soc. Am. J., 45, 13-20.

Kutilek, M., Zayani, K., Haverkamp, R., Parlange, J.Y., and Vachaud, G., 1991. Scaling of the

Richards equation under invariant flux boundary conditions. Water Resources Research, Vol. 27, No.9, p. 2181-2185.

Lomen, D.O., and Warrick A.W., 1978. Time-dependent solutions t oteh one-dimensional linearized moisture flow equation with water extraction. Journal of Hydrology, 39, p.59-67.

Pavese, J., 1999. Tesis de Maestría en Ingeniería de los Recursos Hídricos. Facultad de Ingeniería y Ciencias Hídricas, Universidad Nacional del Litoral, Santa Fe (en preparación).

Pavese, J., y Rodríguez, L., 1998. Análisis de la influencia del riego sobre el ascenso del nivel freático. XVII Congreso Nacional del Agua. Santa Fe, 3-7 agosto 1998. p. 90-98.

Rodríguez, L., Pavese, J., and Vionnet, C., 1998. The effect of excess irrigation water on the water table fluctuations in the Neuquen River Valley, Argentina. Hydrology in a changing environment. Vol. II, pp. 163-164. Proceedings of the British Hydrological Society International Conference, Exeter, U.K., 6-10 July.

Vionnet C. and Rodríguez L., 1998. A simple time-stepping strategy for coupled river-aquifer systems. In: Computational Mechanics, New trends and applications. Editors: S. Idelsohon, E. Oñate and E. Dvorkin. Proceedings of the IV World Conference on Computational Mechanics. Buenos Aires, Argentina.

Yortsos, Y.C., 1987. The relationship between immisible and miscible displacement in porous media. Am. Inst. Chem. Eng. J., 33(11): 1912-1915.