Ecuaciones de base

En el programa de cómputo se resuelve la forma completa o simplificada de la ecuación de Boussinesq unidimensional del drenaje agrícola. En general, se puede escribir como:

donde H = H(x, t) es la carga hidráulica contada a partir de un estrato impermeable o nivel de referencia (L) y es función de la coordenada horizontal x y del tiempo t; T(H) es la transmisibilidad del acuífero (L ² T ^-1); para el caso de acuíferos libres es T(H) = K _s H; K _s es la conductividad hidráulica a saturación del suelo (LT ^-1); μ(H) es el coeficiente de almacenamiento del acuífero; R _w (t) es el volumen de recarga o descarga en la unidad de tiempo por unidad de área del acuífero (LT ^-1).

La recarga o descarga R _w (t) es una variable difícil de determinar, ya que depende de las condiciones de flujo existentes en la superficie de suelo (lluvia, evaporación, tirante de agua), el flujo del agua y el régimen de humedad en la zona vadosa, la estratigrafía de la zona no saturada del suelo, la macroporosidad, grietas, etcétera. En esta nueva versión del modelo de simulación numérico se incorporó el término R _w (t), representándolo mediante un polinomio que depende del tiempo:

donde a _r , b _r , c _r y d _r son coeficientes a calcular, considerando mediciones experimentales o estimaciones que se tengan de la evolución en el tiempo de la recarga para un evento específico.

Se consideró en el programa la relación entre el coeficiente de almacenamiento y la curva de retención de humedad establecida por ^{Fuentes et al. (2009)}:

donde θ _s es el contenido volumétrico de agua a saturación (L ³ L ^-3); y H _s es la elevación de la superficie del suelo (L). Es posible obtener una representación analítica explícita del coeficiente de almacenamiento a partir de la Ecuación (3) si la curva de retención de humedad del suelo es conocida. Se seleccionó la relación clásica de ^{Van Genuchten (1980)}, ampliamente aceptada en estudios de campo y de laboratorio por su flexibilidad descriptiva:

donde θ _r es el contenido volumétrico residual (L ³ L ^-3); ψ _d < 0 es un parámetro de escala del potencial de presión del agua en el suelo ψ, (L); m y n son parámetros de forma adimensionales. En esta nueva versión del programa de cómputo se incluyeron relaciones mecanicistas entre m y n que satisfacen la teoría de la infiltración, la restricción clásica de ^{Burdine (1953)}m = 1 - 2/n, así como las relaciones fractales de ^{Fuentes, Brambila, Vauclin, Parlange & Haverkamp (2001)}, que derivan en su estudio de la conductividad hidráulica de los suelos no saturados. Las restricciones fractales son las denominadas poro neutro m = (1 - 4s/n)/s, poro geométrico m = (1 - 2s/n)/s y poro grande m = (1 - 4s/n)/2s; donde s es la dimensión cociente del objeto fractal, tal como se considera el suelo en ese estudio, siendo s la razón entre la dimensión fractal del objeto D _f (^{Mandelbrot, 1982}; ^{Falconer, 2014}) y la dimensión del espacio de Euclides E, s = D _f /E.

Al introducir la Ecuación (4) en la Ecuación (3) se obtuvo la siguiente relación:

Se retuvo en esta segunda versión del programa la condición de radiación fractal para los drenes derivada por ^{Zavala et al. (2007)}:

donde q _d es el flujo de drenaje; γ es un coeficiente de conductancia adimensional; K _in = √(K _s K _d ) es la conductividad de la interfaz suelo-dren; K _s es la conductividad hidráulica a saturación del suelo; K _d la conductividad de la pared del dren; P, la profundidad de los drenes; H _d , la carga hidráulica en la posición del dren; D _o , el espesor del acuífero; s- = (1/2)(s ₁+s _d ) es la dimensión cociente en la interfaz suelo-dren; s ₁ es la dimensión cociente del suelo y s _d es la dimensión cociente de la pared del dren. La relación entre la porosidad areal (μ _areal ) y la porosidad volumétrica (φ) se obtiene de acuerdo con la idea probabilista como µ _areal = φ ^s φ^s= φ ^2s , se tiene que la porosidad volumétrica:

y la porosidad areal:

La conductividad de la pared del dren se obtiene con una fórmula basada en la ley de Poiseuille: K _d = (1/2)(g/ⱱ)µ _{areal_d} (R _HD )², donde g es la aceleración de la gravedad; ⱱ, la viscosidad cinemática del agua; R _HD , el radio hidráulico del dren. Se introdujo en el versión 2.0 de DRENAS la opción de manejar directamente la condición de radiación lineal (s- = 0.5 en la Ecuación (6)) sin necesidad de utilizar las rutinas de cálculo de s ₁ y s _d .

