La ecuación de Green y Ampt

La ecuación de Green y Ampt se establece a partir de la ecuación de continuidad y la ley de ^{Darcy (1856)} con las siguientes hipótesis: a) el perfil de humedad inicial en una columna de suelo es uniforme: θ = θ_o, b) la presión del agua en la superficie del suelo es hidrostática: ψ = h ≥ 0 (donde h es la lámina de agua sobre la superficie del suelo), c) existe un frente de humedecimiento bien definido caracterizado por una presión negativa: ψ = -h_f () 0 (donde h_f es la succión (1) en el frente de humedecimiento), y d) la región entre la superficie del suelo y el frente de humedecimiento está completamente saturada (flujo en pistón): θ = θ_s y K = K_s (donde K_s es la conductividad hidráulica a saturación; es decir, el valor de la conductividad hidráulica de la ley de Darcy correspondiente al contenido volumétrico de agua a saturación). La ecuación diferencial ordinaria resultante es la siguiente:

donde I es la lámina infiltrada acumulada, t es el tiempo, y θ_s y θ_o son los contenidos de humedad a saturación e inicial, respectivamente.

Si el tirante de agua sobre la superficie se considera independiente del tiempo, la Ecuación (1) se integra analíticamente con la condición inicial I = 0 en t = 0, lo cual resulta en lo siguiente:

Cuando se desprecia la cantidad de agua que fluye en el infinito, la lámina infiltrada como lo suponen Green y Ampt es igual al volumen por unidad de área almacenado en el pistón: I(t) = z_f(t)Δθ con Δθ = θ_s - θ_o, y z_f(t) es el frente del pistón (^{Fuentes et al., 2012}).

Optimización de parámetros

En los modelos de regresión no lineal, cada observación y_i se escribe en términos de una función de respuesta f(x_i; β). Una diferencia importante de los modelos de regresión no lineal es que el número de parámetros β de la regresión no se relaciona directamente con el número de variables x_i en el modelo (^{Cornejo-Zuniga & Rebolledo-Vega, 2016}). Para estimar los parámetros K_s y h_f a partir de una prueba de infiltración, se utiliza el algoritmo Levenberg-Marquardt (^{Moré, 1978}), el cual ha sido una técnica estándar para problemas de mínimos cuadrados no lineales (^{Chávez et al., 2022}; ^{Fuentes, et al., 2022}; ^{Šimůnek & Hopmans, 2018}). La actualización de parámetros se realiza de manera iterativa con la siguiente expresión:

donde J es la matriz Jacobiana relacionada con las variaciones de la función de infiltración respecto a cada parámetro a optimizar, I_d es la matriz identidad, r es el vector de diferencias entre la lámina infiltrada medida y la calculada con el algoritmo de optimización, y μ es el parámetro de amortiguación que se ajusta en cada iteración. La matriz Jacobiana se calcula de la siguiente manera:

donde I_m es la ecuación de Green y Ampt, y m es el número de datos medidos. En este algoritmo los valores de la matriz Jacobiana se aproximan numéricamente mediante derivadas centradas para mejorar el tiempo de cómputo.

La optimización de los parámetros se codificó en un programa con funciones y subrutinas, las cuales comienzan con un par de valores iniciales, para posteriormente calcular la matriz Jacobiana y el vector de errores. Con esto es posible resolver la Ecuación (3), que da como resultado un par de valores nuevos. Lo anterior permite calcular el vector de errores y revisar si es menor a la iteración anterior dada una tolerancia. En caso de ser mayor a la tolerancia, el algoritmo inicia de nuevo con el par encontrado en la última iteración, y si este vector de errores es menor a la tolerancia, el programa termina e imprime los resultados. El diagrama completo se muestra en la Figura 1.

Figura 1. Diagrama de flujo de Levenberg-Marquardt para la optimización de parámetros.  

Parámetros de entrada

Los datos de entrada que se requieren son los contenidos de humedad inicial y a saturación, así como los datos obtenidos de una prueba de infiltración (tiempo y lámina infiltrada). Para obtener el contenido de humedad inicial (θ_o), las variables de acceso son la densidad aparente del suelo (ρ_t), el contenido gravimétrico de humedad (ω_o) y la densidad del agua (ρ_w):

La porosidad (φ) es el volumen del espacio vacío de los medios porosos, la cual se calcula a partir de la densidad aparente (ρ_t) y la densidad de las partículas de cuarzo (ρ_s):

donde ρ_s = 2.65 g∙cm^-3.

El contenido de humedad a saturación (θ_s) es el volumen de agua contenido en el espacio poroso, normalmente asimilada a la porosidad volumétrica (φ) por la siguiente desigualdad: 0 ≤ θ_s ≤ φ; sin embargo, es importante denotar que, en un suelo aparentemente saturado con agua, generalmente queda atrapada una cierta cantidad de aire. Por lo tanto, el contenido volumétrico de humedad a saturación puede ser tomado como θ_s = 0.9φ (^{Haverkamp et al., 2016}; ^{Rogowski, 1971}). En este trabajo, θ_s será tomado como la porosidad total: θ_s ≈ φ.

