El paquete para flujo en zona no saturada de ModFlow es uno de sus complementos más recientes y todavía se encuentra en fase experimental pues necesita de modelos probados en más partes del mundo para comprobar su validez (que hasta ahora es innegable).

El flujo vertical en zona no saturada es gobernado por diversos fenómenos difíciles de cuantificar, por lo que su cálculo es también muy complejo. Uno de los mejores enfoques usados actualmente para modelar el flujo en zona no saturada es la ecuación de Richards.

Sin embargo, ésta es altamente no lineal debido a la mutua dependencia entre el contenido volumétrico de agua, la conductividad hidráulica y la succión, además de ser muy complicada de resolver matemáticamente debido a que requiere la solución de ecuaciones en las tres dimensiones (cosa que para esta situación normalmente no es necesario porque el flujo es predominantemente vertical). Por todo lo anterior se opta por una simplificación del procedimiento. Así, la forma unidimensional de la ecuación de Richards se aproxima por una ecuación de onda cinemática que se resuelve por el método de las características.

Despreciando el término de difusividad hidráulica y combinando lo restante con la ecuación de variación de Abbott, se llega a un sistema de ecuaciones en el que se debe imponer que las tasas de variación del contenido volumétrico de agua respecto al tiempo y al espacio sean indeterminadas a través de las líneas características.

Así, se obtienen tres ecuaciones características que definen la velocidad de onda:

, donde v(θ) es la velocidad característica restringida al eje Z.

La onda representa un frente húmedo en la zona no saturada que se origina por el incremento repentino de la infiltración del suelo.

Un aumento de la tasa de cambio volumétrico del agua produce un frente de humectación que se representa con una onda de conducción y un decremento, un frente de secado, que se representa mediante una onda de arrastre. Si una onda de arrastre alcanza a una de humectación, la atenúa, haciendo que adopte su contenido de agua y reduciendo su velocidad de propagación. Al contrario, cuando una onda de humectación alcanza a otra más lenta, o a una onda de arrastre, la remueve, y prevalece su contenido de agua, resultado en una rehumectación. Cada vez que una onda alcanza a otra se produce un shock y esto ocasiona discontinuidad en propiedades hidráulicas como contenido de agua y flujo.

En contraste con un frente de humectación que se mantiene fuerte, las ondas de arrastre se dilatan con el tiempo debido a la gravedad. En consecuencia, las ondas de arrastre deben ser divididas en una serie de ondas incrementales o ser representadas por una función que describa el drenaje interno con el tiempo.

La relación de velocidades de dos puntos en una onda de arrastre es proporcional a su relación de profundidades:

Ahora, debido a que el efecto de la ET crea una pendiente no linear en el contenido de agua sobre el perfil de onda de conducción, y ésta es discontinua al despreciar el término de la difusividad hidráulica, no se pueden utilizar ecuaciones de conductividad hidráulica -K(θ)- como la función Brooks-Corey, por lo que se debe recurrir a la integración sobre un volumen de control que contiene un solo frente de humectación equivalente de igual masa para poder resolver la ecuación 1 que considera la difusión hidráulica.

Así, se llega a demostrar la siguiente relación:

Por otro lado, la pendiente es continua a lo largo de una onda de arrastre, por lo que es totalmente aplicable la función Brooks-Corey:

Donde el contenido residual de agua se puede aproximar a la retención específica y ésta a su vez se puede hallar restando la porosidad menos el rendimiento específico. El exponente Brooks-Corey es un valor adimensional que incluye los efectos de la distribución del tamaño de poros y, por lo tanto, de la succión capilar.

La velocidad del punto más profundo del frente de secado se puede hallar con:

Una relación entre el punto más profundo y los demás puntos de la onda de arrastre es determinada por la siguiente ecuación:

Sin embargo, las ecuaciones ya no son aplicables cuando se interceptan las ondas. Para estos casos, un enfoque muy simple es dividir la onda de arrastre en incrementos y calcular sus velocidades en cada incremento por diferencias finitas.

