Introducción
El objetivo del diseño de riego por melgas es conseguir la aplicación de una lámina de riego requerida por el cultivo de la manera más uniforme posible, conservando una eficiencia de aplicación alta. El diseño del riego consiste en determinar el gasto de aporte y el tiempo durante el cual se aplica dicho gasto en la cabecera de la melga para lograr la mayor uniformidad posible; es decir, en determinar el gasto óptimo para una longitud de melga específica y una caracterización hidrodinámica del suelo dada. Para realizar dicha caracterización hidrodinámica puede utilizarse, por ejemplo, información de una prueba de avance de riego, a partir de la cual se realiza el procedimiento de ajuste u optimación de los parámetros, buscando minimizar el error obtenido entre las observaciones realizadas y las obtenidas mediante el modelo empleado para la simulación de la fase de avance del riego por gravedad.
El objetivo del presente estudio fue desarrollar un modelo que permita hacer de forma automática la optimación de los parámetros de conductividad hidráulica a saturación del suelo y de la presión en el frente de humedecimiento, empleando una forma precisa de modelar el flujo del agua a superficie libre, para lo cual se hace uso de las ecuaciones de Saint-Venant; mientras que para la descripción del flujo del agua en el suelo, se utiliza la ecuación de Green y Ampt, la cual es una forma sencilla de modelar el fenómeno manteniendo bases físico-matemáticas en la representación.
Flujo del agua sobre la superficie del suelo
El flujo del agua con una superficie libre es modelado con las ecuaciones de Saint-Venant que resultan de la aplicación de las leyes de conservación de masa y cantidad de movimiento. En una melga, la relación entre su ancho y el tirante de agua permite considerar las ecuaciones correspondientes al escurrimiento sobre una superficie de ancho infinito (Woolhiser, 1975). La ecuación de continuidad se escribe como:
la ecuación de momentum se escribe en la forma recomendada por Saucedo, Fuentes y Zavala (2005):
donde q(x,t) = U(x,t)h(x,t) es el gasto por unidad de ancho de melga [L2T-1]; x, la coordenada espacial en la dirección principal del movimiento del agua en la melga [L]; t, el tiempo [T]; U, la velocidad media; h el tirante de agua [L]; Jo, la pendiente topográfica de la melga [LL-1]; J, la pendiente de fricción [LL-1]; VI = 31/ dt, el flujo de infiltración [LT-1], es decir, el volumen de agua infiltrado en la unidad de tiempo por unidad de ancho y por unidad de longitud de la melga; I, la lámina infiltrada [L]; g, la aceleración gravitacional [LT-2]; el parámetro adimensional β = Ulx/U, siendo α = 1-U¡x / U, donde UIx es la proyección en la dirección del movimiento de la velocidad de salida de la masa de agua debido a la infiltración.
La relación entre las variables hidráulicas q y h con la pendiente de fricción, denominada ley de resistencia hidráulica, es adoptada de acuerdo con Fuentes, De León, Saucedo, Par-lange y Antonino (2004), haciendo uso de una ley potencial de resistencia:
donde ν es el coeficiente de viscosidad cinemática del agua [L2T-1] y k es un factor adimensional; d es un parámetro adimensional que varía entre 1/2 < d < 1 en función del tipo de flujo; los valores extremos d = 1/2 y d = 1 corresponden respectivamente al régimen de Chézy y al flujo laminar de Poiseuille.
Las condiciones inicial y de frontera que deben sujetar a las ecuaciones de Saint-Venant para modelar la fase de avance son:
donde xf(t) es la posición del frente de onda para el tiempo t y qo es el gasto de aporte en la entrada de la melga.
Para cerrar el sistema es necesario conocer la forma en que evoluciona en el tiempo la lámina infiltrada en toda posición sobre la melga, es decir, la ley de infiltración.
Flujo del agua en el suelo
El modelo de Green y Ampt (1911) se establece a partir de la ecuación de continuidad y la ley de Darcy 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 el tirante de agua; c) existe un frente de humedecimiento bien definido caracterizado por una presión negativa: ψ = -hf < 0, donde hf es la succión en el frente de humedecimiento; d) la región entre la superficie del suelo y el frente de humedecimiento está completamente saturada (flujo en pistón): θ = θs y K = Ks, donde K 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 Δθ = θ s - θo es la capacidad de almacenamiento; I es el volumen infiltrado acumulado por unidad de superficie de suelo o lámina infiltrada.
Solución numérica
Para simular la fase de avance del riego por melgas se hace uso del modelo desarrollado por Saucedo, Zavala y Fuentes (2011). La forma discreta de la ecuación de continuidad para la fase de avance se escribe como:
La ecuación de momentum guarda la misma forma discreta para las cuatro fases del riego:
En las ecuaciones (7) y (8), δt es el paso de tiempo, ω y φ son factores de peso en espacio y tiempo. El cálculo de los coeficientes se realiza con base en los valores pertenecientes al nivel de tiempo anterior
La ecuación de Green y Ampt (ecuación (6)) se resuelve numéricamente usando un método de diferencias finitas; el procedimiento se encuentra bien documentado en la literatura y puede consultarse, por ejemplo, en Burden y Faires (1985).
Se ha utilizado un paso de tiempo constante Δt = 1.0 s para el acoplamiento de las ecuaciones de Saint-Venant, y Green y Ampt. La discretiza-ción utilizada para la solución de la forma completa de las ecuaciones de Saint-Venant guarda semejanza con las reportadas en la literatura: Katopodes y Strelkoff (1977): Δtmín = 5 s, Akanbi y Katopodes (1988): Δtmáx = 1 s, Playán, Walker y Merkley (1994) Δtmín = 2.12 s.
Optimización de parámetros
Para estimar los parámetros hidrodinámicos Ks y hf , a partir de los datos registrados durante una prueba de avance, puede utilizarse el método Levenberg-Marquardt, en el cual el valor estimado de los parámetros en una iteración dada se calcula teniendo en cuenta una expresión de la forma:
donde J(p) es la matriz Jacobiana relacionada con las variaciones de la función de posición del frente de avance xf = xf (t) respecto a cada uno de los parámetros a optimar (Ks y hf ); I es la matriz identidad; K(p), el vector de diferencias entre las posiciones del frente de avance observadas en campo y aquellas obtenidas con el método Levenberg-Marquardt en la iteración anterior; λ se denomina parámetro de amortiguación y es importante en la convergencia del método, su valor puede determinarse, por ejemplo, con alguno de los procedimiento sugeridos por Griva, Nash y Sofer (2009).
La matriz Jacobiana J(p) se calcula como:
en esta ecuación, i = 1, 2,..., n es el número de puntos que se registraron durante la prueba de avance, los cuales se recomienda sean diez a lo más, para reducir el esfuerzo de cómputo. La función de posición del frente de avance xf = xf(t) es resultado de la aplicación del acoplamiento numérico interno de las ecuaciones de Saint-Venant, y Green y Ampt. La estimación de las derivadas se realiza numéricamente, a partir de la función xf = xf(t) calculada para dos valores cercanos de los parámetros a optimar.
A manera de ejemplo, se toma el caso de una prueba de avance reportada en la literatura, registrada en el suelo franco de Montencillo por Fuentes (1992). Se dispone de la siguiente información: gasto unitario qo = 0.0032 m3/s/m; pendiente topográfica Jo = 0.002 m/m; longitud de la melga L = 100 m; parámetros para la ley de resistencia d = 1 y k = 1/54; parámetro en la ecuación de cantidad de movimiento β = 2; valor inicial del contenido volumétrico de agua θ0 = 0.2749 cm3/cm3. El valor de la porosidad se asimila al contenido volumétrico de agua a saturación, de forma tal que θί, = 0.4865 cm3/ cm3. En la Figura 1 se muestra la gráfica de los valores registrados del frente de avance.
En las Figuras 2 a 4 (3)se muestra la aproximación que se tiene para las iteraciones primera, tercera y séptima de la aplicación del método Levenberg-Marquardt. En la Figura 5 se muestra la evolución del error cuadrático medio conforme avanzan las iteraciones del método. Los valores de los parámetros hf = 40.4 cm y Ks = 1.35 cm/h fueron obtenidos mediante la aplicación del procedimiento de optimización para reproducir datos de la prueba de riego efectuada en el experimento mencionado.
Puede apreciarse que el número de iteraciones para lograr la convergencia es relativamente bajo, pero debe tenerse en cuenta que el mayor esfuerzo de cómputo se destina a la estimación de la matriz Jacobiana que aparece en las ecuaciones (9) y (10).
Conclusiones
Se desarrolló un método para la estimación de parámetros de infiltración basado en el empleo de las ecuaciones de Saint-Venant para describir el flujo del agua sobre el suelo, y la ecuación de Green y Ampt para representar el flujo del agua en el suelo. La estimación de los parámetros hidrodinámicos de conductividad hidráulica a saturación y presión en el frente de humedad se realiza aplicando el método LevenbergMarquardt. El modelo así obtenido permite el ajuste de los parámetros hidrodinámicos a partir de datos de pruebas de avance de riego por melgas y de la textura del suelo.