Modelaciones numéricas simplificadas del drenaje agrícola se pueden llevar a cabo imponiendo en los drenes condiciones de frontera tipo Dirichlet (carga hidráulica en la posición espacial del dren). En la nueva versión del programa de cómputo se incorporó esta condición de frontera como una alternativa a la condición de radiación (6); se seleccionó la siguiente función, que depende de la variable independiente tiempo:

donde x _dren es la coordenada horizontal donde se ubica el dren; y a _d , b _d , c _d y d _d son coeficientes que se deben calcular a partir de mediciones que se tengan de la evolución de la carga hidráulica sobre el dren. La Ecuación (9) tiene la flexibilidad de describir comportamientos extremos de la carga sobre el dren, como el abatimiento y ascenso.

Por último, para realizar la simulación numérica con la Ecuación (1) debe definirse el estado inicial de la carga hidráulica en el sistema. En este trabajo se incluyó la opción de usar hasta un polinomio de tercer grado, que es función de la coordenada horizontal x:

donde a _p , b _p , c _p y d _p son coeficientes a determinar a partir de mediciones de la carga hidráulica en el sistema en el instante inicial de cálculo o simulación.

En la solución numérica del sistema (1)-(10) se empleó el método del elemento finito tipo Galerkin para la discretización espacial; un esquema de diferencias finitas para la discretización temporal; el método iterativo de Picard para la linealización del sistema resultante, y un método de gradiente conjugado precondicionado de para la solución del sistema de ecuaciones algebraicas (^{Zavala et al., 2014}; ^{Zienkiewicz,‎ Taylor, & Zhu, 2013}; ^{Noor & Peters, 1987}). El sistema de ecuaciones resultante de la discretización se programó en Visual Basic 2017.

Es conocido que cálculos simplificados de separación entre drenes y abatimiento del manto freático pueden ser realizados aplicando soluciones analíticas obtenidas para formas reducidas de la ecuación de Boussinesq (Ecuación (1)); por ejemplo, para régimen de flujo permanente se tiene la fórmula de ^{Hooghoudt (1940)}, y en régimen transitorio las relaciones de Glover-Dumm (^{Dumm, 1954}) y ^{Fuentes et al. (1997)}, las cuales ya se incluían en la versión original del programa de cómputo y se retienen en la actualización.

Desarrollo de la interfaz gráfica

En la segunda versión del programa de cómputo DRENAS, se reprogramaron en Visual Basic 2017 todos los módulos de captura de información de la versión original; se desarrollaron bases de datos internas independientes de Microsoft Access; se desvinculó de Microsoft Excel el tratamiento de las gráficas, y se generó el archivo ejecutable para instalar DRENAS 2.0 en cualquier computadora que disponga de sistema operativo Windows de 64 bits, e incluso se generó también la opción para sistemas Windows 32 bits.

El programa DRENAS 2.0 tiene dos módulos de cálculo: uno denominado soluciones analíticas y otro modelo numérico (Figura 1). El módulo soluciones analíticas tiene cuatro secciones: dos para realizar cálculos para condiciones de flujo de agua en régimen permanente usando la solución de ^{Hooghoudt (1940)} y dos para condiciones de flujo en régimen transitorio; uno para realizar cálculos con la solución de ^{Dumm (1954)} y el otro con la solución de ^{Fuentes et al. (1997)}. En estas cuatro secciones se puede calcular la separación entre drenes o el módulo de drenaje. Respecto a la versión original del modelo DRENAS, en la versión 2.0 estas secciones fueron mejoradas, al incorporar las siguientes opciones: a) almacenamiento de las simulaciones realizadas en una base de datos interna; b) carga de ejemplos previos; c) generación de un reporte de la simulación y alternativa para exportarlo a un archivo pdf; d) La gráfica de resultados también puede ser exportada como imagen. Un ejemplo de una sección de cálculo del módulo de soluciones analíticas se presenta en la Figura 2.

Figura 1 Pantalla general del programa DRENAS 2.0. 

Figura 2 Una sección de cálculo del módulo soluciones analíticas. 

El módulo de la solución numérica se programó en una forma Tab Index (en inglés), que contiene seis pestañas o secciones: las primeras cuatro están desarrolladas para la captura de datos referente al sistema de drenaje, características del dren, propiedades del suelo y datos necesarios para la simulación numérica, como son los relacionados con el control de tiempo, discretización espacial, condiciones límite y visualización gráfica de resultados; la quinta pestaña sirve para almacenar los datos del ejemplo en una base de datos interna, y la última pestaña contiene el botón para ejecutar la simulación numérica.

