Modelo de conductividad hidráulica dual para el movimiento del agua en suelos Macroporosos

Dual permeability model to water flow in Macroporous soils

Manuel Zavala¹, Heber Saucedo², Carlos Fuentes³, Carlos Bautista¹

¹ Universidad Autónoma de Zacatecas. Jardín Juárez 147, Centro Histórico, 98000, Zacatecas, Zacatecas, México, (mzavala73@yahoo.com.mx), (baucap@uaz.edu.mx).

³ Universidad Autónoma de Querétaro. Cerro de las Campanas, 76010, Santiago de Querétaro, Querétaro, México. (cfuentes@uaq.mx).

Recibido: mayo, 2011.
Aprobado: febrero, 2012.

Resumen

La alta presencia de macroporos en el suelo origina el desarrollo de flujos preferenciales que alteran significativamente las condiciones hidráulicas en el medio. El estudio detallado del flujo del agua para estas condiciones se aborda con la teoría de conductividad dual, que representa al suelo mediante dos sistemas porosos interconectados de propiedades hidráulicas contrastantes, uno representa los macroporos del suelo y el otro su matriz porosa. Sin embargo, este enfoque no considera las leyes de Laplace y Poiseuille para representar adecuadamente la influencia del tamaño de los poros de cada medio en la capacidad de retención de humedad y en la conductividad hidráulica. En este estudio se presenta un modelo numérico unidimensional para simular el flujo del agua a través de suelos con fuerte presencia de macroporos, que describe los procesos de transferencia desarrollados en ambos medios con dos ecuaciones de Richards acopladas. El modelo incorpora representaciones analíticas para las características hidrodinámicas del suelo, que consideran el efecto del tamaño característico de los poros de cada medio en la resistencia al flujo del agua y en la capacidad de retención. La discretización espacial de las ecuaciones diferenciales se realizó con el método de los elementos finitos tipo Galerkin y la integración en el tiempo con un método de diferencias finitas. La no linealidad de los sistemas resultantes de las discretizaciones se trata con un método iterativo de aproximaciones sucesivas. El modelo se aplica a la simulación de escenarios de infiltración con condiciones típicas del riego por aspersión y por gravedad para mostrar su capacidad de descripción.

Palabras clave: conductividad fractal, curva de retención, ley de Laplace, ley de Poiseuille, poro grande, poro geométrico.

The high presence of macropores in the soil causes the development of preferential flows significantly altering the hydraulic conditions in the medium. The detailed study of water flow for these conditions is discussed in the theory of dual conductivity, which represents the soil through two interconnected pore systems of contrasting hydraulic properties, one representing soil macropores and the other its porous matrix. However, this approach does not consider the laws by Laplace and Poiseuille to adequately represent the influence of the pore size of each medium in the water retention capacity and hydraulic conductivity. This study presents a one-dimensional numerical model to simulate water flow through soils with a strong presence of macropores, which describes the transfer processes developed in both media with two coupled Richards equations.The model incorporates analytical representations for the soil-hydrodynamic characteristics that consider the effect of the characteristic pore size of each medium on the resistance to flow water and water retention capacity. The spatial discretization of differential equations was performed using the Galerkin finite element method and time integration with a finite difference method. The nonlinearity of the systems resulting from the discretizations is treated with an iterative method of successive approaches. The model is applied to the simulation of infiltration scenarios with typical conditions of sprinkler and surface irrigation to show its power of description.

Key words: fractal conductivity, water retention curve, Laplace law, Poiseuille law, large pores, geometric mean pore.

INTRODUCCIÓN

Los suelos agrícolas pueden presentar fuerte heterogeneidad estructural por la presencia de grietas y fracturas cuya distribución espacial es altamente irregular. Estos macroporos, que separan la matriz del suelo (partículas sólidas), son altamente permeables y a través de ellos el agua puede moverse a velocidades considerablemente altas, generando condiciones de desequilibrio en la distribución de la presión del agua en la matriz y los macroporos. Tal condición limita la aplicación del enfoque tradicional de medio poroso simple al análisis de la dinámica del agua en el suelo (Gerke y van Genuchten, 1993; Šimunek y van Genuchten, 2008).