Experimento de laboratorio

Se tomaron muestras de suelo de parcelas agrícolas de la región de San Juan del Río, Querétaro, México, siguiendo la metodología propuesta por ^{Reynolds y Topp (2007)}. El suelo se tamizó con malla número 10 (2 mm) y se secó al aire libre durante una semana. Antes de colocar el suelo en la columna, se tomó una muestra para enviarla al laboratorio con el fin de obtener el contenido de humedad inicial por el método gravimétrico.

Para realizar la prueba de infiltración, se utilizaron dos columnas de acrílico de diferente longitud y sección transversal circular. Ambas columnas tenían una placa porosa cubierta por un filtro y colocada en la parte inferior con el propósito de retener el suelo y permitir la salida del agua y el aire (Figura 2). El interior de las columnas se recubrió parcialmente con cera para crear rugosidad entre el suelo y la columna de acrílico. El suelo se colocó en la columna en capas de 5 cm de espesor a una densidad similar a la obtenida en campo. Durante la prueba de infiltración se mantuvo un tirante de agua constante.

Figura 2. Representación esquemática de una columna de suelo.  

En el Cuadro 1, adaptado por ^{Saucedo et al. (2013)}, se presentan los valores medios de algunos parámetros del suelo (θ_o, θ_s, K_s y h_f) a partir de la textura del suelo, partiendo de las relaciones porpuestas por ^{Rawls et al. (1991)}. Para conocer estos valores, los suelos utilizados se analizaron en el laboratorio con el hidrómetro de Bouyoucos utilizando los métodos de la norma mexicana NOM-021-SEMARNAT-2000 (^{Secretaría del Medio Ambiente y Recursos Naturales [SEMARNAT], 2002}) y el triángulo de texturas del Departamento de Agricultura de Estados Unidos de América (USDA). Dentro del funcionamiento del algoritmo de optimización se agregó la siguiente condición: si la diferencia de parámetros en la iteración actual con la anterior es menor a 0.00001 %, o el número de iteraciones es igual a 20, termine la ejecución del algoritmo y muestre los resultados.

Codificación del algoritmo

El algoritmo de optimización se programó en el lenguaje Visual Basic.net en el entorno de desarrollo integrado Microsoft Visual Studio 2019. Este programa se puede instalar en cualquier plataforma de Windows 10, y la velocidad en la memoria RAM que se requiere es mínima, por lo que cualquier ordenador con un procesador de gama media y 4 GB de RAM lo soporta. El programa consta de cuatro secciones principales: a) menú principal, b) datos desplegados de una prueba de infiltración, c) resultados gráficos de la lámina de riego medida vs la optimizada y d) resultados de la optimización (Figura 3).

Figura 3. Ventana principal del programa de cómputo desarrollado.  

En el Cuadro 2 se muestran las rutinas utilizadas para el adecuado funcionamiento del programa de cómputo desarrollado. Además, se describe la acción que realizan cada una de ellas.

Pruebas de infiltración

Las condiciones iniciales del experimento se describen en el Cuadro 3, donde se muestran los resultados de los análisis obtenidos en el laboratorio, así como las longitudes de las columnas (L) rellenas de suelo, los diámetros (D) y los tirantes de agua en la superficie (h).

Con los datos de las condiciones iniciales de los experimentos y la lámina infiltrada en el tiempo, se optimizaron los parámetros K_s y h_f de la ecuación de Green y Ampt en función de la textura del suelo. Los resultados de la optimización se muestran para cada sitio en la Figura 4, donde se aprecia el ajuste del modelo de Green y Ampt a los valores medidos.

Figura 4. Comparación entre la lámina infiltrada medida y la calculada con la ecuación de Green y Ampt.  

El Cuadro 4 muestra la optimización obtenida, así como la raíz del error cuadrático medio (RECM), que es resultado de la lámina infiltrada medida y la lámina infiltrada teórica. Se observan errores menores a 0.45 cm, lo cual indica que la optimización de parámetros con valores iniciales relacionados con la textura del suelo funciona de manera adecuada. Es importante resaltar que, a partir de la gráfica anterior, se puede observar el buen ajuste que presenta el modelo optimizado con los datos medidos en laboratorio.

Validación numérica

