Solución en diferencias finitas de la ecuación de Boussinesq con porosidad drenable variable y condición de radiación fractal en la frontera

Finite difference solution of the Boussinesq equation with variable drainable porosity and fractal radiation boundary condition

Carlos Chávez^1*, Carlos Fuentes¹, Manuel Zavala², Felipe Zataráin³

¹ Facultad de Ingeniería, Universidad Autónoma de Querétaro. 76010. Cerro de las Campanas, Santiago de Querétaro, Querétaro, México. (chagcarlos@gmail.com) (cbfuentesr@gmail.com). *Autor responsable.

² Facultad de Ingeniería, Universidad Autónoma de Zacatecas. 98000. Jardín Juárez Núm. 147, Centro Histórico. Zacatecas, México. (mzavala73@yahoo.com.mx).

Recibido: mayo, 2010.
Aprobado: septiembre, 2011.

Resumen

El drenaje subterráneo es utilizado para eliminar excedentes de agua en la zona radical y suelos salinos para lixiviar las sales. La dinámica del agua es estudiada con la ecuación de Boussinesq, sus soluciones analíticas son obtenidas asumiendo que la transmisibilidad del acuífero y la porosidad drenable son constantes y que la superficie libre se abate de manera instantánea sobre los drenes. La solución en el caso general requiere de soluciones numéricas. Se ha mostrado que la condición de frontera en los drenes es una condición de radiación fractal y la porosidad drenable es variable y relacionada con la curva de retención de humedad, y ha sido resuelta con el método del elemento finito, que en un esquema unidimensional puede hacerse equivalente al método de diferencias finitas. Aquí se propone una solución en diferencias finitas de la ecuación diferencial considerando la porosidad drenable variable y la condición de radiación fractal. El esquema en diferencias finitas propuesto ha resultado en dos formulaciones: en una aparecen de manera explícita la carga y la porosidad drenable, variables ligadas con una relación funcional, que se ha denominado esquema mixto; en la otra aparece sólo la carga hidráulica, denominada esquema en carga. Los dos esquemas coinciden cuando la porosidad drenable es independiente de la carga. Los esquemas han sido validados con una solución analítica lineal, y para la no linealidad se ha mostrado que la convergencia numérica es estable y concisa. La solución numérica es útil para la caracterización hidrodinámica del suelo a través de una modelación inversa, y para un mejor diseño de los sistemas de drenaje agrícola subterráneo ya que las hipótesis consideradas en las soluciones clásicas han sido eliminadas.

Palabras clave: formulación mixta, formulación en carga, curva de retención, modelación inversa.

Abstract

Key words: mixed formulation, head formulation, retention curve, inverse modeling.

Introducción

La dinámica del agua en los sistemas de drenaje ha sido estudiada con la ecuación de Boussinesq de los acuíferos libres. En su forma unidimensional, es una de las bases para construir soluciones analíticas aproximadas de la dinámica del agua en un sistema de drenaje en régimen permanente o transitorio (Hooghoudt, 1940; Dumm, 1954). Estas soluciones no representan adecuadamente las condiciones reales, ya que han sido construidas bajo supuestos como: porosidad drenable y transmisibilidad del acuífero son constantes y que la superficie libre se abate de manera instantánea sobre los drenes. Sin embargo, considerando condiciones más representativas conduce a dificultades analíticas, razón por la cual la utilización de métodos numéricos es necesaria para construir soluciones de la ecuación de Boussinesq.

En una línea de investigación Zavala et al. (2004, 2007), basados en los conceptos de la geometría fractal y en experiencias de drenaje, recomiendan una condición de radiación fractal, la cual incluye la radiación lineal utilizada por Fuentes et al. (1997). En cuanto a la porosidad drenable, Fuentes et al. (2009), basados en los conceptos de lámina drenable y lámina drenada así como en experiencias de drenaje, proponen una expresión analítica en la cual interviene la curva de retención de humedad de los suelos.