Una alternativa para considerar el efecto de la macroporosidad en el movimiento del agua en el suelo, conocida como aproximación multicontínua (Gwo et al., 1995), representa los macroporos y la matriz del suelo mediante dos continuos o regiones. En esta aproximación se usan diferentes formas de las ecuaciones de conservación de masa y movimiento para describir los procesos de transferencia de agua en cada región, acoplándose ambos medios a través del término de fuente de las ecuaciones de continuidad. Entre los modelos multicontinuos se destacan los modelos de porosidad dual (Zimmerman et al., 1996; Schwartz et al., 2000) y los modelos de conductividad dual (Gerke y van Genuchten, 1993; Vogel et al., 2000; Kodešová et al., 2010). Los primeros consideran que el flujo del agua sólo ocurre a través de los macroporos y que la matriz de suelo es exclusivamente una fuente de almacenamiento o abastecimiento, mientras que los modelos de conductividad dual toman en cuenta que los bloques de matriz de suelo forman un sistema interconectado en el que también ocurre transporte de agua.

Las componentes esenciales de los modelos de conductividad dual son los términos que gobiernan la transferencia de agua entre los sistemas que representan los macroporos y la matriz del suelo, y en la literatura se reportan expresiones empíricas (Othmer et al., 1991) y semi-empíricas (Gerke y van Genuchten, 1993; 1996) para representar estos términos de transferencia. En estos modelos se asumen las mismas relaciones funcionales de conductividad hidráulica y de retención de humedad para el sistema de macroporos y para el sistema de la matriz del suelo. Sin embargo, físicamente esta hipótesis no es correcta pues de acuerdo con las leyes de Laplace y Poiseuille el tamaño característico de los poros del medio define la capacidad de retención y la resistencia al flujo del agua, lo cual implica que los modelos de conductividad hidráulica y la característica de humedad no pueden ser idénticos para la matriz y los macroporos.

El objetivo de este estudio fue desarrollar un modelo numérico unidimensional para describir los procesos de transferencia de agua en suelos macroporosos, considerando la teoría de conductividad dual e incorporando funciones de conductividad hidráulica y retención de humedad adecuadas al tamaño característico de los poros de cada continuo o región.

En la teoría de medios porosos duales las propiedades del suelo se calculan como la suma de dos componentes locales, una asociada con el sistema de la matriz de suelo que se representará con el subíndice m, y la otra a los macroporos que se representará con el subíndice f.

Ecuaciones de base

De acuerdo con Gerke y van Genuchten (1993), las propiedades hidráulicas globales del suelo están relacionadas con las propiedades locales (matriz y macroporos), de la siguiente manera:

donde ø es la porosidad volumétrica [L³ L^-3]; θ es el contenido volumétrico de agua [L³ L^-3], que en un suelo parcialmente saturado es función del potencial de presión del agua en el suelo ψ; q es el flujo de Darcy [LT^-1]; w_f es la proporción volumétrica de macroporos respecto del volumen total de suelo; y w_m es el porcentaje de bloques de la matriz del suelo (agregados y microporos) w_m = 1 - w_f .

Considerando las relaciones (1.1)-(1.3) y la teoría de conductividad dual, el flujo del agua en la matriz y los macroporos se describe mediante dos ecuaciones de Richards (1931) acopladas:

donde ∇ es el operador gradiente [L^-1]; ψ es el potencial de presión del agua en el medio poroso [L]; K es la conductividad hidráulica [L^-1], que es función del potencial de presión y; z es la coordenada vertical orientada positivamente en dirección descendente [L]; t es el tiempo [T]; y Γ_w es el término que representa la transferencia de agua entre macroporos y matriz del suelo.