La validación del algoritmo numérico se realizó para eliminar errores de programación y revisar la consistencia de soluciones correctas. Se usaron datos reportados en la literatura, cuyos resultados se obtuvieron mediante la solución unidimensional de la ecuación de Richards (^{Chávez et al., 2016}; ^{Fuentes et al., 2020}). Los datos utilizados corresponden a un suelo franco-arenoso de Montecillo, México (^{Zataráin et al., 2003}). Los valores de los parámetros de este suelo son: ρ = 1.3607 g∙cm^-3, θ_o = 0.1391 cm³∙cm^-3, h = 1.50 cm, φ = 0.4865 cm³∙cm^-3, θ_s = φ y L = 70 cm. Mediante el programa de cómputo se optimizaron los parámetros K_s y h_f de la ecuación de Green y Ampt, en función de la textura del suelo (Figura 5).

Figura 5. Comparación de la lámina infiltrada obtenida con la ecuación de Richards y la modelada con la ecuación de
Green y Ampt.  

Al comparar el resultado obtenido y el reportado en la literatura (Cuadro 5), se observa que no existe una diferencia considerable, ya que se considera igual la raíz del RECM (0.1953 cm), y ésta nos indica que existe poca variación del resultado medido contra el modelo de Green y Ampt optimizando K_s y h_f .

El parámetro h_f es característico de la ecuación de Green y Ampt; por ello, no puede ser comparado con ningún otro parámetro de otras ecuaciones. Sin embargo, el valor se ajusta dependiendo de las proporciones de arena, limo y arcilla del suelo, por lo que se puede considerar inversamente proporcional a K_s, ya que cuando éste disminuye, el valor de h_f aumenta.

Si únicamente se considera la textura del suelo para realizar el diseño del riego por gravedad, el valor medio de K_s para el suelo franco-arenoso analizado sería de 2.9 cm·h^-1 (Cuadro 1). Sin embargo, en este tipo de suelo predominan las características de la arena, lo cual conlleva a que se tengan valores de K_s mayores a 1.5 cm·h^-1. El rango de variación entre ambos valores es muy amplio, lo cual genera un cambio sustancial en los valores característicos de la infiltración.

Por lo anterior, se recomienda realizar pruebas de infiltración en columnas de suelo para representar unidimensionalmente el flujo de manera adecuada, ya que, en otro tipo de pruebas como el método de doble cilindro, normalmente se cuenta con un flujo bidimensional y, además, ocurre un flujo preferencial por las paredes de los cilindros que en la mayoría de las veces se le atribuyen a la velocidad de infiltración. La ventaja al utilizar columnas con muestras de suelo alteradas en laboratorio es que permite tener control de todas las variables, como la densidad, la profundidad, la lámina de agua, la humedad inicial, entre otras.

Sensibilidad temporal de los parámetros optimizados de la prueba de infiltración

Se realizó un análisis de sensibilidad de los parámetros de Green y Ampt para verificar cómo se propagan los errores de estimación de los dos parámetros a las estimaciones de la infiltración. Se hicieron optimizaciones con los resultados experimentales en diferentes tiempos para cada sitio. En S1 se utilizó un incremento acumulado de Δt = 1 h, en el S2 un Δt = 4 h y en S3 un Δt = 1 h. En la Figura 6 se observa que los parámetros presentan la misma tendencia, mientras el K_s disminuye gradualmente, el valor de h_f aumenta. El proceso inverso ocurre con la misma tendencia; esto indica que el algoritmo optimiza los valores de acuerdo con la textura del suelo y entre más datos de infiltración se tienen, los parámetros optimizados pueden variar considerablemente, ya sea disminuyendo o aumentando sus valores de forma gradual e inversamente proporcional entre ellos. Asimismo, se muestra la importancia de realizar una prueba de infiltración en tiempos largos, ya que, si se realiza en tiempos cortos, los poros no se encuentran completamente llenos de agua y solo se tendría un valor de conductividad hidráulica y no la conductividad hidráulica a saturación.

Figura 6. Evolución de los parámetros K_s y h_f.  

El programa está delimitado para que cada par de nuevos valores que se propongan estén dentro de los límites de los parámetros; es decir K_s )( 0 y 0 () h_f ≤ 200. Los parámetros se ajustan en función de la textura del suelo, haciendo énfasis en la Ecuación (7):

donde S es la sorbilidad del medio poroso (^{Philip, 1957}). Sin el uso de parámetros iniciales relacionados con la textura del suelo, pueden existir una infinidad de parámetros de ajuste de los datos medidos en una prueba de infiltración.

En este trabajado se utilizó un ordenador con las siguientes características: Procesador Intel® CoreTM i7-4710MQ CPU @ 2.50 GHz y memoria RAM de 32 GB. El tiempo de cómputo requerido para encontrar los valores óptimos de K_s y h_f en los experimentos fue de 2 s. Esto debido a que la solución de Green y Ampt no requiere del conocimiento de valores en el nivel de tiempo anterior; por lo tanto, es posible resolverla exclusivamente en el tiempo requerido, como en los tiempos exactos donde se compara con los datos medidos. Por su parte, la ecuación de Richards requiere los valores de presión en el nivel de tiempo anterior, lo cual retrasa, en al menos una iteración, el cálculo de la lámina infiltrada.