Existe un esquema en diferencias finitas, basado en el esquema de Laasonen, propuesto por Zataráin et al. (1998) para resolver numéricamente la ecuación de Richards aplicada al fenómeno de la infiltración del agua en los suelos, con excelentes resultados. Aparte de su alta precisión, estabilidad y convergencia, el esquema tiene la ventaja adicional de su naturaleza intuitiva ya que está basado en un balance local de masa. Este esquema puede ser utilizado para resolver la ecuación unidimensional de Boussinesq del drenaje agrícola. De esta manera, el presente trabajo tiene como objetivo resolver la ecuación unidimensional de Boussinesq con el método de diferencias finitas, considerando que la porosidad drenable es variable y que las condiciones de frontera en los drenes son de radiación fractal.

Ecuación de Boussinesq

En el estudio de la dinámica del agua en sistemas de drenaje con la ecuación de Boussinesq, se asume generalmente que las variaciones de la carga hidráulica a lo largo de los tubos de drenaje (dirección y) son despreciables con respecto a las variaciones de la carga en un corte transversal (dirección x). Así, la ecuación a resolver en el dominio mostrado en la Figura 1 se reduce a una dimensión, que resulta de la ecuación de continuidad ((vH)/t + (Hq)/x = R_w) y de la ley de Darcy (q = —K_s H/x) :

donde, H = H (x,t) es la carga hidráulica contada a partir de un estrato impermeable, y es una función de la coordenada horizontal (x) y del tiempo (t); q es el flujo de Darcy o caudal por unidad de área; K_s es la conductividad hidráulica a saturación; v = v(H) es la porosidad drenable como una función de la carga; R_w es el volumen de recarga en la unidad de tiempo por unidad de área del acuífero. El caudal unitario de agua Q_u = Hq es proporcionado por:

donde T(H) es la transmisibilidad.

La capacidad de almacenamiento está definida por:

donde W = vH es la lámina de agua drenable. La igualdad μ = v se da cuando la porosidad drenable es independiente de la carga.

Una expresión de la capacidad de almacenamiento es dada por Fuentes et al. (2009):

donde θ_s es el contenido de humedad a saturación; y θ(H — H_s) representa la evolución del contenido de humedad en la posición z = H_s mientras la superficie libre desciende.

La porosidad drenable se deduce a partir de la igualdad de las ecuaciones (3) y (4):

donde es la variable de integración.

La porosidad drenable

Para calcular la capacidad de almacenamiento y la porosidad drenable es necesario proporcionar la curva de retención de humedad del suelo. En la literatura es bastante común representarla con la ecuación de van Genuchten (1980): , donde θ_r es el contenido residual de humedad, ψ_d es un parámetro de escala de la presión, m y n son dos parámetros de forma positivos. Introduciendo esta ecuación en las ecuaciones (4) y (5) proporciona la capacidad de almacenamiento y la porosidad drenable siguientes:

La porosidad drenable no tiene una forma analítica cerrada, pudiendo ser calculada mediante integración numérica. Una forma cerrada puede ser construida a partir de la difusividad de Fujita (1952) y de la relación conductividad hidráulica–difusividad hidráulica de Parlange et al. (1982) (Fuentes et al., 1992), definidas por:

donde Θ = (θ – θ_r) / (θ_s – θ_r) es un grado efectivo de saturación; α y β son parámetros de forma adimensionales tales que 0<α<1 y 0<β<1; λ_c es la escala de Bouwer (1964). De la definición de la difusividad hidráulica D(θ) = K(θ) dψ / dθ, considerando la condición θ = θ_s cuando ψ = 0, se deduce:

donde ψ_c = —λ_c.

La porosidad drenable se obtiene de las ecuaciones (5) y (10):

