## Articulo

• Similares en SciELO

## versión On-line ISSN 2007-2422

### Tecnol. cienc. agua vol.9 no.1 Jiutepec ene./feb. 2018  Epub 24-Nov-2020

#### https://doi.org/10.24850/j-tyca-2018-01-06

Artículos

Modelación bidimensional de la infiltración del agua en surcos aplicando el gradiente conjugado

1Universidad Autónoma de Querétaro, Querétaro, México, chagcalos@uaq.mx

2Universidad Autónoma de Querétaro, Querétaro, México, juan_mota_0@live.com.mx

3Instituto Mexicano de Tecnología del Agua, Jiutepec, México, cfuentes@tlaloc.imta.mx

4Universidad Autónoma de Querétaro, Querétaro, México, toniokv2@gmail.com

Resumen

Palabras clave Richards; bidimensional; diferencias finitas; riego por surcos

Abstract

Knowledge of the bulb moisture distribution during an irrigation event helps to understand the depth of wet soil better. The flow of water in soil is modeled with Richards’ equation, however, most solutions require a lot of computational time or the solution is bounded to a specific domain. In this study, the two-dimensional Richards’ equation is solved using two schemes: one implicit in the entire depth and a mixed scheme: explicit in the unsaturated zone and implicit in the saturated zone. The implicit scheme was solved with the conjugate gradient for non-linear problems using a sinusoidal function on the surface. The comparison was performed between the two schemes and the results of wetting profiles and infiltrated depth were similar, given an error criterion. The infiltration process was simulated in a furrow irrigation event through a constant depth in the three furrows, and in a second event they were imposed a different depth to each one. The moisture distribution profiles obtained during the simulated time can help us to give recommendations for the disadvantages of establishing a different depth in the furrows during an irrigation event.

Keywords Richards; two dimensional; finite differences; furrow irrigation

Introducción

Las dos modalidades más comunes para el riego por gravedad son riego por melgas y riego por surcos. Estos procesos ocurren, como es natural, en el espacio 3D. Cuando la geometría del espacio y las condiciones de frontera son invariantes respecto a algún eje coordenado, el proceso de infiltración puede simplificarse y considerarse como un proceso unidimensional; este es el caso del riego por melgas, en el cual la geometría del espacio es parecida a un paralelepípedo. En el riego por surcos, la geometría del proceso es más compleja y requiere que se consideren al menos dos variables espaciales para su simulación correcta (López-Avendaño, Palacios-Vélez, Fuentes-Ruíz, Rendón-Pimentel, & García-Villanueva, 1997).

La ecuación que representa el movimiento del agua en un suelo parcialmente saturado es obtenida mediante la combinación del principio de continuidad y la ley de Darcy. Si se introduce la capacidad específica como la derivada de la curva de retención CΨ=θ , la ecuación resultante queda expresada en términos de una variable θ(Ψ), denominada ecuación de Richards (1931), que en su forma bidimensional adquiere la forma:

CΨΨt=xKΨΨx+zKΨΨz-KΨz (1)

donde θ es el contenido de humedad; Ψ, la presión; x, una variable espacial (horizontal); z, la profundidad; K, la conductividad hidráulica y t es el tiempo.

Para resolver esta ecuación se incorporan las características hidrodinámicas del suelo: la curva de retención de humedad y la curva de conductividad hidráulica. Las funciones más utilizadas para representarlas son la ecuación de van Genuchten (1980), y la ecuación de Brooks y Corey (1964), respectivamente:

θΨ-θrθs-θr=1+ΨΨdn-m (2)

Kθ=Ksθ-θrθs-θrη (3)

donde θ s es el contenido de humedad a saturación, generalmente asimilada a la porosidad total del suelo (ϕ); θ r , el contenido de humedad residual; Ψ, la presión; Ψ d , un valor característico de la presión; K s , la conductividad hidráulica saturada; η, m y n son parámetros de forma adimensionales y positivos.

Los parámetros m y n son relacionados a partir de la restricción propuesta por van Genuchten (1980) al modelo de la conductividad hidráulica de Burdine (1953): m = 1-2/n. El parámetro η es relacionado con m y n, y la porosidad total mediante (Fuentes, Zavala, & Saucedo, 2009):

η=2s2mn+1 (4)

