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:

1

la ecuación de momentum se escribe en la forma recomendada por ^{Saucedo, Fuentes y Zavala (2005)}:

2

donde q(x,t) = U(x,t)h(x,t) es el gasto por unidad de ancho de melga [L²T^-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]; J_o, la pendiente topográfica de la melga [LL^-1]; J, la pendiente de fricción [LL^-1]; V_I = 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 β = U_lx/U, siendo α = 1-U¡_x / U, donde U_Ix 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:

3

donde ν es el coeficiente de viscosidad cinemática del agua [L²T^-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:

4

5

donde x_f(t) es la posición del frente de onda para el tiempo t y q_o 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: ψ = -h_f < 0, donde h_f 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 = K_s, 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:

6

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:

7

La ecuación de momentum guarda la misma forma discreta para las cuatro fases del riego:

8

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 q- = (1 - φ)q_m + φq_j,, h- = (1 - φ) h_m + φh_j (^{Saucedo et al., 2005}). Se introducen las variables pequeñas (^{Strelkoff & Katopodes, 1977}): δh_r, δq_r, δh_l y δq_l, tales que para las celdas interiores se tiene: h_l, = h_j. + δh_l, h_r = h_m + δh_r q_l = q_j. + δq_l y q_r = q_m + δq_r . Se sustituyen en las formas discretas de las ecuaciones de continuidad y momentum y se obtiene un sistema de ecuaciones algebraicas, cuya solución permite avanzar en el tiempo los valores de las variables hidráulicas. El detalle de las entradas al sistema matricial puede consultarse en ^{Saucedo et al. (2011)}.

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)}: Δt_mín = 5 s, ^{Akanbi y Katopodes (1988)}: Δt_máx = 1 s, ^{Playán, Walker y Merkley (1994)} Δt_mín = 2.12 s.

Optimización de parámetros

Para estimar los parámetros hidrodinámicos K_s y h_f , 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:

9

donde J(p) es la matriz Jacobiana relacionada con las variaciones de la función de posición del frente de avance x_f = x_f (t) respecto a cada uno de los parámetros a optimar (K_s y h_f ); 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:

10

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 x_f = x_f(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 x_f = x_f(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 q_o = 0.0032 m³/s/m; pendiente topográfica J_o = 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.2749 cm³/cm³. El valor de la porosidad se asimila al contenido volumétrico de agua a saturación, de forma tal que θ_ί, = 0.4865 cm³/ cm³. En la Figura 1 se muestra la gráfica de los valores registrados del frente de avance.

Figura 1 Valores observados en una prueba de avance reportada en la literatura. 

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 h_f = 40.4 cm y K_s = 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.

Figura 2 Resultado de la primera iteración del método Levenberg-Marquardt, h_f = 33.4 cm y K_s = 1.58 cm/h. 

Figura 3 Resultado de la tercera iteración del método Levenberg-Marquardt, h_f = 38.3 cm y K_s = 1.49 cm/h 

Figura 4 Resultado de la séptima iteración del método Levenberg-Marquardt, R² = 0.997, h_f = 40.4 cm y K_s = 1.35 cm/h. 

Figura 5 Evolución del error cuadrático medio durante la aplicación del método Levenberg-Marquardt. 

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.