Se debe notar que la función θ(ψ) es implícita en la ecuación (10), y en consecuencia la función v(H). Estas funciones pueden ser explicitadas en función de la presión si se acepta α = β; en este caso la conductividad en función de la presión corresponde a la ecuación de Gardner (1958) K(ψ) = K_s exp (ψ / λ_c) ampliamente utilizada en estudios teóricos. La curva θ(ψ) correspondiente es la siguiente:

De la ecuación (4) se obtiene la capacidad de almacenamiento:

y de la ecuación (11) la porosidad drenable:

La lámina drenada, de acuerdo con Fuentes et al. (2009) es:

El contenido de humedad a saturación puede ser asimilado a la porosidad total (Φ), la cual es estimada a partir de la densidad total del suelo seco (ρ_t) y de la densidad de las partículas (ρ_o) con la fórmula Φ = 1 — ρ_t/ρ_o; el contenido de humedad residual puede ser asumido igual a cero.

Condiciones inicial y de frontera

La carga hidráulica contada a partir del estrato impermeable H(x,t) está relacionada con la carga h(x,t) contada a partir de los drenes, de acuerdo con la Figura 1, por:

donde D_o es la altura de los drenes a partir del estrato impermeable. La variación transversal de h al inicio del proceso de drenaje es considerada como la condición inicial:

En cuanto a las condiciones de frontera o condiciones en los drenes ubicados en x=0 y x=L, se han asumido formas diversas: abatimiento instantáneo en los drenes o condición de Dirichlet (Dumm, 1954), flujo proporcional a la carga o condición de radiación lineal (qh) (Fuentes et al., 1997; Fragoza et al., 2003), y de radiación fractal (Zavala et al., 2004, 2007). La condición de radiación fractal se escribe como sigue:

donde el signo positivo corresponde al dren posicionado en x= 0; el negativo al posicionado en x= L; q_s es el flujo correspondiente a h_s y depende de las características de la interfaz suelo–dren. En cuanto al exponente s, los autores argumentan que está definido por s=D/E, donde D es la dimensión fractal efectiva de la interfaz suelo–dren y E=3 es la dimensión de Euclides del espacio físico. La relación entre s y la porosidad efectiva de la interfaz es proporcionada por la ecuación presentada por Fuentes et al. (2001): (1—Φ)^s + Φ^2s =1 . La ecuación (18) contiene como casos particulares la condición de radiación lineal cuando s=1/2 y una condición de radiación cuadrática cuando s= 1. El gasto de agua que fluye a través de la frontera por unidad de longitud de dren en un sistema de drenes paralelos se obtiene de las ecuaciones (2), (16) y (18).

La evolución temporal de la lámina drenada se calcula con la siguiente expresión:

donde es una variable de integración.

Materiales y Métodos

La adaptación de la ecuación de Richards a la ecuación de Boussinesq requiere de la discretización del dominio como se muestra en la Figura 2. Para plantear la resolución numérica de la ecuación (1), se introducen los parámetros de interpolación definidos por:

tales que 0<γ<1 y 0<ω<1; i=1,2,... y j=1,2,... son los índices para el espacio y el tiempo respectivamente.

La variable dependiente (H) en un nodo intermedio i+γ para todo j se estima como:

mientras que en el tiempo intermedio j+ω para todo i se estima como:

La ecuación de continuidad se aplica en el tiempo t_j+ω la derivada temporal puede ser discretizada de acuerdo con las dos formulaciones siguientes:

Para identificar en lo sucesivo a los dos esquemas numéricos resultantes, el primero es denominado esquema mixto y el segundo esquema en la carga, ya que en el primero aparecen de manera explícita la carga y el volumen de agua que se drena, mientras que en el segundo aparece explícitamente sólo la carga; las dos formulaciones coinciden cuando la porosidad drenable es independiente de la carga y la formulación en la carga no requiere de la integración numérica en la ecuación (7), para calcular la porosidad drenable. La discretización de la derivada espacial alrededor del nodo i–ésimo es la siguiente:

De la ecuación (2) se obtiene el caudal unitario definido en los nudos intermedios:

Las cargas en los diferentes nodos y en el tiempo intermedio se obtienen de la ecuación (23), las cuales son introducidas en las ecuaciones (27) y (28) y éstas a su vez en la ecuación (26). Luego, las ecuaciones (24) y (26) se llevan a la ecuación de continuidad y se asocian términos semejantes, resultando el sistema de ecuaciones algebraicas siguiente:

donde

Para el esquema en la carga, ecuación (25), los coeficientes B_i y E_i deben ser redefinidos reemplazando en las ecuaciones (31) y (33) y por . Una vez especificadas las condiciones iniciales y de frontera, la ecuación (29) puede ser resuelta de manera eficiente mediante el algoritmo de Thomas (ver Zataráin et al., 1998). Se debe notar que el sistema no es lineal, puesto que los coeficientes (ecuaciones 30–33) dependen de la propia solución. La resolución para cada paso de tiempo es, por lo tanto iterativa

La condición de radiación fractal

Para linealizar las condiciones de frontera, se introduce una generalización del coeficiente de conductancia, es decir:

Es necesario señalar que K depende de la propia solución, sin embargo como el proceso de solución del sistema (29) es iterativo, este parámetro se calcula en función del estimador precedente.

Selección de los incrementos espacial (Δx) y temporal (Δt) en la solución numérica

De acuerdo con Zataráin et al. (1998) la discretización del dominio se realiza de modo que el incremento x_i — x_i–1 = δx sea constante para i = 4,5... N — 2 excepto en la vecindad de los drenes, es decir para x₁ = 0 : i) x₂ — x₁ = 0.4 δx, x₃ — x₂ = 0.6 δx, Δx₁= 0.1 δx ; Δx₂ = 0.6 δx y ii) x_N = L, x_N — x_N–1 = 0.4 δx, Δx_N–1 = 0.6 δx, Δx_N = 0.1 δx. El valor de interpolación en el espacio se toma como γ = en el dominio excepto, como se puede inferir, en la primera y última celdas. En cuanto a la discretización del tiempo, dada la del espacio, se sigue el enfoque clásico de escribir las ecuaciones del movimiento en forma adimensional, válido en medios homogéneos, para obtener relaciones entre la escalas espaciales y temporales características. Introduciendo variables adimensionales en la ecuación de Boussinesq, ecuación (1), definidas como x_* = x / L, t_* = t / τ, H_* = H / H_s, μ_* = μ / v_s, R_w* = R_wL² / T_sH_s, donde v_s= v(H_s) y T_s = K_sH_s, se obtiene la misma ecuación de Boussinesq con variables con asteriscos si τ = v_sL² / T_s. Dada la naturaleza parabólica de la ecuación diferencial se define el parámetro M = (Δx_*)² / Δt_*, que se encuentra comparando la solución en diferencias finitas con soluciones analíticas. El valor del parámetro para los tiempos cortos recomendado por Zataráin etal. (1998) es del orden de M0.1.

Resultados y Discusión

Comparación de la solución con una solución analítica

Para definir valores primeros de los parámetros de interpolación en el espacio y en el tiempo (γ y ω), la solución numérica se compara con una solución analítica obtenida de la ecuación de Boussinesq en un caso particular. Esta solución ha sido construida por Fuentes et al. (1997) para una linealización de la ecuación diferencial representada por una transmisibilidad constante pero con una condición de radiación lineal en los drenes y que incluye la ecuación clásica de Glover–Dumm (Dumm, 1954), a saber:

donde τ = L² / ; α_n son los valores propios y se obtienen de las raíces positivas de f (α) = α / γ_c – γ_c / α –2cot(α) = 0 ; A_n son las amplitudes correspondientes calculadas con la expresión A_n = 2 {α_n sin(α_n) + γ_c[1–cos (α_n)]} / (α_n² +γ_c²+ 2γ_c) ; y γ_c es un coeficiente de conductancia de la interfaz suelo–dren.