donde s = D f /E, en la cual D f es la dimensión fractal del medio poroso y E = 3 es la dimensión de Euclides del espacio físico, ligado a la porosidad total con la ecuación implícita siguiente:

1-ϕs+ϕ2s=1(5)

Existen algunos trabajos en la literatura que estudian la ecuación de Richards en dos dimensiones: Zavala-Trejo, Fuentes-Ruíz y Saucedo-Rojas (2005) estudiaron el drenaje agrícola subterráneo, para lo cual utilizaron la ecuación de Richards bidimensional (2D). Saucedo, Zavala, Fuentes y Castanedo (2013) acoplaron la ecuación de Barré de Saint Venant unidimensional con la ecuación de Richards 2D para modelar el riego por melgas. López-Avendaño et al. (1997) analizaron la infiltración de agua en un dominio tipo surco utilizando la ecuación de Richards 2D. Sin embargo, la solución que aplican para resolverlo es con la discretización de la derivada de manera convencional: agregando un factor de peso o parámetros de interpolación para el tiempo y el espacio. En este estudio se resuelve la ecuación de Richards 2D ―aplicado al riego por surcos― usando dos esquemas: un esquema implícito para todo el perfil y un esquema mixto: explícito en la zona no saturada e implícito en la saturada. El esquema implícito se resolvió con el gradiente conjugado para problemas no lineales usando en la superficie una función sinusoidal. El modelo propuesto está adaptado para simular la variación del tirante en el tiempo (condición de Dirichlet) y en la zona seca del surco se puede imponer una condición de tipo Neumann, por lo que pueden aplicarse condiciones de evaporación o una función de transpiración de acuerdo con las necesidades que se requieran.

Materiales y métodos

Creación del dominio espacial

El dominio espacial para resolver la ecuación de Richards se creó con las siguientes etapas: 1) selección y discretización de las fronteras del dominio; 2) relleno de las regiones interiores con puntos aleatorios distribuidos uniformemente, y 3) triangulación con el algoritmo de Delaunay (1934). Los dos tipos de fronteras existentes son exteriores (donde los puntos de interés están en el interior) e interiores (donde los puntos de interés están en el exterior).

Para modelar la infiltración en los surcos se supuso una frontera cuadrada y el lado superior se sustituyó por una función periódica, en este caso, se optó por una función sinusoidal, f(x) = 15 sin (πx/45), que da un surco de 30 cm de altura y 90 cm de lomo a lomo (figura 1). Una vez que se obtuvo la frontera, el siguiente paso fue obtener la discretización de la misma.

Para rellenar el interior primero se obtuvo el cuadrado más pequeño que contiene a todos los puntos frontera. Una vez obtenido este cuadrado C: [a, b] x [c, d] se buscó el punto P k usando una distribución uniforme dentro del cuadrado C. Para saber si P k estaba dentro de la frontera, se encontró la suma de los ángulos desde P k de cada par de puntos contiguos de la frontera (figuras 2 y 3).

Para un punto en el interior de una frontera, la suma de los ángulos mencionados es de ±360 grados y para un punto en el exterior de la frontera la suma es de cero grados. Al descartar los puntos exteriores queda un dominio (figura 3(a)). La gran cantidad de puntos en la figura 3(a) que se generaron en este primer paso es porque se requiere que las regiones interiores estén lo más densas posibles. Posteriormente, se procedió a eliminar puntos que estuvieran a una distancia que se consideró conveniente, por ejemplo, todos aquellos que estuvieron a una distancia menor de 5 mm uno del otro (figura 3(b)), para finalmente proceder a la unión de los puntos restantes, obteniéndose como resultado una triangulación del dominio, como el mostrado en la figura 3(c). El dominio con más puntos tendrá una mayor cantidad de nodos, lo que implica un costo computacional mayor para encontrar la solución.

No obstante que la función que se utilizó para representar un surco es periódica y que con una mitad sería suficiente para modelar la infiltración, en este trabajo se hace toda la ejemplificación usando tres surcos; esto debido a que se realiza una simulación con cargas diferentes en cada uno de los surcos con la finalidad de observar el comportamiento del perfil de humedecimiento.

Solución en diferencias finitas para la ecuación de Richards 2D

La ecuación (1) evaluada en algún punto t en el tiempo se escribe de la siguiente forma:

CΨkΨkt=xKΨkΨkx+zKΨkΨkz-KΨkz(6)

En esta ecuación, Ψ(k) es una función de tres variables Ψ (x, z, t). La ecuación (1) con derivada temporal discreta adquiere las siguientes dos formas (además de otras, dependiendo del esquema y precisión elegidos):

CΨkΨk+1-Ψktk+1-tk=xKΨkΨkx+zKΨkΨkz-KΨkz(7)

CΨkΨk-Ψk-1tk-tk-1=xKΨkΨkx+zKΨkΨkz-KΨkz(8)

La ecuación (7) corresponde al modelo explícito y la ecuación (8) corresponde al modelo implícito. Así, se define el lado derecho de la ecuación (8) como el operador Op (Ψ k ); las ecuaciones (7) y (8) toman una forma más simple:

CΨkΨk+1-Ψktk+1-tk=Op(Ψk)(9)

CΨkΨk-Ψk-1tk-tk-1=Op(Ψk)(10)

Cuando C [Ψ (k)] ≠ 0, la solución que se obtiene es una solución recursiva:

Ψk+1=tk+1-tkCΨkOpΨk+Ψk(11)

Sin embargo, en la zona saturada C k ) = 0, por lo que en esos puntos no puede aplicarse la ecuación (11) . Para estos casos hay que optimizar la ecuación (10) . El algoritmo que se utilizó en este trabajo es el del gradiente conjugado para problemas no lineales, también llamado “gradiente conjugado no lineal”. El método requiere que se defina una función a minimizar, en nuestro caso requerimos que se cumpla la ecuación (10) , por lo que naturalmente una función de error es la siguiente:

EΨk=CΨkΨk-Ψk-1tk-tk-1-OpΨk(12)

Siguiendo el procedimiento descrito por Dai y Yuan (1999), a partir de la ecuación (12) , definida como la función de error, y con el operador definido a partir de la ecuación (8) , el procedimiento consta de los siguientes pasos: 1) se escogió x0=Ψ(k) ; 2) x0=-xE(x0) ; 3) α0=arg minEx0+αx0 ; 4) x1=x0+α0x0 y 5) S0=x0. Del paso 1 al paso 5 es la parte donde inicia x0, los siguientes pasos (6 - 10) son iterativos y el algoritmo termina cuando Ex0 es pequeño: 6) xn=-xE(xn) ; 7) β=ΔxnTΔxn/Δxn-1TΔxn-1 ; 8) Sn=Δxn+βSn-1 ; 9) αn=argminExn+αΔxn y 10) xn+1=xn+αnSn. Por último se escoge xn+1 como la solución a las presiones en el tiempo t. Aunque es posible aplicar el algoritmo para todo el vector de presiones, por razones de rendimiento los cambios se realizaron sólo en la zona saturada, debido a que, en la zona no saturada, es posible aplicar la regla de recurrencia. La condición inicial para x0 es Ψ(k) en la zona no saturada con el método explícito y Ψ(k-1) en la zona saturada.

Método para resolver el operador OpΨk

Para simplificar la notación se hicieron las siguientes simplificaciones: se define Ψi como la presión en el nodo i en cualquier tiempo; la coordenada espacial del nodo i se define como Ci=xi,zi , donde xi es la coordenada en la longitud y zi en la profundidad; Ci,j es la coordenada espacial del j-ésimo nodo relacionado con el nodo i, es decir Ci,j = Ci,j=xi,j,zi,j ; mj es la cantidad de nodos relacionados con el nodo i, y Pi es la terna xi,zi,Ψi=Ci,Ψi. Utilizando derivación de orden uno, la derivada en Ci es estimada derivando el plano ax+bz+c, que mejor ajusta a los puntos Pi,P1,i,P2,i,,Pmi,i. Una forma de aproximar dicho plano es utilizar la matriz pseudoinversa definida por Cj,i y resolviendo para a, b y c; los puntos Pi,P1,i,P2,i,,Pmi,i satisfacen ax+bz+c, por lo que axi+bzi+c=Ψi y axi,j+bzi,j+c=Ψi,j. Las ecuaciones anteriores definen el siguiente sistema de ecuaciones:

Aabc=xizi1xi,1zi,11xi,mizi,mi1abc=ΨiΨi,1Ψi,mi(13)

