text in 
Spanish (pdf)
Article in xml format
Article references
Automatic translation
Send this article by e-mail
Cited by SciELO 
Similars in
SciELO 
Permalink
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
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:
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):
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:
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.
Figura 1 Izquierda: frontera con la parte superior parecida a un surco de riego. La parte superior f es una función sinusoidal trasladada y escalada para que concuerde con la parte superior de los lados a y c. Derecha: discretización de la frontera izquierda, puntos frontera.
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).
Figura 2. Punto Pk dentro de una frontera tipo “surco”. Obsérvese que para un punto que está en el interior de una frontera, el ángulo cubierto es de 360 grados.
Figura 3 (a): dominio tipo “surco” con 10,000 puntos interiores; (b): dominio (1) después de eliminar los puntos que tienen una distancia menor a 5 mm; (c): triangulación del dominio (2).
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:
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):
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:
Cuando C [Ψ (k)] ≠ 0, la solución que se obtiene es una solución recursiva:
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:
El gradiente conjugado no lineal para Richards 2D
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ó
Método para resolver el operador
Para simplificar la notación se hicieron las siguientes simplificaciones: se
define
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:
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
donde aplicando la misma consideración que en la solución de la ecuación (14), la matriz A queda definida como:
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.
Cálculo de la lámina infiltrada
El volumen infiltrado por unidad de longitud de cauce en la unidad de tiempo
(
donde
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
Figura 4 Comparación a los 240 minutos entre el esquema implícito (líneas blancas) y el esquema mixto (líneas negras).
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.
Figura 5 Primera hora de infiltración con un tirante de 10 cm desde la base del surco. A los 60 minutos el tirante cambia a 5 cm.
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.
Además, el contenido de humedad en el lomo del surco con carga de 9 cm es mayor que en el lomo del surco con 3 cm; en este último, el principal problema es que se retrasa el mojado hacia las capas inferiores (figura 9). Este problema es muy común en los eventos de riego por gravedad, donde el caudal aplicado en cada surco no es homogéneo debido principalmente a que los sifones no están calibrados, que hay fallas en la apertura de compuertas o, como en la mayoría de los casos, que el gasto proporcionado en cada surco es a criterio del regador. Lo anterior trae consigo varias complicaciones, entre la más importante está que la lámina aplicada en cada surco no es homogénea, lo que conduce a que en algunas partes de la parcela se aplique más agua de la que necesita o, en el peor de los casos, que se aplique una cantidad menor.
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.