Los valores utilizados para la simulación son los reportados por Fuentes et al. (1997), también utilizados por Fragoza et al. (2003), y corresponden a una parcela con drenaje subterráneo ubicada en el Distrito de Riego 076 Valle del Carrizo, Sinaloa, México: L=50 m, K_s=0.557 m/d, =0.1087 m³/ m³, =2.5065 m²/d, D_o=3.5 m, H_s= 5.0 m y /e=1.5. Para observar el efecto que tiene el parámetro M en la solución numérica, se realizaron varias simulaciones y se calculó la suma de los cuadrados de los errores (SCE) asumiendo Ax constante y variando los valores de ω y M. Los resultados se muestran en la Figura 3.

El cálculo de la SCE se realizó a intervalos de 1 día hasta un total de 60 días para evitar errores ocultos en los tiempos cortos. Se puede ver que la SCE es menor con valores de ω cercanos a 1, y con valores de M menores a 0.5. Por otra parte al disminuir el valor de ω los errores en la solución aumentan aún con valores menores de 0.5 en el parámetro M. De forma general puede apreciarse que conforme el valor de M aumenta y el valor de ω disminuye, la SCE crece de manera significativa y viceversa. De esta manera, se pudo apreciar que mientras el valor (M>0.4), en la solución se presentaban inestabilidades en las fronteras, como las que se observan en la Figura 4. En ésta, se muestra la variación del abatimiento de la superficie libre, en todo el dominio de solución y un acercamiento sobre el dren, con diferentes pasos de interpolación (ω). La discretización en el espacio se realizó con Δx=1.00 m y en el tiempo con Δt=0.01 d, que corresponden a un valor de M4.33, valor que es superior al recomendado por Zataráin et al. (1998) para los tiempos cortos.

Para las simulaciones realizadas y la mostrada en la Figura 4 se pudo apreciar que el paso de interpolación óptimo que hace que la solución numérica coincida con la solución analítica, dado un criterio de error, es ω=0.98, mismo resultado que se obtiene cuando ω= 1.00. Por otra parte, cuando el valor de M<0.4 las diferencias entre las solución numérica y la solución analítica con diferentes valores de ω son mínimas. Un ejemplo de lo anterior es el que se muestra en la Figura 5, donde se puede observar el abatimiento del perfil para dos valores diferentes de M: M4.33 (Δx=1.00 m y Δt=0.01 d) y M0.04 (Δx=0.01 m y Δt=0.0001 d).

Con el primer valor de M puede verse que aún existen diferencias con los diferentes valores de c pero son mínimas, y utilizando la última opción se ve el buen acuerdo entre la solución analítica y la solución propuesta, dado un criterio de error. Con los valores de Δx= 0.01 m y Δt= 0.0001 d obtenidos con anterioridad, se procedió a realizar una simulación por un periodo de tiempo mayor, el cual se muestra en la Figura 6, donde se aprecia el abatimiento de la superficie libre y la evolución del volumen drenado por unidad de área de suelo. Los resultados muestran que no existen diferencias significativas entre la solución analítica y la solución en diferencias finitas.

Comparación de los dos esquemas numéricos

