Luis Guarracino (1) , Lucila Candela Lledó (2) , Juan E. Santos (1,3) .

(1) CONICET, Departamento de Geofísica Aplicada, Facultad de Cs. Astronómicas y Geofísicas, Universidad Nacional de La Plata (Argentina).

(2) Departamento de Ingeniería del Terreno, Universidad Politécnica de Cataluña (España).

(3) Department of Mathematics, Purdue University (USA).

Monitoring of field experiments using tracers and numerical simulation are both indispensable tools for the study of water and contaminant transport in the unsaturated zone. This paper presents a numerical model for analyzing and interpreting field experiments which were carried out in an experimental patch over the Maresme aquifer (Barcelona, Spain). Tests were designed to characterize the migration of glyphosate herbicide in the unsaturated zone and to determine the risk of contamination in the aquifer. The munerical model is one-dimentional and it is based on the simultaneous resolution of the Richards’ equation and the advection-diffusion equation by means of mixed finite element techniques. The model icludes linear adsorption proceses and first order decay of the contaminant substance. Model’s parameters and boundary condictions were defined using laboratory and field data.

El monitoreo de ensayos con trazadores y la simulación numérica constituyen herramientas imprescindibles para el estudio del transporte de agua y de sustancias contaminantes en la zona no saturada. En este trabajo se presenta un modelado numérico para analizar e interpretar ensayos realizados en una parcela experimental ubicada sobre el acuífero del Maresme (Barcelona, España). Los experimentos fueron diseñados para caracterizar la migración del herbicida glifosato en la zona no saturada y determinar el riesgo de contaminación en el acuífero. El modelo numérico utilizado es unidimensional y se basa en la resolución simultánea de las ecuaciones de Richards y de transporte advectivo-difusivo mediante técnicas de elementos finitos mixtos. El modelado incluye procesos de adsorción lineal y de degradación de primer orden de la sustancia contaminante. Los distintos parámetros y las condiciones de borde se definieron en base a los datos medidos en la parcela y en laboratorio.

Las prácticas agrícolas constituyen una importante fuente de contaminación compuesta por sustancias químicamente complejas susceptibles de alcanzar la zona saturada del terreno.

Los pesticidas y fertilizantes que son utilizados en forma periódica sobre grandes extensiones areales hacen altamente vulnerables a la contaminación a todos aquellos acuíferos cuya superficie de impluvio está sometida a la explotación agrícola.

El potencial contaminante de un plaguicida está dado por su movilidad y persistencia. El plaguicida debe ser lo suficientemente móvil y persistente como para alcanzar y eliminar el organismo específicamente atacado minimizando su impacto sobre las aguas subterráneas. Estas dos características están controladas principalmente por los procesos de adsorción y de degradación.

La adsorción se manifiesta como un retardo del movimiento del contaminante respecto de la velocidad del agua. Este proceso no afecta la cantidad total de plaguicida presente en el suelo pero puede disminuir e incluso eliminar la cantidad disponible para el transporte. La magnitud de este fenómeno dependerá de la composición físico-química del suelo, en particular de su contenido en materia orgánica y arcillas. Por otra parte, los plaguicidas pueden sufrir la rotura de su estructura molecular en fragmentos que dan lugar a compuestos inorgánicos como productos finales de la reacción. Este proceso se denomina degradación, siendo ésta la única vía para que un plaguicida sea totalmente eliminado del medio ambiente.

El objetivo de este trabajo es presentar un modelo numérico que describa los procesos anteriormente enunciados y realizar una simulación numérica para analizar e interpretar ensayos realizados en una parcela experimental. Los ensayos fueron diseñados para determinar el riesgo de contaminación del acuífero del Maresme (Barcelona) por el uso del herbicida glifosato.

El glifosato es uno de los herbicidas más ampliamente usados a nivel mundial. Se trata de un herbicida no selectivo de post emergencia utilizado para controlar una amplia variedad de hierbas y arbustos tanto en zonas no cultivadas como en distintos tipos de cultivos. El uso del glifosato, en condiciones normales, es considerado como de bajo riesgo para la salud humana.

Sin embargo su único metabolito conocido, el AMPA (aminomethylphosponic acid), resulta más persistente y tóxico que el producto original. Restos de AMPA han sido detectados en acuíferos someros de Holanda y Alemania.

