INFLUENCIA DE LAS CURVAS DE RETENCIÓN DE HUMEDAD EN LA ESTIMACIÓN DE LA EVAPORACIÓN DE ACUÍFEROS SOMEROS
El objetivo de este artículo es determinar la influencia que tienen las curvas características del suelo no saturado en la determinación de la evaporación en una columna de suelo desnudo. Se utiliza el código VS2DTI (Lappala et al., 1987) para determinar la evaporación en función de la profundidad del agua subterránea y se efectúa un análisis de sensibilidad de los distintos parámetros de las curvas características. Adicionalmente se compara el efecto de usar la curva h(θ) obtenida en laboratorio con las propuestas a partir de la curva granulométrica.
MARCO TEÓRICO
Ecuación del flujo en la zona no saturada
La ecuación que rige el flujo en la zona no saturada (Richards, 1931), se obtiene combinando la ley de Darcy con la ecuación de continuidad, que para un flujo vertical en una dimensión se escribe como:
Donde C(h) es la capacidad específica de humedad (L-1); K es la conductividad hidráulica (LT-1); z es la coordenada vertical (L), positiva hacia abajo; h es la altura de presión (L), función del contenido volumétrico de agua θ (L3L-3). Para obtener la solución de la ecuación de flujo vertical en una dimensión, se requiere especificar las condiciones iniciales y de borde apropiadas. La condición inicial está dada por la altura de presión h en todo el perfil y las condiciones de borde que pueden ser aplicadas son la condición de potencial o condición de flujo.
Solución numérica de la ecuación de Richards
Existen diversos códigos que permiten resolver la ecuación del flujo en la zona no saturada. La discretización espacial y temporal del dominio se realiza a través de esquemas de diferencias finitas o elementos finitos. Mediante diferencias finitas es posible obtener soluciones adecuadas para la mayoría de los problemas, y presenta limitaciones cuando la geometría del dominio es compleja y cuando las propiedades de los elementos no tienen un comportamiento homogéneo en todas sus direcciones. El método de los elementos finitos, permite ajustarse a geometrías de todo tipo, para lo cual requiere implementar métodos de discretización adecuados.
El código VS2DTI (Variably Saturated 2 Dimensional Transport Interface) resuelve la ecuación de Richards para el flujo en la zona no saturada (Healy, 1990), utilizando un esquema de diferencias finitas que permite analizar casos de una o dos dimensiones con coordenadas cartesianas o radiales. Las curvas de succión y de conductividad hidráulica se representan mediante las relaciones propuestas por van Genuchten (1980), Brooks y Corey (1964), Haverkamp et al. (1977) o puntos experimentales. El código permite establecer condiciones de borde de contorno sobre la carga hidráulica (H), presión (h), flujo, infiltración, evaporación, transpiración de plantas y superficies de filtración.
La condición de evaporación se puede imponer en dos etapas. La primera etapa, se produce cuando el suelo cercano a la superficie está húmedo y el agua deja el sistema a una tasa equivalente a la demanda evaporativa de la atmósfera, conocida como tasa potencial. Esta tasa se mantiene en la medida en que el suelo sea capaz de satisfacer dicha demanda. Cuando el suelo se va secando se da inicio a la segunda etapa en que la tasa evaporativa cae rápidamente debido a la disminución del valor de la conductividad hidráulica no saturada (Tyler et al., 2006).
El flujo evaporativo (EVPa) se establece debido al gradiente de presión existente entre la atmósfera (celda imaginaria) y la primera celda del dominio representado por la celda 1, tal como se presenta en la Figura 1.(Lappala et al., 1987; Liu et al., 2005) y se escribe como:
donde EVPa es la evaporación real en la superficie [LT-1]; EVP(t) es la evaporación potencial del suelo [LT-1]; Ks es la conductividad hidráulica saturada [LT-1]; Kr es la conductividad hidráulica relativa (-); SRES es la resistencia de la superficie que corresponde al recíproco de la distancia desde el nodo hasta la superficie [L-1]; HA es el potencial de presión de la atmósfera [L]; y h es el potencial de presión en el primer nodo de la superficie del volumen modelado [L].
Figura 1. Esquema para la condición de borde de evaporación
El parámetro SRES se utiliza para modelar la existencia de una costra superficial menos permeable que el suelo (e.g. costra salina). HA es calculado mediante la ecuación de Kelvin (Lappala et al., 1987) como:
donde R es la constante de los gases (8.31 g m2 s−2 K−1 mol−1); T es la temperatura del aire (K); Mw es el peso molecular del agua (0,018 Kg mol−1); g es la aceleración de gravedad (9,81 m s−2); y Rh es la humedad relativa de la atmósfera (-).
METODOLOGÍA
Para lograr los objetivos propuestos se estimaron las curvas características a partir de ensayos de laboratorio y características granulométricas de una muestra de arena fina. La evaporación se determinó mediante el código numérico VS2DTI.
Estimación de Curvas Características
Las curvas características hidrodinámicas del suelo se determinaron por dos métodos: 1) a partir de ensayos de laboratorio, y 2) a partir de modelos empíricos o funciones de pedotransferencia (PFT). En el laboratorio se determinó la curva de succión h(θ) por los métodos de la columna colgante (ASTM, 2003) y de la olla de presión (ASTM, 2003) que consistieron en la aplicación de diferentes succiones a una muestra de suelo y el registro del contenido de humedad gravimétrico. Primero se realizó el método de la columna colgante, saturando la muestra con agua y aplicando succiones entre 0 y 80 kpa, manteniendo cada succión durante un día y luego se trasladó la muestra a la olla de presión y se aplicó succiones entre 80 y 1500 kpa. Posteriormente se secaron las muestras en un horno a 110+/-5°C. La conductividad hidráulica saturada Ks se midió con un permeámetro de carga constante, y se determinó la conductividad hidráulica no saturada K(h) mediante el ajuste de los parámetros de las ecuaciones de van Genuchten (1980) , que se escriben como:
Donde Se es la saturación efectiva [-]; θh, θr y θs son el contenido de humedad volumétrico, residual y saturado respectivamente [L3L-3]; Ks es la conductividad hidráulica saturada [LT-1]; m y n son constantes [-]; y m = 1-1/n.
La estimación de las curvas características por medio de las PFT se realizó con el modelo ROSETTA (Schaap et al., 2001) que estima los parámetros del ajuste propuesto por van Genuchten (1980), utilizando redes neuronales (Schaap and Leij, 1998) mediante 5 modelos jerárquicos que dependen de los datos de entrada (clasificación del suelo USDA, textura, densidad del suelo, capacidad de campo o contenido de humedad a 330 cm de succión y punto de marchitez o contenido de humedad a presión de succión15 bar). En este trabajo se usaron los modelos denominados H2 y H3. El modelo H2 utiliza el porcentaje de arena, limo y arcilla como datos de entrada, y el modelo H3 incorpora, además de la textura, la densidad de suelo seco como predictor. Las características granulométricas del suelo fueron determinadas por los ensayos de tamizado y sedimentación.
Discretización espacial y condiciones de contorno
Para la simulación del flujo evaporativo se utilizó el código VS2DTI con una discretización en coordenadas radiales para el caso de una columna circular de 15.3 cm de diámetro y 150 cm de altura, y un refinamiento variable con celdas más finas en la zona superior por donde se produce el flujo evaporativo. En total se utilizaron 230 celdas. La celda superior posee un espesor de 10−4 cm y las siguientes crecen en progresión geométrica a una razón de 5% más que el grosor de la anterior.
Como condición de borde superior del dominio se utilizó una condición de evaporación que permite fijar un valor para la evaporación potencial (EVP) y para el potencial de presión de la atmósfera (HA). Se indicó que no había costra en la superficie que limitara el flujo evaporativo y por lo tanto la variable SRES se tomó como 2/Dz, considerando el valor de Dz de la celda superior. Para el borde inferior del dominio se fijó una condición de borde de carga constante (H = cte) para definir la profundidad del nivel freático, que corresponde al caso en que la presión de succión del suelo es nula (h = 0). Para los bordes laterales se impuso una condición de borde impermeable.
RESULTADOS Y DISCUSIÓN
Estimación de las curvas características h(θ) y K(θ) El modelo fue aplicado para estimar la evaporación desde el agua subterránea en una columna de suelo cuya granulometría, densidad, contenido de humedad y conductividad hidráulica son conocidos. En la Figura 2 se presenta la curva granulométrica, que corresponde a una curva típica de arena con material fino menor al 10%, y de acuerdo con la clasificación de la USDA (United States Department of Agriculture) el suelo estudiado tiene un 95% de arena, 2% de limo, y 3% de arcilla.
Los puntos de la curva de succión determinados experimentalmente fueron adecuados para describir su forma en el rango de presiones de succión entre 0 y -10 m.c.a. El método de la columna colgante resulta efectivo para definir el primer tramo de la curva (0 a 1 m.c.a.) al igual que resulta efectivo el método de la olla de presión para succiones más elevadas. Es recomendable que ambos métodos se apliquen en un rango intermedio de presiones para corroborar que los valores medidos por ambos métodos concuerden.
En la Tabla 1 se presentan los parámetros de van Genuchten de las ecuaciones (4) y (5) para h(θ) y K(h). La curva de succión ajustada a los datos experimentales de la muestra de arena fina es capaz de mantener sus poros con un contenido de humedad cercano a la saturación (q » qs) en una franja capilar amplia y cercana a 60 cm. Para una profundidad de la napa inferior a 60 cm el perfil de humedad se mantiene saturado y el flujo evaporativo estará dominado por la capacidad de la atmósfera para inducirlo. El suelo será capaz entonces, de acuerdo al código VS2DTI, de proporcionar en todo momento el flujo demandado dentro de este rango de profundidad del nivel freático.
En la Figura 3 se presentan las curvas de succión h(q) y de conductividad hidráulica K(q) ajustadas con van Genuchten (1980). Las curvas h(θ) estimadas a partir de la granulometría muestran claras diferencias principalmente en la zona cercana a la saturación (Figura 3a). El contenido de humedad saturado (qs) obtenido mediante los modelos H2 y H3 resultaron mayores al medido en laboratorio. La curva medida en laboratorio muestra una franja capilar que se extiende hasta una succión cercana a los 50 cm, mientras que para los ajustes mediante los modelos H2 y H3 esta franja no sobrepasa los 30 cm. El contenido de humedad residual (qr) obtenido con los modelos H2 y H3 es muy similar al medido en laboratorio.
En la Figura 3b se observa que las diferencias entre la curva ajustada a partir del valor de Ks medido en el laboratorio y las estimadas mediante los modelos H2 y H3 a partir de la granulometría son importantes. El dominio para el cual está definida la curva K(q) es diferente, debido a que el extremo superior (qs) es diferente para cada uno de los ajustes. Estas diferencias son mayores cerca de la condición de saturación y van disminuyendo hacia el contenido de humedad residual (q » qr). La pendiente de la curva K(q) asociada a cada ajuste se ve también modificada debido a estas diferencias ya mencionadas, pero además por la variación del parámetro n en el ajuste de van Genuchten que afecta el valor de m en la ecuación (5).
Las mayores diferencias se encuentran en los parámetros p, n y Ks. Los modelos H2 y H3 sobreestiman los valores de p y Ks, y subestiman los valores de n. Una sobrestimación del parámetro p produce una subestimación de la presión de burbujeo o de entrada de aire (hb =1/p) y por lo tanto una disminución en la altura de la franja capilar, mientras que un aumento (disminución) en el parámetro n produce un aumento (disminución) en la altura de la franja capilar (Figura 4). Las diferencias existentes entre las estimaciones de las curvas de succión y Conductividad hidráulica experimental y las curvas ajustadas a partir de la granulometría se explican principalmente por la metodología empleada en los modelos empíricos. El método de redes neuronales utilizado por Schaap y Leij (1998) depende fuertemente de los datos utilizados en el proceso de entrenamiento de dichas redes, principalmente de la cantidad y distribución textural de los datos, y a pesar de que el entrenamiento para este modelo se realizó con cerca de 10.000 puntos h(θ) y 1500 puntos K(h), la incertidumbre reside en el uso para tipos de suelos distintos a los usados en el entrenamiento de las redes neuronales.
Estimación de la evaporación diaria en función de la profundidad de la napa
En la Figura 5 se presenta el flujo evaporativo EV (mm/día) estimado al resolver la ecuación (1) con el código VS2DTI, donde se observa profundidades en las cuales la evaporación calculada corresponde a a la evaporación potencial, denominada profundidad de desacoplamiento, entre 60 y 80 cm, con una fase de decaimiento abrupta y un flujo de evaporación despreciable bajo los 140 cm de profundidad. No se observan grandes diferencias en los valores de evaporación obtenidos.
La influencia de los parámetros p y n sobre el flujo de evaporación que se establece para diferentes profundidades del nivel freático se refleja directamente con la profundidad de desacoplamiento. Mientras mayor es el espesor de la franja capilar mayor es la profundidad hasta la cual el suelo es capaz de proporcionar el flujo demandado por la atmósfera. Esto se ve reflejado en las curvas de evaporación obtenidas mediante el modelo numérico en donde la curva de evaporación asociada a las características hidrodinámicas ajustadas a partir de los datos experimentales presenta una profundidad de desacoplamiento mayor a la obtenida para las curvas asociadas a los ajustes mediante los modelos H2 y H3. La franja capilar de la curva de succión ajustada a partir de los datos experimentales es del orden de 65 cm, mientras que el espesor de la franja capilar de las curvas de succión obtenidas a partir de la granulometría mediante los modelos H2 y H3 son del orden de 30 cm.
Análisis de sensibilidad de parámetros de VG en el flujo de evaporación
Como caso base de flujo de evaporación se utilizó el calculado con las curvas h(q) y K(q) ajustadas a las mediciones experimentales y se aplicó una variación del 5% sobre cada parámetro y se obtuvo el flujo evaporativo para diferentes profundidades de la napa: 80, 100 y 140 cm de profundidad. En la Tabla 2 se presenta la variación porcentual del flujo evaporativo obtenida según cada parámetro perturbado, en donde los valores negativos representan una disminución porcentual del flujo evaporativo con respecto a la situación base considerada.
Los parámetros qr y qs no afectan el comportamiento del modelo sobre el flujo evaporativo calculado cuando la napa se encuentra sobre los 100 cm de profundidad, aunque muestra una variación del orden de 0,2% cuando está a una profundidad de 140 cm. La sensibilidad asociada al parámetro HA no parece ser relevante. En todos los casos se obtuvo una variación menor al 0.1%, de modo que su influencia se considera despreciable ante variaciones de esa magnitud (±5%).
La conductividad hidráulica saturada presenta una sensibilidad similar para las distintas profundidades del nivel freático consideradas. Asimismo, se observa una variación en el flujo evaporativo calculado que es equivalente a la impuesta al parámetro Ks, ya que cuando se aumenta o disminuye el valor de Ks se produce un aumento o disminución de igual magnitud porcentual en el flujo evaporativo calculado., Esto se produce cuando la evaporación actual es menor que la evaporación potencial ya que el flujo evaporativo es proporcional a la conductividad hidráulica saturada y HA no varía ya que la temperatura es constante, tal como aparece en la ecuación 2.
El modelo se muestra muy sensible ante cambios en los coeficientes p y n. La mayor variación porcentual para estos parámetros se registró cuando la napa se encuentra a 140 cm de profundidad, con un aumento del flujo en un 193.11% para el parámetro de p y 119.20% para el caso de n con respecto a la situación inicial. Al igual que ocurre con el parámetro p, cuando n disminuye, el flujo evaporativo calculado aumenta y disminuye cuando n aumenta. Esto muestra que existe una relación negativa entre el flujo evaporativo y los parámetros p y n.
CONCLUSIONES Y RECOMENDACIONES
A pesar de las grandes diferencias observadas en la estimación de las curvas características h(q) y k(q) obtenidas a partir de ensayos experimentales y a partir de las PFT, estas no generan grandes diferencias en los valores del flujo de evaporación estimados en la modelación.
Las mayores diferencias observadas en las curvas características y en el flujo de evaporación estimado, la producen los parámetros p, n y Ks. Los modelos H2 y H3 sobreestiman los valores de p y Ks, y subestiman los valores de n.
El código VS2DTI es capaz de predecir el flujo evaporativo pero presenta algunas limitaciones. La condición de borde en la superficie produce una profundidad de desacoplamiento que parece estar determinada por el espesor de la franja capilar, suponiendo un flujo no saturado que no es real; el modelo resuelve las ecuaciones sólo para transporte de agua, para lo cual impone una condición de borde en función de la evaporación potencial generando una distancia de desacoplamiento que no existe en la realidad y que se debe a que el movimiento del agua líquida, vapor de agua y calor están fuertemente interrelacionados. Por lo tanto, las simulaciones y las soluciones analíticas que se basan sólo en la ecuación de Richards no son capaces de describir la totalidad del perfil.