Los esquemas mixto y en carga se comparan entre sí aceptando los valores γ =0.5 y ω = 0.98. El suelo utilizado es el caracterizado por Saucedo et al. (2003), con los valores θ_s=0.5245 cm³/cm³, θ_r =0 cm³/cm³, K_s = 0.446 m/d; los valores de los parámetros de las características hidrodinámicas son: i) para Fujita y Parlange λ_c=0.521 m y α = 0.98; ii) para van Genuchten con la restricción m=1—2/n, m = 0.066 y ψ_d = — 0.15 m. Para comparar los esquemas, se propone un distanciamiento entre drenes L=25 m y una profundidad de drenes H_s =1.5 m. Los resultados de la simulación numérica obtenidas con las dos características hidrodinámicas son mostradas en las Figuras 7 y 8. En la Figura 7 se muestra la evolución del abatimiento de la carga y la lámina drenada para tiempos de 60 y 250 d, respectivamente, utilizando las características hidrodinámicas de Fujita y Parlange, y en la Figura 8 se aprecian los resultados obtenidos con las características hidrodinámicas de van Genuchten para los tiempos mencionados con anterioridad. Puede verse que no existen diferencias aparentes entre los esquemas mixto y en carga con cada una de las características hidrodinámicas utilizadas, los resultados muestran que utilizar cualquiera de los dos esquemas conduce al mismo resultado.

Problema inverso

Para evaluar la capacidad de la condición de radiación fractal con capacidad de almacenamiento variable, se hace uso de la información experimental presentada por Zavala et al. (2003). Las características del módulo de drenaje y los parámetros del suelo usados en la simulación se muestran en el Cuadro 1.

Las características hidrodinámicas utilizadas son las de Fujita y Parlange, para lo cual se fija el valor de α = 0.98 y los parámetros λ_c y κ₀ se estiman mediante la minimización de la suma de los cuadrados de los errores entre la lámina drenada medida y la lámina drenada calculada con la solución numérica en el transcurso del tiempo. En la Figura 9 se presentan los resultados obtenidos del problema inverso para 24 h y 240 h.

Se usaron los dos esquemas numéricos y el resultado fue el mismo, sin embargo, sólo se muestra el esquema en carga con la condición de radiación fractal. La mejor aproximación a los datos experimentales es obtenida con λ_c =0.5152 m y κ₀ = 0.10 que proporcionan un ECM= 0.1283 cm. La condición de radiación fractal y la capacidad de almacenamiento variable reproducen de buena manera los datos medidos en laboratorio para el intervalo de tiempo menor a 6 h, posteriormente hay una ligera subestimación de los mismos hasta un tiempo de 150 h, sin embargo, el buen acuerdo entre la lámina drenada medida y la lámina drenada obtenida con el modelo es evidente.

Conclusiones

Se ha resuelto la ecuación unidimensional de Boussinesq del drenaje agrícola con el método de diferencias finitas basada en un balance local de masa. Resultaron dos esquemas de discretización de la derivada temporal, la cual representa el cambio de almacenamiento en este balance. En uno aparecen de manera explícita la carga y la porosidad drenable, variables ligadas con una relación funcional, que se ha denominado esquema mixto; en el otro aparece sólo la carga hidráulica, denominado esquema en carga. Los dos esquemas coinciden cuando la porosidad drenable es independiente de la carga. La validación parcial de ambos esquemas fue realizada mediante la comparación de las soluciones numéricas obtenidas con una solución analítica de la literatura construida para condiciones de linealidad. Las evoluciones de la carga de agua en el perfil y la lámina drenada calculadas con la solución analítica son similares, bajo un criterio de error, a las calculadas con la solución numérica para todo tiempo. La ausencia de fluctuaciones tanto en el tiempo como en el espacio de la carga y de la lámina drenada permite recomendar los esquemas numéricos de la ecuación unidimensional de Boussinesq propuestos, para el estudio de la dinámica del agua en los sistemas de drenaje agrícola subterráneos. En particular, la solución numérica construida puede ser utilizada para la caracterización hidrodinámica del suelo a través de una modelación inversa, es decir, a partir de datos experimentales se pueden inferir los parámetros del sistema. Además, la solución numérica propuesta puede ser utilizada para un mejor diseño de sistemas de drenaje agrícola subterráneo ya que las hipótesis consideradas en las soluciones clásicas han sido eliminadas.

Literatura Citada