El acuífero del Maresme es una franja costera que se extiende al norte de la ciudad de Barcelona, limitada hacia el interior por la Cordillera Costero Catalana. Prácticamente toda el área de recarga del acuífero esta ocupada por actividades agrícolas de regadío detectándose una elevada contaminación por nitratos. Curiosamente, y a pesar de su uso intenso, se han detectado muy pocos restos de plaguicidas en el acuífero. Este último punto abrió una serie de cuestionamientos acerca de la dinámica de los plaguicidas en la zona no saturada (ZNS) y de los metabolitos que deberían analizarse. En este sentido, durante los años 1994 y 1995 se realizaron ensayos en dos parcelas experimentales no cultivadas para estudiar la migración del herbicida glifosato en la ZNS. En el presente trabajo se analizan los datos obtenidos en una de estas parcelas en la que se realizó un ensayo de 100 días de duración que consistió en la aplicación simultánea de glifosato y de un trazador conservativo.

Desde el punto de vista matemático el problema consiste en la resolución de un sistema de ecuaciones diferenciales acopladas que describen el movimiento del agua y del contaminante en un medio poroso de saturación variable. Se asume que el movimiento del agua en la ZNS obedece la ecuación de Richards (Richards, 1931) y que la concentración del soluto (contaminante) está gobernada por la ecuación de transporte advectivo-difusivo. Para describir los procesos de adsorción se utilizaron isotermas lineales mientras que la degradación de la sustancia contaminante se modeló mediante un término de primer orden en la ecuación de transporte advectivo-difusivo.

La parcela no cultivada, de 168 m ² de superficie, fue instrumentada en los terrenos de un centro de investigación agraria de Cabrils (Barcelona). La litología de la zona es básicamente arenosa y homogénea con escasa fracción de materia orgánica y arcillas, encontrándose el nivel freático a 5.5 m de profundidad.

El experimento de campo consistió en la aplicación simultánea del herbicida glifosato y de un trazador conservativo (Br) sobre la superficie de la parcela. Una vez realizada la aplicación se procedió a un monitoreo de la migración de ambas sustancias mediante un muestreo periódico.

Para el ensayo de trazador conservativo se utilizó una disolución de NaBr preparada con 6000 gr del producto comercial disueltos en 300 litros de agua. Dicha disolución equivale a una concentración de Br de 15940 mgr/l. En cuanto a la estimación del volumen óptimo de aplicación del glifosato se basó en las recomendaciones del fabricante en cuanto a dosis y volumen de caldo a utilizar, resultando una concentración de 7200 mgr/l. En ambos casos la aplicación se realizó mediante pulverización standard sobre la superficie de la parcela.

La presión del agua en la zona no saturada se midió in situ mediante tensiómetros instalados a 30, 60, 90 y 120 cm de profundidad. El contenido de humedad y las concentraciones en el suelo se analizaron mediante ensayos destructivos a distintas profundidades. En cuanto a las variables meteorológicas se dispuso de valores diarios de temperatura, humedad relativa, velocidad media del viento, radiación neta y precipitaciones, proporcionados por una estación situada a 200 m del área experimental.

Las muestras inalteradas de suelo se tomaron mediante barrena con martillo percutor en dos puntos diferentes a los 14, 23, 34, 62 y 92 días de iniciados los ensayos. Previo a la realización del experimento se determinaron las concentraciones de Br y glifosato presentes en el suelo con el fin de establecer el estado inicial de contaminación.

En esta sección se describen las ecuaciones que gobiernan la migración de una sustancia contaminante sometida a procesos de adsorción y degradación, junto con los métodos numéricos adoptados para su resolución. El modelo propuesto es unidimensional pues se asume que el flujo en la ZNS es exclusivamente vertical.

Para modelar el flujo de agua se utilizó la ecuación de Richards expresada mediante el siguiente sistema:

donde h es la altura piezométrica (cm); . , el contenido de agua (cm ³ /cm ³ ); q , el flujo de agua (cm/s); K , la conductividad hidráulica (cm/s); S , un sumidero que representa el agua absorbida por las raíces (cm ³ /cm ³ s); z , la coordenada vertical (cm); y t , la variable temporal (s).