Se decide con el parámetro NSTRAIL la cantidad de incrementos, pero con más de 10 se predice un error en el balance de masa de menos del 0.05 %. Además, se debe tener en cuenta que debido al recálculo ineludible de los cambios en las propiedades de las ondas ocasionados por los shocks para mantener la continuidad, la discretización temporal del paquete UZF es específica del problema resuelto, por lo que los escalones de tiempo usados son diferentes a los usados en el esquema principal de ModFlow.

Cabe resaltar que se desprecian los efectos de gradiente hidráulico negativo producidos por la succión ejercida por las plantas dentro de su profundidad de extinción de evapotranspiración (ET) debido a que pueden originar fuerzas difusivas importantes por la absorción por las raíces; asegurando que las ecuaciones de onda cinemática puedan seguir siendo válidas. Sin embargo, sí se sustrae agua de acuerdo a los valores de tasa de demanda de ET mediante la pérdida instantánea de agua a lo largo de un intervalo de profundidades igual a la profundidad de extinción de ET. Si la demanda no es satisfecha, el agua puede ser extraída desde el acuífero siempre que la napa freática esté por encima de dicha profundidad de extinción.

ET se especifica como una tasa (longitud por unidad de tiempo) dentro del paquete UZF1 (variable de entrada PET) y se convierte internamente a una tasa por unidad de profundidad dividiendo la tasa de ET especificada por la profundidad de extinción (EXTDP).

Las últimas dos ecuaciones dadas sirven para hallar el contenido de humedad en la cabeza del frente húmedo y la profundidad de los puntos a lo largo del intervalo de análisis.

Límite de infiltración

Se pueden obtener las siguientes ecuaciones para modelar la infiltración:

El contenido de agua se establece en el contenido de agua saturada cuando la tasa de infiltración que se especifica en el input excede la conductividad hidráulica vertical saturada. Si este es el caso, entonces la diferencia se multiplica por el área de vista en planta de la celda del modelo y esta tasa volumétrica de agua se puede añadir a un tramo de río específico o a un lago, estableciendo la variable IRUNFLG mayor que cero.

Un algoritmo fue escrito en FORTRAN utilizando todas las ecuaciones anteriores para simular el flujo vertical en zona no saturada y ET. Este algoritmo mantiene un registro de la ubicación, contenido de agua, la velocidad y el flujo de todas las ondas en la zona no saturada a través del tiempo y la interacción de las ondas de adelantamiento entre sí. Además, los cambios en el espesor de la zona no saturada, debido a los cambios en la elevación de la napa freática, se consideran con respecto a las ondas que llegan a la napa freática y contribuyen al almacenamiento del agua subterránea.

Acoplamiento a MODFLOW

El flujo de la zona no saturada como la recarga a la napa freática en un acuífero no confinado se resta del lado derecho del sistema de ecuaciones que se resuelven por MODFLOW-2005:

UZF1 no provoca inestabilidades adicionales o lenta convergencia de MODFLOW-2005, cuando el agua subterránea no está descargando a la superficie. Los niveles de agua subterránea que se elevan por encima de la superficie terrestre añaden una dependencia de carga adicional al utilizar el paquete UZF1 debido al cálculo de la descarga de agua subterránea a la superficie y, por lo tanto, la convergencia puede ser significativamente más lenta cuando el agua subterránea está descargando a la superficie desde muchas celdas del modelo.

Una tasa de infiltración en la superficie se especifica para cada columna vertical en el dominio del modelo activo. La tasa de infiltración se multiplica por el área de vista en planta de la celda para producir una tasa de flujo volumétrico. Se utilizan tres opciones en el UZF1 para determinar donde el flujo en la zona no saturada recarga al agua subterránea (variable NUZTOP):

La primera opción especifica que toda la recarga se concentra sólo en la capa superior del modelo,
La segunda opción permite especificar la recarga en una capa determinada en cada columna vertical y
La tercera opción permite la recarga en la celda activa más superficial en una columna vertical.

Cuando la última opción se utiliza, el programa continuará encauzando el agua a través de la zona no saturada hasta alcanzar la napa freática de la celda activa más superficial. Sin embargo, las propiedades hidráulicas de la zona no saturada permanecerán constantes, incluso si la zona no saturada se extiende por más de una celda en la dirección vertical. Si no hay celdas activas en una columna debajo de la superficie, entonces, la infiltración y la recarga no se producen y un estado de alerta aparecerá como output. Si una celda se activa durante una simulación y el elemento correspondiente a la matriz UZAF, IUZFBND, es distinto de cero, entonces, la infiltración y la recarga empezarán a producirse.