El término de transferencia Γ_w se considera proporcional a la diferencia de presiones entre el medio que representa los macroporos y el que representa la matriz del suelo:

donde a_w es un coeficiente de transferencia de agua [L^-1 T^-1].

El intercambio local de agua entre los medios causa un decremento en la presión en el dominio drenado y un incremento en la presión en el dominio receptor conforme a las curvas de retención de humedad respectivas. Gerke y van Genuchten (1993) analizan este proceso de transferencia en la interfaz matriz-macroporo considerando varias formas geométricas de los agregados y proponen la siguiente relación funcional para el coeficiente de transferencia:

donde a es una longitud característica que puede asumirse igual a la mitad del ancho de la estructura de la matriz de suelo [L], b es un coeficiente geométrico adimensional, g_w es un parámetro adimensional, y K_in (ψ) es la conductividad de la interfaz matriz-macroporo [LT^-1].

La conductividad K_in en la interfaz matriz-macroporo puede calcularse como el promedio aritmético de la conductividad hidráulica en cada medio:

donde K_sin es la conductividad hidráulica a saturación en la interfaz matriz-macroporo [LT^-1]; K_m y K_f son, respectivamente, la conductividad hidráulica en la matriz y macroporos [LT^-1]; K_sm y K_sf son la conductividad hidráulica a saturación en la matriz y macroporos [LT^-1].

Modelos funcionales para la conductividad hidráulica y retención de humedad

En los modelos de conductividad dual reportados en la literatura (Gerke y van Genuchten 1993; Vogel et al., 2000; Kodešová et al., 2010), se asumen las mismas relaciones funcionales de conductividad hidráulica y retención de humedad para los macroporos y la matriz del suelo. Sin embargo, de acuerdo con las leyes de Laplace y Poiseuille, la presión del agua en el suelo y la velocidad poral media dependen del tamaño de poro, por lo que no es adecuado usar el mismo modelo funcional para macroporos y microporos.

De acuerdo con Fuentes et al. (2001), el modelo del poro de la media geométrica es:

mientras que el modelo del poro grande es:

donde J es la variable de integración y Θ es la presión del agua como una función del grado de saturación efectivo definido por:

en la cual Θ_s y Θ_r son los contenidos volumétricos de agua a saturación y residual.

El parámetro s=D/E, donde D es la dimensión fractal del suelo y E=3 es la dimensión euclidiana del espacio físico, es una función de la porosidad volumétrica total del suelo (ø) definida de manera implícita por:

en las ecuaciones (7) y (8) considerando restricciones particulares entre m y n, permite derivar los modelos fractales para la conductividad hidráulica de la media geométrica y del poro grande:

con m y n parámetros de forma adimensionales, y ψ_d un parámetro de escala de la presión [L].

En el presente estudio se propone el uso de los modelos para la curva de retención de humedad y conductividad hidráulica (11) y (12) en el medio que representa la matriz del suelo; y usan las relaciones (11) y (13) para el medio que representa los macroporos.

Esquema numérico

Se desarrolla un modelo numérico para simular la transferencia de agua en columnas de suelos duales, teniendo como base las formas unidimensionales de las ecuaciones de Richards (2) y (3) así como las relaciones (11), (12) y (13).

La forma unidimensional de las ecuaciones (2) y (3) considerando (4), puede ser escrita:

Al sustituir las variables θ y ψ por las expresiones dadas en (16) y (17), se genera un error que se minimiza formando una integral ponderada sobre el dominio de solución y requiriendo que la integral sea cero. El uso de funciones de interpolación lineales permite obtener:

donde M es la matriz de masa, B es la matriz de rigidez, T es la matriz de transferencia, Q el vector de flujos en la frontera, y G es el vector de la fuerza de gravedad. Los coeficientes de estas matrices a nivel elemento son:

donde es el tamaño del elemento finito, G es la frontera del dominio, es la conductividad hidráulica promedio en el elemento, son funciones de interpolación del sistema de masa concentrado (Zienkiewicz et al., 2005).