Para resolver numéricamente la ecuación (1) es necesario contar con expresiones del contenido de agua y de la conductividad hidráulica en función h . En este modelado se utilizaron las curvas propuestas por van Genutchen (1980):

donde m n =-1 1/; . r y . s son los contenidos de agua residual y máximo; K s , la conductividad hidráulica cuando el medio está saturado; y finalmente, a y n son parámetros del modelo que se determinan experimentalmente.

La condición de borde en la superficie del terreno es de tipo Neumann y está asociada a los valores de la precipitación (o irrigación) y de la evaporación de la humedad del suelo. En el extremo inferior del dominio la condición de borde es de tipo Dirichlet y corresponde a la altura piezométrica calculada a partir de los valores de los niveles freáticos.

Para modelar la concentración de una sustancia contaminante sometida a procesos de adsorción y de degradación de primer orden se utilizó la ecuación de transporte advectivo-difusivo expresada mediante el siguiente sistema:

donde la variable c es la concentración volumétrica del soluto (gr/cm ³ ); R, el factor de retardación (adimensional); D, el coeficiente de dispersión-difusión (cm ² /s); u , el flujo advectivo-difusivo del contaminante (gr/cm ² s); y ｵ, el coeficiente de degradación de primer orden.

Si suponemos que la concentración del soluto en las fases sólida (matriz porosa) y líquida (agua) están relacionadas por una isoterma de adsorción lineal, el factor de retardación R y el coeficiente de degradación ｵ, toman las siguientes expresiones (Jury et al.,1991):

siendo . la densidad de la matriz porosa (gr/cm ³ ); K D , el coeficiente de distribución o de reparto (cm ³ /gr) que caracteriza a la isoterma lineal; ｵW y ｵS , las constantes de degradación del contaminante en las fases líquida y sólida (1/s) , respectivamente. El coeficiente de dispersión D se calcula mediante la siguiente expresión (Bear,1972):

donde . es la dispersividad del medio (cm); D 0 es el coeficiente de dispersión iónica o molecular en el agua (cm ² /s) y t es el factor de tortuosidad (adimensional). Para evaluar dicho factor se utiliza el modelo propuesto por Millington y Quirk (1961):

Para la resolución del sistema (2), se utilizaron condiciones de borde de tipo Neumann en ambos extremos del dominio correspondientes al valor del flujo advectivo-dispersivo del contaminante. Por otra parte, es necesario establecer un perfil de concentración que describa el estado de contaminación inicial a partir del cual se comienza la simulación.

Los sistemas (1) y (2) deben resolverse en forma secuencial. En primer lugar se determina el campo de velocidades y el contenido de agua mediante la resolución de la ecuación de Richards. A partir de estos valores es posible evaluar los coeficientes de la ecuación de transporte y proceder a su aproximación numérica.

La discretización temporal del sistema (1) se realizó mediante un esquema backward Euler combinado con un método de Picard modificado para tratar las no-linealidades de la ecuación, según las ideas presentadas en el trabajo de Celia et al. (1990). Para la discretización espacial se utilizó un método mixto híbrido de elementos finitos que permite aproximar en forma simultánea y con igual precisión tanto la altura piezométrica como el flujo de agua (Guarracino and Santos, 1998). Para la discretización del sistema (2), nuevamente se utilizó un esquema backward Euler junto con un método mixto híbrido de elementos finitos para aproximar la concentración y el flujo advectivo-dispersivo del contaminante.

El algoritmo resultante fue escrito en lenguaje FORTRAN 77. Las resultados numéricos de la ecuación de Richards fueron testeados usando el algoritmo desarrollado por Celia et al. (1990), en tanto que los resultados de la ecuación de transporte fueron validados utilizando soluciones analíticas conocidas (Guarracino et al., 1998).

Para calcular los parámetros de la curva de retención del modelo de van Genutchen, se utilizó el programa de dominio público denominado RETC (van Genutchen, 1991). Dicho programa estima los coeficientes del modelo a partir de la curva observada mediante un método de optimización paramétrica no lineal por mínimos cuadrados. La curva de retención se obtuvo a partir de la correlación de los valores de la altura piezométrica y del contenido de agua a los 30, 60 y 90 cm de profundidad. En la tabla 1 se listan los parámetros obtenidos y en la figura 1 se muestra el ajuste alcanzado.