Si la napa freática se eleva en la zona no saturada, el acuífero producirá una cantidad de agua igual a la elevación de la napa freática multiplicada por la diferencia entre el contenido de agua y la retención específica. Alternativamente, si la napa freática desciende debido a la descarga del agua subterránea o bombeo a algún otro lugar, entonces el contenido de agua en la zona no saturada entre la napa freática antigua y nueva se fija igual a la de la retención específica.

La recarga que ocurre cuando se eleva la napa freática en un escalón de tiempo se calcula por ondas de enrutamiento a través de la zona no saturada, usando como base de la zona no saturada la elevación de la napa freática al fin del escalón de tiempo anterior. El flujo que cruza la base en el escalón de tiempo se acumula y el volumen de agua en la zona no saturada a través del cual la napa freática se elevó durante el escalón de tiempo se añade al volumen que fluyó pasando la base. Esta suma se divide por el escalón de tiempo para obtener una tasa de recarga volumétrica. Una vez que la ecuación de flujo ha sido resuelta para el escalón de tiempo, una nueva base de la zona no saturada en cada celda activa se calcula. La recarga que ocurre durante una disminución de la napa freática se calcula basándose en el flujo de agua desde la zona no saturada durante el escalón de tiempo usando como base de la zona no saturada a la predicha en la iteración o escalón de tiempo anterior.

Descarga a la superficie

La tasa volumétrica de descarga de agua subterránea a la superficie terrestre se calcula sobre la base de la siguiente ecuación:

El paquete UZF1 incluye una opción para especificar una matriz bidimensional de valores (variable IRUNBND) que identifica donde la descarga de agua subterránea a la superficie se agregará a un arroyo o lago si los paquetes SFR2 y LAK están activos. Además, la variable IRUNDFLG debe ser mayor que 1 para que el agua se añada a los arroyos o lagos. La descarga de agua subterránea a la superficie de la tierra y la infiltración en exceso de la conductividad hidráulica vertical saturada se añade instantáneamente a los ríos y lagos especificados.

Un presupuesto de agua para la zona no saturada es seguido independiente del presupuesto de agua subterránea en MODFLOW-2005. El balance de masa de la zona no saturada representa toda el agua entre la superficie terrestre y el manto freático. El balance de masa para la zona saturada se tiene en cuenta en el presupuesto de agua subterránea. Cuando la napa freática se eleva por encima de la parte superior de una celda del modelo, entonces el balance de masa de la zona no saturada deja de ser calculada para esa celda, y el almacenamiento no saturado y el cambio en la zona de almacenamiento no saturado para esa celda se ponen en cero.

No hay limitaciones de procedimiento relativas a períodos de estrés y las longitudes de escalones de tiempo asociadas con el paquete UZF1. Sin embargo, debido a que UZF1 retrasa la recarga de la napa freática, los cambios de estrés no corresponden necesariamente con el comienzo de un período de estrés. Por lo tanto, se debe tener precaución cuando se simula escalones de tiempo que aumentan durante un período de estrés (especificando la variable TSMULT > 1,0 en la entrada de datos para el archivo de discretización) ya que el agua que se filtra a través de la zona no saturada puede alcanzar el nivel freático en el final de un período de estrés, cuando la longitud del escalón de tiempo es máxima.

Simulaciones en estado estacionario

El estado estacionario se puede especificar como el período de estrés inicial de una simulación cuando se utiliza el paquete UZF1. Cuando se da este caso, la tasa de recarga de agua subterránea calculada en UZF1 se fija igual a la tasa de infiltración y la ET no se remueve de la zona no saturada.

El contenido de agua es constante dentro de cada celda del modelo durante el estado estacionario. La tasa de ET especificada se resta de la tasa de infiltración para calcular una tasa de recarga de estado estacionario. La ET también se puede producir a partir de agua subterránea cuando el nivel freático está por encima de la profundidad de extinción de ET. Un perfil de contenido de agua uniforme se calcula para cada celda activa, que se utiliza de forma automática como las condiciones iniciales para un siguiente período de estrés transitorio.