La integración en el tiempo de las ecuaciones (18) y (19) se realiza con:

Utilizando un esquema de diferencias finitas implícito (w=1) se obtiene un sistema de ecuaciones no lineal, el cual se resuelve aplicando el método de Celia et al. (1990):

donde c(ψ)=dθ/dψ es la capacidad específica [L^-1], y p es el número de iteraciones en el intervalo de tiempo.

Las ecuaciones (24) y (25) al combinarse forman un sistema de ecuaciones con 2n incógnitas (acoplamiento numérico interno), donde la matriz de coeficientes es heptadiagonal simétrica. El sistema es resuelto utilizando el método de Cholesky para matrices simétricas almacenadas vectorialmente.

RESULTADOS Y DISCUSIÓN

La aplicación del modelo de conductividad dual requiere el conocimiento previo de los parámetros físicos e hidráulicos del medio poroso a analizar. En el presente estudio se utiliza y reinterpreta la información sintética (supuesta) de Gerke y van Genuchten (1993), la cual es considerada para ilustrar la versatilidad y capacidad de un modelo numérico de simulación (Gerke y van Genuchten, 1996; Vogel etal., 2000).

Los parámetros físicos e hidráulicos para la simulación de escenarios de infiltración en medios duales se presentan en el Cuadro 1. Los valores de θ_s, θ_r, K_s y w corresponden a los datos supuestos por Gerke y van Genuchten (1993); la dimensión cociente "s" se obtuvo resolviendo numéricamente la relación (10) con ø=θ_s; mientras que los parámetros m y ψ_d se obtuvieron ajustando las curvas de retención sintéticas de Gerke y van Genuchten (1993) con los modelos del poro de la media geométrica (11) y (12) y el modelo del poro grande (11) y (13), minimizando la raíz del error cuadrático medio (RECM). Los resultados de este proceso de optimización se presentan en las Figuras 1A y 1B, donde se muestra que los modelos del poro de la media geométrica y del poro grande tienen capacidad para describir la retención de humedad en suelos de textura fina y en suelos de textura gruesa considerando el tamaño predominante de los poros del suelo, que es lo físicamente correcto. En las Figuras 2A y 2B se comparan las curvas de conductividad hidráulica sintéticas de Gerke y van Genuchten (1993) con las curvas descritas con los modelos (12) y (13). En ambos casos las curvas descritas con los modelos de poro de la media geométrica y del poro grande indican que la conductividad hidráulica disminuirá mas rápidamente conforme el contenido volumétrico de agua decrece en el suelo, respecto de las tendencias descritas con el modelo de un solo poro característico de Gerke y van Genuchten (1993).

Gerke y van Genuchten (1993) indican que el parámetro g_w que interviene en la relación (5) es prácticamente independiente de la geometría de los agregados, propiedades hidráulicas y de las condiciones iniciales en el medio, por lo que puede usarse el valor constante g_w 0.4. Adicionalmente, para la simulación se asume que los agregados tienen forma rectangular; por tanto β=3 (Gerke y van Genuchten, 1996), y que la medida promedio de los bloques de matriz del suelo es de 2 cm, en consecuencia a=1cm.

Una vez conocidos sus parámetros, el modelo se aplica a la descripción de escenarios de infiltración en suelos macroporosos considerando condiciones asociadas al riego parcelario para ilustrar su utilidad en esta área.

 Escenario 1

Primero se analiza el escenario de infiltración planteado por Gerke y van Genuchten (1993), los quienes consideran una columna de suelo de 100 cm de longitud con una distribución inicial constante de la presión tanto en la matriz de suelo como en los macroporos ψ_m (z, 0)= ψ_f (z, 0)=-1000 cm. Se modela la infiltración originada por una intensidad de precipitación constante en la superficie de la columna, como la aplicada a través de aspersores, de q=50 cm/d; se asume que la infiltración se realiza sólo a través de la frontera superior de la columna (z=0)que representa los macroporos (q_f=q/w_f) y que la frontera superior de la matriz de suelo es impermeable (q_m=0); además las fronteras inferiores de ambos dominios (z=100) se suponen impermeables (q_m=q_f=0).