La matriz A de la ecuación (13) no necesariamente es invertible; sin embargo, partiendo de la hipótesis que todo sistema se puede optimizar (Penrose, 1955), queda definida como:

abc=A+ΨiΨi,1Ψi,mi(14)

donde A+ es la matriz pseudoinversa de A.

Con los valores a, b y c obtenidos, es posible derivar hacia cualquier dirección, en particular, hacia x y z. De manera similar a la aproximación lineal, utilizando derivación de orden dos, la derivada en Ci se estimó derivando la superficie ax2+bxy+cy2+dx+ey+f que mejor ajustara a los puntos Pi,P1,i,P2,i,,Pmi,i. Una forma de aproximar dicha superficie es utilizar la matriz pseudoinversa definida por Cj,i y resolviendo para a, b, c, d, e, y f; así se obtuvieron los puntos Pi,P1,i,P2,i,,Pmi,i que satisfacen ax2+bxy+cy2+dx+ey+f, por lo que también satisfacen axi2+bxiyi+cyi2+dxi+eyi+f=Ψi y axi,j2+bxi,jyi,j+cyi,j2+dxi,j+eyi,j+f=Ψi,j. Estas ecuaciones quedan definidas de la siguiente manera:

Aabc=xizi1xi,1zi,11xi,mizi,mi1abc=ΨiΨi,1Ψi,mi(15)

donde aplicando la misma consideración que en la solución de la ecuación (14), la matriz A queda definida como:

abcdef=A+ΨiΨi,1Ψi,mi(16)

donde A+ es la matriz pseudoinversa de A.

Al final de este procedimiento, con los valores a, b, c, d, e, y f obtenidos, se calculó la derivada.

El volumen infiltrado por unidad de longitud de cauce en la unidad de tiempo (A˙I)se calculó con la siguiente expresión (Fuentes, De León-Mojarro, & Hernández-Saucedo, 2012):

A˙I=AIt,AIx,t=PmqIx,y,z,tdPm(17)

donde qIx,y,z,t es el caudal de infiltración por unidad de superficie de canal o caudal unitario y Pm es el perímetro mojado en cada surco.

Aplicación

Comparación entre el esquema implícito y el esquema mixto

Los dos esquemas (implícito y mixto) se compararon. Los valores que se utilizaron fueron del suelo Montecillo (Fuentes, 1992; Saucedo, Zavala, & Fuentes, 2011; Saucedo et al., 2013), los cuales son θr=0.0000cm3cm-3, θs=0.4865cm3cm-3, Ks=1.8400cmh-1, Ψd=-45.9000cm, m=0.1258 y una profundidad L=90cm. En esta simulación se utilizó una distancia promedio entre nodos en ambos esquemas de Δz=1 cm, y una discretización temporal de Δt=0.01 h para el método mixto, mientras que para el implícito se usó un Δt=0.00128 h. Los resultados para 240 min (figura 4) indican que no hay diferencias entre ambos esquemas, pero el esquema mixto requiere una discretización temporal mucho más fina que el esquema implícito. Sin embargo, el esquema implícito requiere más cálculos para cada tiempo, lo que hace la simulación de un evento de infiltración un poco más lenta.

El tiempo de cómputo requerido para simular los 240 min, usando una computadora de 6 GB en RAM, para el método mixto fue de 18 min, mientras que para el método implícito fue de 24 min aproximadamente, es decir, el método implícito requirió de un 30% de tiempo adicional que el método mixto, sin embargo, con los recursos computacionales de hoy en día, el tiempo de simulación puede acortarse un poco más.

Es importante mencionar que este tiempo de cómputo corresponde a la discretización de los tres surcos; si esta simulación se realizara únicamente en la mitad de un surco, los tiempos para obtener la solución se reducirían con la misma discretización espacial y temporal a 4.5 y 6 minutos, respectivamente. El resultado de las curvas de infiltración de las soluciones son prácticamente las mismas.

Simulación de la infiltración en surcos en 2D: tirantes constantes en los surcos