No se requieren valores del contenido inicial de agua cuando se especifica el período de estrés inicial como el estado estacionario. Sin embargo, una opción para el desarrollo de los perfiles iniciales de contenido de agua es ejecutar el modelo para varios años mientras se repite un ciclo de tasas de infiltración. Esta etapa se puede llevar a cabo después de una simulación en estado estacionario.

Limitaciones

Gradientes de presión negativas también pueden dar lugar a la redistribución lateral y vertical en la zona no saturada. Debido a que se ignoran los gradientes de presión negativa, UZF1 probablemente subestima el avance del frente húmedo en los primeros tiempos.

No hay ninguna consideración para el período anterior a la formación de charcos cuando la infiltración debe estar representada por un flujo constante igual a la intensidad de lluvia, más que el contenido de agua constante asignado al límite. Esto implica cierto error ya que diversos autores muestran la dependencia entre la infiltración y la infiltración acumulada antes del tiempo de encharcamiento. Así, para la mayor precipitación o las tasas de deshielo, el modelo subestima la infiltración durante el inicio de ésta. UZF1 no simula una franja capilar encima del nivel freático y supone que el agua subterránea es liberada de forma instantánea desde el almacenamiento.

El paquete UZF calcula la recarga y la evapotranspiración, por lo que el RCH, EVT y los paquetes de ETS no deben ser utilizados normalmente junto con el paquete UZF. Sin embargo, MODFLOW no impide que UZF sea utilizado en conjunto con el RCH, EVT y los paquetes de ETS.

El paquete UZF no tiene en cuenta los cambios en la superficie terrestre debido a subsidencia.

Método de cálculo de la profundidad y el ancho, los caudales de entrada, las propiedades de cauces, las dimensiones y los coeficientes de rugosidad de Manning para los segmentos de los arroyos

Corriente o arroyo: Las propiedades de cauces y dimensiones flujo para cada tramo del arroyo se especifican para el primer y último alcance en cada segmento. Una interpolación lineal se utiliza para determinar cada valor en el punto medio de cada tramo.

El método determina cómo se calcula la profundidad y la anchura de corriente para cada tramo de un segmento.

La profundidad y la anchura del flujo se especifican cuando ICALC es 0,
La profundidad de la corriente se calcula a partir de la ecuación de Manning asumiendo un amplio canal rectangular cuando ICALC es 1;
La profundidad y la anchura del flujo se calculan a partir de la ecuación de Manning con una sección transversal de ocho puntos cuando ICALC es 2;
La profundidad y la anchura del flujo se calculan a partir de una función de potencia cuando ICALC es 3, y
La profundidad y la anchura se calculan a partir de una tabla de valores cuando ICALC es 4.

Caudal de entrada: Flujo de entrada calculado a partir de 1 o más segmentos tributarios. El coeficiente de rugosidad de Manning no se utiliza para el cálculo de la profundidad o la anchura de la corriente cuando ICALC es 0, 3, o 4, y el coeficiente de rugosidad de Manning para el desbordamiento del río no se utiliza cuando ICALC es 1.

Ejemplo:

Conclusiones hechas sobre modelos realizados

Los patrones de recarga son el resultado de los efectos acumulados de ET, variación espacial de recarga, conductividad hidráulica, topografía y espesor de la zona no saturada.
La recarga de agua subterránea sumada sobre toda la cuenca indica dos patrones distintos: una respuesta casi inmediata a los picos de la infiltración en las tierras bajas del valle, donde la zona no saturada es delgada, y una respuesta difusa que varía lentamente durante varios meses en las zonas de montaña, donde la zona no saturada es relativamente gruesa.
La tasa de ET volumétrica acumulada de la cuenca está limitada por la disponibilidad de agua en la zona no saturada y la profundidad a la napa freática, excepto durante los períodos de alta infiltración. La cantidad de pérdida de ET es sensible a la profundidad de extinción y contenido de agua a dicha profundidad, de manera que la atención se debe poner en la estimación de la profundidad de extinción y contenido de agua de extinción para reducir la incertidumbre en las estimaciones de recarga.
Un año de muy por encima del promedio de recarga puede tener un efecto prolongado en la descarga de agua subterránea a la superficie de la tierra y a los arroyos.
Los perfiles de contenido de agua por debajo de las celdas individuales del modelo indican que el contenido de agua sigue una envolvente estacional que se hace más estrecha con la profundidad debido al patrón estacional de la infiltración y ET.
La marcada diferencia entre la infiltración, recarga y caudal sugiere que la zona no saturada podría tener un fuerte efecto de las fluctuaciones en los niveles de agua subterránea y caudal durante y después de las tormentas. Los años en los que las tasas de infiltración superan la tasa anual media pueden aumentar el almacenamiento de agua subterránea, incluso en rocas relativamente poco permeables y dar lugar al aumento del flujo de corrientes que pueden extenderse mucho más allá del período de infiltración.