El dominio de solución se discretiza generando una malla uniforme de 1501 nodos y 1500 elementos (Δz=0.1 cm). Se usa un paso de tiempo inicial de Δt_i=1×10^-3 s que se incrementó en el transcurso de la simulación hasta un máximo de Δt_max=30 s en función del número de iteraciones realizadas para resolver el sistema acoplado (24) y (25). Se simula la aplicación de una intensidad de lluvia por aspersión durante 4 h.

En las Figuras 3, 4 y 5 se presentan los resultados de la simulación. El modelo de conductividad dual describe que el potencial de presión de agua en los macroporos aumenta rápidamente mientras que en la matriz de suelo el aumento de la presión es más lento (Figura 3). Como consecuencia se genera una diferencia de presión en la interfaz matriz-macroporo que origina la transferencia de agua de los macroporos a la matriz (Figura 4). En la Figura 5 se observa que la transferencia de agua entre ambos medios es mayor en los puntos próximos al frente de mojado y decrece gradualmente al alejarse del mismo. Además, en la Figura 4 muestra que el valor máximo de la transferencia de agua Γ_w disminuye en el tiempo debido a la atenuación de la diferencia de presión en el frente de mojado. La simulación de este escenario de infiltración permite mostrar la importancia del flujo de agua a través de los macroporos y su influencia en la distribución de humedad en el suelo; en el riego esto cobra mayor importancia en suelos finos agrietados.

Escenario 2

Se describe el flujo del agua en un medio dual originado por la aplicación de un riego por gravedad; se considera que la infiltración de agua por la superficie se da por la fractura y por la matriz del suelo. Se retienen los parámetros hidráulicos del suelo previamente determinados, variando las condiciones de frontera en la superficie de las columnas para representar el riego superficial (ψ_m (0, t)=ψ_f(0, t)=0), y su longitud z=150 cm.

Los resultados de la simulación se presentan en la Figuras 6, 7 y 8 donde se observa que los perfiles del potencial de presión, contenido volumétrico y término de transferencia de agua presentan irregularidades en su curvatura, sobre todo en los tiempos iniciales (1 y 2 h). En los tiempos cortos de la infiltración, la condición de Dirichlet homogénea origina que la presión en la matriz de suelo se equilibre rápidamente con la presión en el macroporo en las regiones próximas a la superficie del suelo, minimizando en esta zona el efecto de la gravedad y con un intercambio de agua gradual entre los medios. Sin embargo, conforme la profundidad de la columna aumenta el efecto de la condición de frontera se atenúa, acentuándose la diferencia entre el movimiento del agua a lo largo de la columna macroporosa y el movimiento en la columna de matriz del suelo. Esta diferencia en el movimiento origina que la presión a lo largo de la columna de macroporos aumente más rápido que en la matriz, generando diferenciales de presión que inducen transferencia de agua del macroporo a la matriz, teniendo tasas de transferencias de agua mayores que las registradas en zonas próximas a la superficie. Al transcurrir el tiempo, los perfiles de humedad y presión se estabilizan porque se reduce el efecto de la condición de frontera en la matriz del suelo por la baja conductividad hidráulica de este medio. Al igual que en el escenario anterior, el flujo de agua a través de los macroporos es trascendental para definir la distribución de humedad en el suelo.

CONCLUSIONES

Se desarrolló un modelo de conductividad dual que incorpora funciones de retención de humedad y conductividad hidráulica diferentes para la matriz del suelo y sus macroporos. Este modelo no incorpora más parámetros respecto a los de la propuesta tradicional de usar funciones idénticas para ambos medios. Este modelo se aplicó a la simulación escenarios de riego en suelos macroporosos obteniendo soluciones consistentes y libres de oscilaciones, lo cual muestra su confiabilidad numérica. El modelo es una herramienta útil para analizar alternativas de aplicación del agua en suelos macroporosos que minimicen la percolación, mejoren la uniformidad del riego y minimicen el movimiento de los fertilizantes más allá de la profundidad de las raíces de los cultivos.