Tomando como base las características hidrodinámicas del suelo Montecillo se realizó una simulación del proceso de infiltración en tres surcos, en donde se aplicó un tirante de 10 cm en la primera y la tercera horas de infiltración, y en la segunda y cuarta hora se aplicaron 5 cm de tirante. Se trazaron curvas con los siguientes contenidos volumétricos de agua: 0.47 cm 3 cm -3, 0.43 cm 3 cm -3, 0.39 cm 3 cm -3, 0.35 cm 3 cm -3, 0.31 cm 3 cm -3. En la figura 5 se muestran dos perfiles de humedad correspondientes a la primera hora de infiltración, en el minuto 60 el tirante cambia de 10 a 5 cm. Como puede verse, el frente de humedecimiento avanza más rápido en la primera hora (figura 5) que en la segunda (figura 6) debido a que en la segunda hora se tiene un tirante de 5 cm desde la base del surco, por lo que, en dicha hora, los efectos de la capilaridad aportan más que los efectos de presión. El efecto inverso puede observarse en la imagen superior de la figura 7, en el cual la presión hacia abajo es mayor que en los minutos anteriores, debido a que el tirante cambió a 10 cm en el minuto 120.

La figura 6 muestra dos perfiles de humedad correspondientes a la segunda hora de infiltración y la figura 7 muestra dos perfiles de humedad correspondientes a la tercera y cuarta hora. Como puede verse en la serie de imágenes, cuando se hace el cambio de tirante de 10 a 5 cm en el surco se favorece el humedecimiento en los lomos del surco, mientras que con un tirante mayor la redistribución del perfil de humedad es mayor en la base del surco que en el lomo. Con un tiempo de simulación de 60 minutos se puede apreciar que ya se ha logrado llevar el contenido de humedad a 0.31 cm3 cm-3 en una profundidad de 35 cm y los lomos del surco.

El contenido de humedad en el surco a los 180 minutos de simulación puede verse que está prácticamente saturado, es decir, tiene una cantidad de agua ligeramente mayor a 0.47 cm3 cm-3 y la máxima cantidad de agua que puede tener es de 0.4865 cm3 cm-3. Después de los 135 min, la parte superior de los surcos está saturada, y a partir de ese momento, para cualquier contenido de humedad en el perfil del suelo tiende a ser horizontal, que se asemeja a una infiltración constante (figura 7).

Simulación de la infiltración en surcos en 2D: tirantes diferentes en los surcos

Usando las mismas características hidrodinámicas del suelo Montecillo del inciso anterior se realizó una simulación en tres surcos, en donde en cada uno de ellos las cargas ahora fueron diferentes: 3, 9 y 6 cm, de izquierda a derecha. Los resultados de la simulación mostraron que la redistribución de humedad en los lomos del surco y en el perfil del suelo durante el proceso de infiltración es mejor en el área entre los surcos con cargas 9 y 6 cm (área 1), que el que está entre los surcos con 3 y 9 cm de carga (área 2) (figuras 8 y 9).

La distribución de la humedad en los lomos de los surcos, la unión de los bulbos de humedecimiento y área de mojado es mejor en la primera área que en la segunda. En el área 1, el surco que tiene 9 cm de carga tiene que aportar mayor cantidad de agua al surco de la izquierda (carga de 3 cm), por lo que la distribución de la humedad hacia capas inferiores es menor que entre los surcos de 9 y 6 cm.

Por tratarse un problema en dos dimensiones, la lámina infiltrada es ahora representada como el área infiltrada, que si se multiplica por una longitud se obtiene el volumen infiltrado. Así, el área infiltrada (cm2) se calculó en cada surco con la ecuación (17) y al hacer las comparaciones (figura 10) se puede ver que para el tiempo de 240 min, los surcos de 9, 6 y 3 cm de tirante infiltran en promedio 512, 467 y 401 cm2, respectivamente. Es decir, se estarían infiltrando 110 cm2 menos en el surco con carga de 3 cm; además de que en los surcos con mayor tirante tenderán, por diferencia en los tirantes, a contribuir al humedecimiento de los surcos con menor tirante, lo que también trae consigo una disminución en el patrón de mojado en cada uno de los surcos que les corresponde humedecer.

Conclusiones

El fenómeno de infiltración en un dominio bidimensional se modeló y resolvió mediante el método del gradiente conjugado. La naturaleza de la ecuación de Richards hizo que no fuera suficiente el método explícito ―a diferencia de otras ecuaciones no lineales―, por lo que se combinó con el método implícito. La metodología presentada mostró buenos resultados, además de que muestra la bondad de incorporar un tirante variable en el tiempo.