Comentarios hechos sobre modelos realizados:

Una constante Brooks-Corey exponente de 4 y un contenido de agua saturada constante de 0,3 fueron asignados a todas las células.
El almacenamiento específico se establece en 5 × 10-7 m-1 y el rendimiento específico se especificó como 0,05 en las crestas y 0,25 en los valles cerca de los riachuelos.
Las corrientes estuvieron representados por 15 segmentos formados por 201 tramos usando el paquete SFR2 (Streamflow-Routing).
La conductividad hidráulica vertical de la zona no saturada tenía los mismos valores que la conductividad hidráulica horizontal y vertical utilizado en el paquete LPF (Layer-Property Flow).
El período de simulación se divide en 365 periodos de estrés de un día. El primer período de estrés es el estado estacionario, y los períodos de estrés posteriores fueron transitorios. Un escalón de tiempo se simuló para cada período de estrés.
Condiciones de no-flujo fueron simuladas a través de la parte inferior y los lados del modelo a excepción de tres celdas por debajo y adyacentes al río en la salida de la cuenca, que se especifica como células de carga constante.
Las tasas de infiltración variaron en el modelo y fueron tres veces mayores en las cimas de cordilleras, en comparación con las tierras bajas de valle.
El rango de las tasas de infiltración aplicados se calcularon sobre la base de la diferencia media entre los dos medidores de precipitación ubicados en las partes superior e inferior de la cuenca.
El paquete PCG (Preconditioned Conjugate-Gradient) se utilizó para resolver la carga de agua subterránea en todo el dominio del modelo.
El criterio de la carga de cierre utilizado en la simulación transitoria era 0,0009 m y el criterio de flujo de cierre fue de 0.005 m3/d. El error de balance de masa absoluta de la zona no saturada varió de 0 a 0,08 por ciento para todos los períodos de estrés y la zona saturada fluctuó entre 0 y 0,22 por ciento. El error de balance de masas acumulada al final de la simulación de 365 días era 0 por ciento de la zona no saturada y 0,10 por ciento para la zona saturada.
La mayor parte del caudal simulado durante los primeros seis meses fue de descarga de agua subterránea a la corriente. El flujo superficial total, incluyendo la descarga de agua subterránea a la superficie, más el exceso de saturación de la tierra, superó la descarga de agua subterránea a las corrientes a partir de mediados de mayo hasta mediados de agosto. Después de mediados de agosto, el flujo total de tierra se redujo y la descarga de agua subterránea a los arroyos dominó el hidrograma.
El aumento de la superficie sobre la que el agua subterránea descargó en la superficie, provocó el aumento de la escorrentía superficial (descarga de agua subterránea a la superficie, más el exceso de saturación de la tierra).
El frente húmedo creado por muchos tipos diferentes de infiltración tardó más de un año para alcanzar el nivel freático en las zonas donde la conductividad hidráulica vertical saturada era baja y la profundidad de agua subterránea era mayor de 20 m por debajo de la superficie terrestre.

Referencias

Niswonger, R.G., Prudic, D.E. & Regan, R.S. Documentation of the Unsaturated-Zone Flow (UZF1) Package for Modeling Unsaturated Flow Between the Land Surface and the Water Table with MODFLOW-2005. United States Geological Survey. A Product of the Ground-Water Resources Program. Chapter 19 of Section A, Ground Water, of Book 6, Modeling Techniques. Techniques and Methods 6-A19. Available at: http://pubs.usgs.gov/tm/2006/tm6a19/pdf/tm6a19.pdf