La calibración experimental del modelo debe realizarse mediante procedimientos e información independiente para maximizar su capacidad predictiva. Las fracciones de macroporos y matriz deben estimar de muestras de suelo representativas del sitio, cuantificando los tamaños de los poros y aplicando la clasificación de macroporosidad. Para minimizar la cantidad de parámetros del modelo el contenido volumétrico residual puede suponerse nulo, aceptando la hipótesis de que la humedad en el suelo tiende a cero al aumentar la succión. Los parámetros m y n de la matriz y macroporos se estiman a partir de datos de porosidad volumétrica total y la curva granulométrica de cada medio. La caracterización hidrodinámica del suelo dual se completa determinando los parámetros de escala ψ_d y K_s de cada medio, iniciando de un evento transitorio de flujo de agua a través del suelo.

LITERATURA CITADA

Celia, M. A., E. T. Bouloutas, and R. L. Zarba. 1990. A general mass conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 7: 1483-1496.         [ Links ]

Fuentes, C., F. Brambila, M. Vauclin, J-Y Parlange, and R. Haverkamp. 2001. Modelación fractal de la conductividad hidráulica de los suelos no saturados. Ingeniería Hidráulica en México 16 (2): 119-137.         [ Links ]

Gerke, H. H., and M. Th. van Genuchten. 1993. A dual porosity model for simulating the preferential movement of water and solutes in porous media. Water Resour. Res. 29 (2): 305-319.         [ Links ]

Gerke, H. H., and M. Th. van Genuchten. 1996. Macroscopic representation of structural geometry for simulating water movement in dual porosity media. Adv. Water Resour. 19 (6): 343-357.         [ Links ]

Gwo, J. P., P. M. Jardine, G. V. Wilson, and G. T. Yeh. 1995. A multiple pore region concept to modeling mass transfer in subsurface media. J. Hydrol. 164: 217-234.         [ Links ]

Kodesová, R., J. Simunek, A. Nikodem, and V. Jirku. 2010. Estimation of the dual-permeability model parameters using tension disk infiltrometer and Guelph permeameter. Vadose Zone J. 9: 213-225.         [ Links ]

Othmer, H., B. Diekkrüger, and M. Kutilek. 1991. Bimodal porosity and unsaturated hydraulic conductivity. Soil Sci. 152 (3): 139-149.         [ Links ]

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

Schwartz., R. C., A. S. R. Juo, and K. J. McInnes. 2000. Estimating parameters for a dual-porosity model to describe non-equilibrium, reactive transport in a fine-textured soil. J. Hydrol. 229: 149-167.         [ Links ]

Šimunek, J., and M. Th. van Genuchten. 2008. Modeling nonequilibrium flow and transport processes using HYDRUS. Vadose Zone J. 7: 782-797.         [ Links ]

Van Genuchten, M. Th. 1980. A closed-form equation for predicting the hydraulic conductivity of the unsaturated soils. Soil Sci. Soc. Am. J. 44: 892-898.         [ Links ]

Vogel, T., H.H. Gerke, R. Zhang, and M.Th. van Genuchten. 2000. Modeling flow and transport in a two-dimensional dual-permeability system with spatially variable hydraulic properties. J. Hydrol. 238: 78-89.         [ Links ]

Zienkiewicz, O.C., R.L. Taylor, and J.Z. Zhu. 2005. The Finite Element Method. Its Basis and Fundamental. Ed. Elsevier, Amsterdam. 733 p.         [ Links ]

Zimmerman, R. W., T. Hadgu, and G. S. Bodvarson. 1996. A new lumped-parameter model flow in unsaturated dual porosity media. Adv. Water Resour. 19 (5): 317-327.         [ Links ]