Las nuevas capacidades del módulo de simulación numérica de la versión 2.0 de DRENAS que se programaron se establecieron de la siguiente forma:

a) Primero de la base de datos de suelos UNSODA 2.0 (^{Nemes, Schaap, Leij, & Wösten, 2001}) se seleccionaron las curvas experimentales de retención de humedad de 208 suelos reportados en la literatura; se ajustaron con la Ecuación (4), considerando las restricciones mecanicistas entre m y n antes descritas (^{Zavala, Saucedo, & Fuentes, 2018}), y los resultados obtenidos se incluyeron en una base de datos interna del programa de cómputo para ponerlos a la disposición de los usuarios del mismo (Figura 3). Si se selecciona la opción de usar datos de esta base del programa, al elegir un suelo, las secciones de propiedades físicas, hidráulicas y coeficiente de almacenamiento se llenan automáticamente con sus parámetros correspondientes. La mejora a resaltar en esta nueva versión es que se pueden manejar las relaciones entre m y n de los modelos de Burdine, poro neutro, poro geométrico y poro grande (Figura 4).

Figura 3 Módulo numérico: base UNSODA ampliada. 

Figura 4 Módulo numérico: propiedades suelo-agua. 

b) En la solución numérica de la ecuación de Boussinesq se incorporó el término de recarga R _w (t) como una función del tiempo (Ecuación (2)); se realizaron las modificaciones de programación necesarias para poder imponer en los drenes la condición de frontera tipo Dirichlet (Ecuación (9)), y se programó también la condición inicial tipo polinomio cúbico (Ecuación (10)). En la Figura 5 se presenta la nueva sección desarrollada para incluir estas tres condiciones.

Figura 5 Módulo numérico: condiciones límite y recarga. 

c) Se desarrolló una nueva sección para registrar los datos de la simulación, cargar datos para comparar simulaciones y para seleccionar el tipo de gráfica a visualizar en tiempo real que no depende de programas externos (Figura 6).

Figura 6 Módulo numérico: gráficas de la simulación. 

Validación 1

Se seleccionó la solución analítica del tipo Hooghoudt presentada en ^{Fragoza et al. (2003)}, que considera radiación lineal (s- = 0.5 en la Ecuación (6)) en los drenes para describir la evolución de la superficie libre en un sistema de drenaje subterráneo con recarga constante (R _o ) y flujo de agua en régimen permanente:

donde H _c es la carga hidráulica al centro entre drenes. Si se define h(x,t) = H(x,t) - D _o y considera la condición de radiación lineal se tiene:

donde h _d es la carga hidráulica sobre el dren (x=0 y x=L); D _o es el espesor del acuífero; y γ el coeficiente de conductancia de la condición de radiación lineal que se usa en los drenes. Analizando la Ecuación (2) se tiene para este caso que a _r =b _r =c _r =0 y d _r =R _o .

Los datos del sistema de drenaje son los reportados por ^{Fragoza et al. (2003)}, derivados de un experimento de campo realizado en el distrito de riego 076 Valle del Carrizo, Sinaloa. El sistema de drenaje evaluado tiene una separación entre drenes L = 50 m, profundidad media de los drenes P = 1.5 m, profundidad del estrato impermeable D _o = 3.5 m, elevaciones del estrato impermeable y de la superficie del suelo H _i =0.0 m y H _s =5.0 m. La parcela tiene suelo de textura arcillosa, con una conductividad hidráulica a saturación de K _s = 0.557 m/d, y el coeficiente de almacenamiento µ = 0.1087. Se asume γ= 1.5 y h _c = 0.5 m. Con la Ecuación (12) se tiene h _d = 0.365 m y el valor de recarga obtenido con la condición para régimen permanente donde la descarga por el dren es igual a la recarga, R _o = 0.000944 m/d. Con estos valores se aplicó la Ecuación (11) para obtener la distribución espacial de la carga hidráulica a lo largo de la separación entre drenes.

Se reprodujo con el programa de cómputo el escenario descrito, modelándose el proceso transitorio flujo de agua en el sistema hasta alcanzar el régimen estacionario, el cual corresponde al descrito con la solución analítica. Se asumió la condición inicial de carga hidráulica constante H(x, t=0) = 4.50 m (por tanto, a _p = b _p = c _p = 0 y d _p = 4.50 m en la Ecuación (10)) y para tener las mismas condiciones de la solución analítica, en el programa se usó la condición de frontera de radiación lineal, s-= 0.5 en la Ecuación (6). Se consideró la hipótesis de igualdad entre la conductividad del suelo y la de la pared del dren K _s = K _d , lo cual genera K _in = K _s = K _d . Se reinterpretó el valor del coeficiente de conductancia reportado por ^{Fragoza et al. (2003)}, para adaptarlo a la forma de radiación que se maneja en DRENAS 2.0, obteniéndose γ= 0.045. La discretización espacial se realizó con el número de nodos N = 1000 (Δx = 0.05 m), y en la discretización del tiempo se usó Δt _ini = 1.157E-05 d, Δt _max = 6.94E-04 d, y Δt _min = 1.157E-06 d. El tiempo de simulación fue de 720 días.