Boussinesq J., 1904. Recherches the'oriques sur l'e'coulement des nappes d'eau infiltre'es dans le sol et sur le de'bit des sources. J. Math. Pure. Appl. 5me. Ser. 10: 5–78.         [ Links ]

Bouwer H. 1964. Rapid field measurement of air entry value and hydraulic conductivity of soil as significant parameters in flow system analysis. Water Resourses Res. 36: 411–424.         [ Links ]

Dumm, L. 1954. Drain spacing formula. Agric. Engi. 35: 726730.         [ Links ]

Fragoza F., C. Fuentes, M. Zavala, F. Zataráin, H. Saucedo, y E. Mejía. 2003. Drenaje agrícola subterráneo con capacidad de almacenamiento variable. Ing. Hidráulica Méx. 18(3): 81–93.         [ Links ]

Fuentes C., F. Brambila, M. Vauclin, J.–Y. Parlange, y R. Haverkamp. 2001. Modelación fractal de la conductividad hidráulica de los suelos no saturados. Ing. Hidráulica Méx. 16(2): 119–137.         [ Links ]

Fuentes C., R. Haverkamp, and J.–Y. Parlange. 1992. Parameter constraints on closed–form soil–water relationships. J. Hydrol. 134: 117–142.         [ Links ]

Fuentes C., R. Namuche, L. Rendón, R. Patrón, O. Palacios, F. Brambila, y A. González. 1997. Solución de la ecuación de Boussinesq del régimen transitorio en el drenaje agrícola bajo condiciones de radiación: el caso del Valle del Carrizo, Sinaloa. Hermosillo, México: VII Congreso Nacional de Irrigación. pp: 3–141 a 3–145.         [ Links ]

Fuentes C., M. Zavala, and H. Saucedo. 2009. Relationship between the Storage Coefficient and the Soil–Water Retention Curve in Subsurface Agricultural Drainage Systems: Water Table Drawdown. J. Irrigation and Drainage Eng. 135(3): 279–285.         [ Links ]

Fujita H. 1952. The exact pattern of a concentration–dependent difussion in a semi–infinite medium, part II. Textil Res. J. 22: 823–827.         [ Links ]

Gardner W. R. 1958. Some steady–state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci. 85: 228–232.         [ Links ]

Hooghoudt S. 1940. Bjidrage tot de kennis van enige natuurkundige grootheden van der grond. Verslag andbouwk Onderzoek 46(7): 515–707.         [ Links ]

Parlange J.–Y., R. D. Braddock, I. Lisley, and R. E. Smith. 1982. Three parameter infiltration equation. Soil Sci. 11: 170–174.         [ Links ]

Richards L. A. 1931. Capillary conduction of liquids trough porous mediums. Physics 1: 313–333.         [ Links ]

Saucedo H., P. Pacheco, C. Fuentes, y M. Zavala. 2003. Efecto de la posición del manto freático en la evolución del frente de avance en el riego por melgas. Ing. Hidráulica Méx. 18(4): 119–126.         [ Links ]

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

Zataráin F., C. Fuentes, V. O. L. Palacios, E. J. Mercado, F. Brambila, y N. Villanueva. 1998. Modelación del transporte de agua y de solutos en el suelo. Agrociencia 32(4): 373–383.         [ Links ]

Zavala M., C. Fuentes, y H. Saucedo. 2003. Sobre la condición de radiación lineal en el drenaje de una columna de suelo inicialmente saturado. Ing. Hidráulica Méx. 18 (2): 121–131.         [ Links ]

Zavala M., C. Fuentes, y H. Saucedo. 2004. Radiación fractal en la ecuación de Boussinesq del drenaje agrícola. Ing. Hidráulica Méx. 19(3): 103–111.         [ Links ]

Zavala M., C. Fuentes, y H. Saucedo. 2007. Nonlinear radiation in the Boussinesq equation of agricultural drainage. J. Hydrol. 332(3): 374–380.         [ Links ]