Tabla 1: Parámetros de la curva de retención.

El coeficiente de reparto K D se define como el cociente de las concentraciones en las fases sólida y líquida. Para su determinación se utilizó la técnica de laboratorio denominada experiencia en batch. Estos ensayos se llevan a cabo sobre suspenciones de suelo con varias disoluciones del herbicida a diferentes concentraciones. Las muestras son continuamente agitadas a temperatura constante. Luego de cierto tiempo, se separa por centrifugación la fase líquida de la sólida y finalmente se determina la concentración en equilibrio del herbicida. En la tabla 2 se listan los valores obtenidos.

Tabla 2: Coeficientes de reparto determinados mediante ensayos de batch.

Para el cálculo de las constantes de degradación ｵW y ｵS se utilizó la siguiente expresión (Sparks, 1989):

donde t 1 /2 es la vida-mitad (half-life) del contaminante. El parámetro t 1 /2 expresa el tiempo necesario para que la masa inicial del contaminante aplicado al subsuelo se reduzca a la mitad.

En el caso del glifosato la vida-mitad se estima en 47 días. (Candela et al., 1998). El resto de los parámetros utilizados se listan en la tabla 3.

Tabla 3: Parámetros del modelo.

El objetivo del trazador de Br es caracterizar el flujo de agua en la ZNS. Cuando se considera una sustancia químicamente no reactiva como el Br, la ecuación (2) se simplifica. Al no producirse degradación de la sustancia contaminante ni procesos de adsorción, el coeficiente ｵ es nulo y el factor de retardación R resulta constante e igual a 1.

Para realizar la simulación numérica es necesario definir las condiciones de borde en base a la información disponible y establecer un estado inicial a partir del cual se comienza la simulación. Esto debe realizarse tanto para la ecuación de Richards como para la ecuación de transporte advectivo-difusivo.

La condición de borde superior para resolver la ecuación de Richards corresponde a un flujo de agua que ingresa o egresa del sistema. Los valores negativos del flujo indican infiltración de agua al subsuelo, en tanto que los positivos expresan evaporación desde la superficie. Los datos meteorológicos fueron utilizados para calcular la evapotranspiración potencial diaria mediante la fórmula de Penman. A partir de estos valores se modeló el efecto de las raíces (termino S en el sistema (1)) y la evaporación superficial. En el borde inferior se especificó la altura piezométrica calculada a partir de la posición del nivel freático. En cuanto a la condición inicial se supuso que el perfil de humedad se encontraba en equilibrio.

Figura 2: perfiles de concentración de Br luego de 14 (a), 34 (b), 62(c) y 92 (d) días de simulación.

Para resolver la ecuación de transporte se especificó el valor del flujo de la disolución ingresada al subsuelo desde la superficie, en tanto que en el extremo inferior del dominio se supuso que no se produce entrada ni salida de Br ( u =0). Las concentraciones medidas antes de iniciar el ensayo no superan los 0.2 microgr/gr por lo que se asumió un estado de contaminación inicial nulo.

En la figura (2) se muestran los perfiles de concentración de Br luego de 14, 34, 62 y 92 días de simulación. En las gráficas también se ilustran los valores de concentración observados a una misma profundidad en dos zonas distintas de la misma parcela. Es importante destacar que la dispersión de los datos puede tener distinto origen: presencia de inhomogeneidades en el terreno, distribución no uniforme de las precipitaciones y/o de la pulverización del Br sobre la superficie de la parcela, errores de medición, etc. Si se tiene en cuenta que estos factores no son en general factibles de cuantificar y de modelar, el ajuste obtenido puede considerarse como muy bueno.

Las condiciones iniciales y de borde para simular el transporte de glifosato son idénticas a las del caso anterior, con la salvedad de que deberá especificarse el valor del flujo de glifosato ingresado al subsuelo en el extremo superior del dominio.