En la Figura 7 (a y b) se presentan los resultados de la evolución del manto freático a la semana y al año de drenaje libre. Se observa que la solución numérica tiende a la solución analítica conforme avanza el tiempo y la reproduce de forma correcta. Esta comparación es una evidencia de que los módulos de captura de información de la solución numérica del programa y en particular el módulo de cálculo que maneja el término de la recarga vertical, condición inicial y condición de frontera, operan bien y están libres de errores de programación; también se observa que las series simuladas son estables y libres de oscilaciones numéricas.

Figura 7. Simulación de un escenario de drenaje agrícola (coeficiente de almacenamiento constante y recarga vertical constante). 

Validación 2

^{Zavala, Fuentes & Saucedo (2004)} realizaron un experimento en un módulo de drenaje en el que se instalaron dos drenes de 30 cm de longitud (l _D ) y 5 cm de diámetro (D _D ), con N _o = 233 perforaciones circulares de diámetro d _o = 0.158 cm, distribuidas de modo uniforme en la superficie de cada tubo. Las dimensiones de cada dren: L = 100 cm, P = 120 cm y D _o = 25 cm. El módulo se rellenó con suelo de textura arenosa, que se saturó aplicando una carga de agua constante en la superficie del suelo hasta eliminar el aire. Con los drenes tapados se removió el exceso de agua en la superficie y se cubrió para eliminar la evaporación. Por último, se destaparon los drenes y se midió durante 10 días (240 horas) el volumen de agua drenado; se debe notar que la condición inicial corresponde a H(x,0) = P + D _o (d _p = 145 cm en la Ecuación (10) y el resto de sus coeficientes son nulos) y la recarga es nula (R _w = 0 en la Ecuación (2)).

La porosidad del suelo fue Φ₁ = 0.5396 cm³/cm³, y con este valor se resolvió la Ecuación (7), obteniéndose s ₁ = 0.7026. El valor de la conductividad hidráulica a saturación reportado es K _s = 18.3 cm/h. El coeficiente de almacenamiento del acuífero descrito con la Ecuación (5) para los parámetros θ _s = Φ₁, θ _r = 0, ψ _d = -41.8 cm y n = 2/(1-m) = 3.19 (restricción de ^{Burdine, 1953}).

La porosidad areal del dren es µ _{areal_d} = A _o /A _d = 0.0097 cm²/cm² y al resolver la Relación (8) se tiene s _d = 0.5688. Con el radio hidráulico de los orificios circulares de la pared del dren R _HD = 0.0395 cm, la viscosidad cinemática del agua ⱱ = 36 cm²/h se aplica la ley de Poiseuille para calcular la conductividad de la pared del dren K_d = 2,670.8 cm/h. Los valores de la conductividad en la interfaz suelo-dren y la dimensión relativa son K_in = 221.08 cm/h y s-= 0.6357. El coeficiente γ de la radiación fractal reportado por ^{Zavala et al. (2004)} se reinterpretó para adaptarlo a la forma de radiación que se maneja en DRENAS 2.0 (Ecuación (6)), obteniéndose γ= 0.0749.

El coeficiente de almacenamiento variable y la condición inicial de saturación total del medio poroso originan que la simulación deba ser ejecutada empleando pasos de tiempo pequeños al inicio de la misma, lo cual minimiza problemas de convergencia y estabilidad numérica. Los pasos de tiempo utilizados son Δt _ini = 2.77E-04 h, Δt _ini = 2.77E-05 h y Δt _max = 0.2 h. El número de nodos para discretizar el espacio es de 200.

En la Figura 8 (a y b) se presenta, respectivamente, el proceso de captura de datos de la simulación y la evolución de la lámina drenada. En la Figura 8b se puede observar que el programa de cómputo reproduce de manera satisfactoria los datos experimentales de lámina drenada, teniéndose una solución numérica estable y libre de oscilaciones, lo cual es un indicador de que el módulo de coeficiente de almacenamiento variable trabaja de forma correcta. Los pasos de tiempo y espacio utilizados en la simulación son adecuados, ya que la lámina drenada total calculada (23.915 cm) difiere de la lámina medida (23.92 cm) en sólo 0.03%.

Figura 8. Simulación de un escenario de drenaje agrícola (coeficiente de almacenamiento variable, recarga nula y condición de radiación fractal en los drenes).