La arbitrariedad en la posición de los puntos del dominio espacial nos permitió crear cualquier discretización espacial; sin embargo, esto ocasiona que las derivadas tengan que ser resueltas usando la superficie que más se acerque a una región de puntos. Con el uso de esta metodología, en la zona seca del surco se puede imponer una condición de tipo Neumann, por lo que pueden aplicarse condiciones de evaporación o una función de transpiración de acuerdo con las necesidades que se requieren.

Finalmente, la solución utilizada puede ayudar a tomar decisiones sobre el tiempo de riego y la profundidad a la cual se desea mojar, ya que en las simulaciones realizadas con tiempos máximos de 4 h se pudo ver que la profundidad saturada alcanzó los 50 cm, tiempo suficiente como para mojar la zona radicular de la mayoría de los cultivos. Sin embargo, en la práctica, estos tiempos son mucho mayores, lo que trae consigo la aplicación de láminas excesivas y bajas eficiencias de aplicación.

Referencias

Brooks, R. H., & Corey, A. T. (1964). Hydraulic properties of porous media and their relation to drainage design. Transactions of the ASAE, 7 (1), 26-28. [ Links ]

Burdine, N. (1953). Relative permeability calculations from pore size distribution data. Journal of Petroleum Technology, 5(3), 71-78. [ Links ]

Dai, Y. H., & Yuan, Y. (1999). A nonlinear conjugate gradient method with a strong global convergence property. SIAM Journal on Optimization, 10, 177-182. [ Links ]

Delaunay, B. (1934). Sur la sphere vide. Izv. Akad. Nauk SSSR, Otdelenie Matematicheskii i Estestvennyka Nauk, 7(793-800), 1-2. [ Links ]

Fuentes, C. (1992). Approche fractale des transferts hydriques dans les sols no-saturés (276 pp.). Tesis de doctorado. Grenoble: Universidad Joseph Fourier de Grenoble. [ Links ]

Fuentes, C., De León-Mojarro, B., & Hernández-Saucedo, F. R. (2012). Hidráulica del riego por gravedad. Capítulo 1 (pp. 1- 60). En: Riego por gravedad. Fuentes, C., & Rendón, L. (eds.). Querétaro, México: Editorial Universidad Autónoma de Querétaro. [ Links ]

Fuentes, C., Zavala, M., & Saucedo, H. (2009). Relationship between the storage coefficient and the soil-water retention curve in subsurface agricultural drainage systems: Water Table Drawdown. J. Irrig. Drain. Eng., 135 (3), 279-285. [ Links ]

López-Avendaño, J. E., Palacios-Vélez, O. L., Fuentes-Ruíz, C., Rendón-Pimentel, L., & García-Villanueva, N. H. (1997). Bidimensional analysis on infiltration in furrow irrigation. Agrociencia, 31(3), 259-269. [ Links ]

Penrose, R. (1955). A generalized inverse for matrices. In: Mathematical proceedings of the Cambridge Philosophical Society, 51(3), 406-413. [ Links ]

Richards, L. A. (1931). Capillary conduction of liquids through porous mediums. Journal of Applied Physics, 1(5), 318-333. [ Links ]

Saucedo, H., Zavala, M., Fuentes, C., & Castanedo, V. (2013). Optimal flow model for plot irrigation. Water Technology and Sciences, 4(3), 127-139. [ Links ]

Saucedo, H., Zavala, M., & Fuentes, C. (2011). Modelo hidrodinámico completo para riego por melgas. Tecnología y ciencias del agua, 2(2), 23-38. [ Links ]

Van Genuchten, M. T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society Of America Journal, 44(5), 892-898. [ Links ]

Zavala-Trejo, M., Fuentes-Ruíz, C., & Saucedo-Rojas, H. E. (2005). Radiación no lineal en la ecuación de Richards bidimensional aplicada al drenaje agrícola subterráneo. Ingeniería hidráulica en México, 20(4), 111-119. [ Links ]

Recibido: 23 de Enero de 2017; Aprobado: 12 de Julio de 2017

Carlos Chávez, chagcalos@uaq.mx

Este es un artículo publicado en acceso abierto bajo una licencia Creative Commons