En la figura 3 se muestran los perfiles de concentración correspondientes a los días en que se realizaron los muestreos. En ellos se observa que el herbicida es fuertemente adsorbido en los primeros 10 cm de suelo y que más allá de esta profundidad la concentración es prácticamente nula. La comparación del avance de los frentes de concentración del glifosato y del trazador de Br ponen de manifiesto la magnitud de este fenómeno (a los 92 días de simulación el Br puede ser encontrado hasta una profundidad de 3 m). La degradación del glifosato se manifiesta como una continua disminución de la amplitud de los perfiles de concentración en los sucesivos días de muestreo.

En la figura 4: se muestra la comparación de los perfiles simulados con los valores observados. El dato observacional más superficial corresponde a la profundidad de 10 cm que prácticamente coincide con el límite predicho por la simulación numérica para la contaminación por glifosato. Este hecho explicaría los bajos valores de concentración observados.

En la figura 5 se ilustra la evolución temporal de la masa de glifosato. El porcentaje de masa recuperada, que se calcula a partir de los perfiles de concentración y de contenido de agua, puede ser considerado como un indicador de la representatividad de los datos. Debido a que el glifosato sufre procesos de degradación la masa disminuye con el transcurso del tiempo. Para realizar la simulación se consideró una vida-mitad de 47 días, lo que implica que transcurrido dicho tiempo la masa inicial se reduce a la mitad. Este hecho es predicho en forma exacta por el modelo según se observa en la figura. La masa recuperada calculada a partir de los datos es sensiblemente menor lo que pone nuevamente de manifiesto los bajos valores de concentración muestreados.

De acuerdo con la simulación numérica los valores máximos de concentración de glifosato se encuentran a profundidades menores de 10 cm. Si se acepta como válido este resultado, las mayores concentraciones del perfil no fueron muestreadas lo que explicaría los bajos valores de concentración medidos y como consecuencia la poca masa recuperada. Por otra parte, la simulación numérica no predice la contaminación encontrada a profundidades mayores. Este fenómeno posiblemente se deba a efectos de macrodispersión que no fueron contemplados en el modelo.

De acuerdo con los resultados obtenidos, la simulación numérica resulta una herramienta válida para el estudio del transporte de sustancias contaminantes en la ZNS. En este sentido, el modelado ha permitido interpretar los resultados experimentales y establecer algunas hipótesis sobre la dinámica del herbicida analizado.

En cuanto a la migración del glifosato se puede concluir que resulta fuertemente adsorbido en los primeros centímetros del suelo y que la contaminación en el acuífero solo podría tener lugar por efectos de macrodispersión (direcciones preferenciales de flujo, macroporosidad).

Finalmente, es de destacar que el modelo propuesto resulta de utilidad para definir futuras campañas de muestreo pues permite definir a priori algunas pautas sobre la evolución de la contaminación en la ZNS.

Bear, J., 1972. Dynamics of Fluids in Porous Media. Esevier, New York.

Candela, L., S. Rao, M. Margiotta and A. Reboucas, 1998. Soil and Groundwater Pollution from Agricultural Activities. IHP-V, Technical Documents in Hydrology, Nro. 19, UNESCO, Paris.

Celia, M. A., E. T. Bouloutas and R. L. Zarba, 1990. A General Mass-conservative Numerical Solution for the unsaturated flow equation.En Water Resources Reserch, 26, 1483 - 1496.

Guarracino, L, C. L. Ravazzoli y J. E. Santos, 1998. Modelado de transporte de contaminantes en las zonas saturada y no saturada utilizando métodos de elementos finitos. Memorias del 4to. Congreso Latinoamericano de Hidrología Subterránea, Volumen 1, pp 156-168.

Guarracino, L. and J. E. Santos, 1998. Numerical Modelling of Unsaturated Flow Using a Hybridized Mixed Finite Element Procedure. En CD-ROM Proceedings of the Fourth World Congress on Computational Mechanics, Part IV, Section 7.1, pp. 1-10.

Jury, W. A., W. R. Gardner and W. H. Gardner, 1991. Soil Physics, Jhon Wiley & Sons Inc., New York.

Richards, L. A., 1931. Capillary conduction of liquids through porous mediums. En Physics 1: 318- 333.

Sparks, D. L., 1989. Kinetics of Soil Chemical Processes, Academic Press, Inc., California.

van Genutchen, M. T., 1980. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. En Soil Science Society of America Journal, 44, 892-898.