Estudios de la Zona No Saturada del Suelo Vol VII. F.J. Samper Calvete y A. Paz González, 2005

MODELIZACIÓN DEL FLUJO TRANSITORIO EN SUELOS FRANCO-LIMOSOS UTILIZANDO HYDRUS-1D

C. Rubio1, P. Llorens1 y M. Th. Van Genuchten21Instituto de Ciencias de la Tierra Jaume Almera (CSIC), Lluís Solé i Sabarís s/n, 08028 Barcelona, España,
crubio@ija.csic.es
2George E. Brown Jr. Salinity Laboratory, USDA-ARS, 450 Big Springs Rd. Riverside, CA, USA.

RESUMEN. El modelo Hydrus-1D (Simunek et al. 1998a),basado en la ecuación de Richards, se ha utilizado para
estudiar los flujos transitorios en profundidad, de una
terraza de cultivo abandonada con una cubierta vegetal
herbácea (cuencas experimentales de Vallcebre, NE
España). Con dicho modelo, se ha realizado un ejercicio decalibración y validación utilizando potenciales mátricos y
humedades volumétricas del perfil. Los resultados
obtenidos han mostrado un buen ajuste a los datos detensiometría medidos en campo cuando el perfil sehumecta, para la profundidad de -20 cm (RMSE=1kPa en
validación), mientras que para el resto del perfil se observa
una reacción más lenta del modelo. Los contenidos hídricos
simulados sobreestiman ligeramente los medidos. Cuando
la demanda evapotranspirativa es elevada, el modelo
sobreestima ligeramente los potenciales mátricos y loscontenidos hídricos, sobre todo para los niveles
superficiales del perfil. Estas diferencias pueden
relacionarse con procesos de contracción-expansión de este
tipo de suelos franco limosos.
ABSTRACT. The Hydrus-1D software package (Simunek
et al. 1998a), based on the Richards equation, was used to
study transient flow in an abandoned terraced area undergrass. The field site was located in the Vallcebre research
catchments (Eastern Pyrenees, North-Eastern Spain). A
split test was carried out to calibrate and evaluate the modelusing observed matric potential and water content profiles.
Results obtained showed a good agreement with the fielddata during wetting at 20cm depth (RMSE=1kPa for thevalidation), but a slow reaction of the simulated potentials
with depths. Observed water contents and pressure heads
were overestimated during periods of high evaporative
demand. Some of the results may have been caused by thepresence of shrinkage/swelling processes in the silt loam
soil, which were only partially accounted for in the model.
1. Introducción

Durante las últimas décadas, se han desarrollado
modelos, empíricos y teóricos, para predecir los flujos

hídricos y el transporte de solutos en la zona no saturada delsuelo. Estos modelos, todavía encuentran en su solución
ciertas dificultades tanto conceptuales como de análisis,
sobretodo a escala de evento transitorio (Abassi et al.
2004). La incertidumbre en la solución de los modelos sedebe a la inestabilidad de los flujos hídricos, causada porlas características intrínsecas del propio flujo, y por las
características estructurales del suelo (Skopp, 1981; Flühler
et al. 1996; Jarvis, 1998), como pueden ser los procesos decontracción-expansión.
Además, en las cuencas de investigación de Vallcebre
uno de los principales problemas encontrados al utilizar losmodelos hidrológicos TOPMODEL y SHETRAN (Gallart
et al. (1993, 1997) y Anderton et al. (2002a, 2002b)), hasido la simulación de la respuesta de las cuencas a losprimeros eventos de lluvia después del periodo de sequía
estival.
Dentro de este marco teórico, el modelo evaluado en este
trabajo es el Hydrus-1D (Simunek et al. 1998a), que simula
la dinámica de flujos unidimensionales solucionando laecuación de Richards (1931).
El objetivo de este trabajo es evaluar HYDRUS-1D, para
simular flujos en suelos franco-limosos que presentan
procesos de contracción-expansión (Sole et al. 1992) yvalorar la respuesta del modelo a los primeros eventos en
condiciones antecedentes secas.
2. Material y métodos

Hydrus-1D puede utilizar el modelo combinado deMualem (1976) y van Genuchten (1980) o puede utilizar laecuación de Vogel y Císlerová (1988), que modificando almodelo anterior confiere flexibilidad al sigmoide de lacurva de retención cerca del área de saturación, mejorando
teóricamente la simulación de flujos hídricos en suelos con
procesos de expansión-contracción.
El modelo se ha aplicado a un perfil de suelo (0-100 cm)
bajo prado mesófilo, en una terraza de cultivo abandonada
con especies graminoides y leguminosas, localizado en lacuenca de investigación de Can Vila (Vallcebre, NE deEspaña). La información de que se dispone son las curvas

257

Rubio et al. Modelización del flujo transitorio en suelos franco-limosos utilizando Hydrus-1D

2.2. Validaciónde retención obtenidas a partir de muestras inalteradas con
los recipientes de tensión hídrica o ‘sand box’ (pF asaturación, 0.4, 1.0 y 1.5) y la membrana a presión (pF 3.0,
3.5 y 4.2). Y la conductividad hidráulica saturada (Kfs)
determinada con permeámetro de Guelph.
La precipitación (P) se ha obtenido a partir de unpluviómetro automático de la estación meteorológica
ubicada en la cuenca; mientras que la evapotranspiración dereferencia (E0) ha sido calculada según Penman-Monteith
FAO (Smith et al. 1992), a partir de las variables
meteorológicas medidas en la misma estación (Tabla 1).
Para la medida de potencial matricial se dispone de unperfil tensiométrico, con 3 tensiómetros (Skye instruments
Ltd.) en continuo (20, 40 y 60 cm), con un registro cada 20
minutos en la parte central de la terraza. Paralelamente, en
un perfil adyacente y en un segundo situado en la parteinterna de la misma terraza (cerca del talud que forma esta
terraza con la superior), se ha medido la humedad
volumétrica una vez por semana a partir del método de
Reflectometría en el Dominio Temporal (TDR), con
mediciones a las profundidades de 0-20, 20-40 y 40-60 cm
utilizando un cable téster Tektronics 1502C.
La distribución de las raíces se ha estimado a partir deobservaciones de campo y se ha descrito con la función
Trapezoidal de Absorción de Gardner (1967), mientras quepara la simulación del efecto del estrés hídrico del suelo
sobre la absorción radicular se ha utilizado la función SShape de van Genuchten (1987).

3. Resultados

2.1. Calibración

3.1. Calibración

Para la calibración del modelo HYDRUS-1D se hautilizado el periodo comprendido entre 5 de Octubre de1996 y 28 de Junio de 1997 a paso de tiempo diario.
El perfil de suelo se ha discretizado en dos capas, desde
superficie hasta los -80 cm de profundidad que presenta una
textura franco-limosa, y desde ésta a -100 cm con una
textura franco-arcillo-limosa.
Utilizando el software RETC (van Genuchten et al. 1991)
se ha ajustado la relación θ(ψ) de los datos observados a lacurva de retención, obteniendo los parámetros θr, θs(humedad volumétrica residual y a saturación), y α y N que
son parámetros de forma de curva.
Las conductividades hidráulicas saturadas de campo (Kfs)
utilizadas en la ecuación de van Genuchten (1980) han sido
los valores medidos en los ensayos infiltrométricos,
mientras que, utilizando la ecuación de Vogel y Císlerová(1988), estos valores se han aumentado 5 veces con objetode mejorar la Kfs para suelos, como en nuestro caso, convalores del parámetro de curva N (van Genuchten, 1980)
comprendidos entre 1.0 y 1.3 (Luckner et al. 1989).
Se han calibrado, a partir del método de ensayo-error, ladistribución del aparato radicular, obteniendo una
profundidad máxima de raíces de 75 cm, y el coeficiente de
reducción del potencial mátrico (h50) con un valor de -35.9
kPa para el modelo de van Genuchten (1987).

La Tabla 2 muestra la eficiencia del modelo en lacalibración, que puede considerarse de aceptable para lasdos ecuaciones utilizadas, con diferencias de subestimaciónde 0.2 kPa para los niveles más superficiales, y de 0.6 kPa
en profundidad.

258

Para la validación del modelo (28 de Septiembre a 22 deDiciembre de 2003), se han utilizado las mismas
propiedades básicas e hidrodinámicas que para lacalibración. Para las condiciones iniciales de potencial
mátrico se han utilizado medidas tensiométricas a -20, -40 y-60 cm.
Se ha validado también el contenido hídrico del sueloutilizando los 2 perfiles de TDR, en la parte interna (I) y enel centro de la terraza (C).
Tabla 1. Totales de precipitación (P) y evapotranspiración de referencia(E0), utilizados en la evaluación del modelo Hydrus-1D.
Año
Nº Obs.
Ptotal
E0 total(días)(mm)(mm)
Calibración
1996-97
267
977.5
429.5
Validación
2003
86
349.3
80.9

2.3. Eficiencia del modelo
Los criterios de eficiencia utilizados han sido, el error
medio (ME) y la raíz cuadrada del error cuadrático medio(RMSE) para evaluar los potenciales matriciales simulados,y el coeficiente de determinación (r2) para la bondad delajuste de los contenidos hídricos.

Tabla 2. Evaluación de la calibración de HYDRUS-1D para el Perfil C.
h= profundidad del tensiómetro (cm), ME= error medio, RMSE= raíz
cuadrada del error cuadrático medio, VG= van Genuchten, VC= Vogel yCíslerová.
Octubre 1996-Junio 1997
Ecuación
ME (kPa)RMSE (kPa)
h-20
h-40
h-60
h-20
h-40
h-60
VG
-1.2
-0.6
0.1
2.4
1.1
0.9
VC
-1.4
-0.8
-0.5
2.5
1.2
0.7

La ecuación de van Genuchten es la que presenta una
menor dispersión de los valores simulados, siendo lasdiferencias con la ecuación de Vogel y Císlerová de 0.1 kPa
a -20 y -40 cm, y de 0.2 kPa a -60 cm.
3.2. Validación. Potencial hídrico.
La simulación utilizando la ecuación de van Genuchten(1980) (Fig. 1), para los niveles profundos (-40 y -60 cm)
presenta un importante desfase en el tiempo de respuestadel modelo a la primera humectación (25.8 mm
precipitación) después de un periodo muy seco (28/09 a
28/10), mostrando en la respuesta, diferencias de sobreRubio et al. Modelización del flujo transitorio en suelos franco-limosos utilizando Hydrus-1D

estimación de alrededor de 10 kPa, mientras que larespuesta del modelo para el nivel más superficial (-20 cm),
se ajusta aceptablemente a los datos medidos.
La estimación de E0 por HYDRUS-1D para el periodo
simulado, presentan un valor similar, con un incremento de
7 mm respecto a la evapotranspiración calculada según
Penman-Monteith FAO.
Al utilizar la ecuación de Vogel y Císlerová (1988), losvalores simulados no mejoran el ajuste a los valores
observados (Fig. 1).
La Tabla 3, muestra la eficiencia en términos de ME y
RMSE para el periodo de validación de HYDRUS-1D. El
modelo de van Genuchten presenta una relativa menor
sobre estimación en profundidad (-60 cm) y una menor
dispersión para todas las profundidades simuladas.
Tabla 3. Evaluación de la validación de HYDRUS-1D para el perfil C. h=profundidad del tensiómetro (cm), ME= error medio, RMSE= raíz
cuadrada del error cuadrático medio, VG= van Genuchten, VC= Vogel yCíslerová.
Octubre-Diciembre 2003 (86 días)
Ecuación
ME (kPa)RMSE (kPa)
h-20
h-40
h-60
h-20
h-40
h-60
VG
-0.1
0.8
0.8
1.0
2.8
2.
VC
-0.1
0.8
1.1
1.3
3.0
3.

3.3. Validación. Humedad volumétrica.
9
La Tabla 4 muestra los ajustes entre los contenidos
hídricos simulados utilizando la ecuación de 3van
Genuchten. Para ello, se han utilizado los dos perfiles dehumedad de suelo, uno en la parte interna (I) y otro en laparte central de la terraza, adyacente al perfil detensiometría utilizado (C).
Tabla 4. Coeficientes de determinación (r2) entre las humedades
volumétricas observadas y simuladas a partir del modelo de van
Genuchten (1980) para los perfiles I y C en validación, p≤ 0.01, p≤ 0.05.
Octubre-Diciembre 2003 (n= 86)

Profundidad (cm)

0-20

20-40

40-60

Perfil I

0.424

0.685

0.065

Perfil C

0.884

0.702

0.471

Tal como era de esperar, los ajustes para los dos perfiles
disminuyen su bondad en profundidad, siendo el perfil C elque muestra los coeficientes de determinación (r2) más
elevados en todos los casos.

4. Discusión y conclusiones

El desfase temporal entre, la rápida respuesta
hidrodinámica del suelo y la más lenta del modelo, seincrementa a medida que las observaciones se hacen más
profundas, presentando problemas para simular una rápida
humectación del perfil en profundidad después de unperiodo seco. Este comportamiento puede ser atribuido a laentrada de agua por vías preferenciales, que se dan ensuelos con macro poros (por ejemplo Jarvis, 1994; Jarvis yLarsson, 1998; Christiansen et al. 2004), en nuestro caso,
principalmente causados por procesos de humectacióndesecación que comportaran la aparición de grietas
observados en el campo. Según Simunek et al. (1998b y
1999) y Ross y Smettem (2000), el hecho de obtener unpeor ajuste con la ecuación de Vogel y Císlerová, puede
deberse a la menor estabilidad temporal que presenta elmodelo cuando el frente húmedo se mueve con rapidez en
profundidad, debido a la existencia de estas vías de flujo
preferencial.
Los resultados obtenidos muestran que, HYDRUS-1D
simula aceptablemente los flujos hídricos en los niveles
superficiales (-20 cm), siendo capaz de simular, tanto larápida humectación como la respuesta a la demanda
evapotranspirativa. El modelo presenta mayor dificultad
en simular los primeros eventos de lluvia después de unperiodo seco, especialmente para los niveles más
profundos, con un tiempo de reacción demasiado largo
frente a la rápida respuesta del suelo.
Se ha observado que, la humedad antecedente en el perfil
determina la respuesta del modelo, sobre todo para losniveles más profundos, donde las diferencias a la respuesta
hidrodinámica han sido más acusadas.
Los mejores ajustes de los contenidos hídricos
lógicamente se obtienen para el perfil de la zona central,adyacente al perfil tensiométrico.
Tal y como ocurre con HYDRUS-1D, la simulación
hidrológica de los primeros eventos de lluvia después delperiodo de sequía estival, es uno de los principales
problemas encontrados por Gallart et al. (1993, 1997) y
Anderton et al. (2002a, 2002b), al utilizar los modeloshidrológicos TOPMODEL y SHETRAN en las cuencas deVallcebre. Este desajuste es debido a la dificultad en
simular correctamente los eventos de transición típicos deestos ambientes sub-húmedos (Latron, 2003), que se
producen en los periodos de humectación de la cuenca.

259

Rubio et al. Modelización del flujo transitorio en suelos franco-limosos utilizando Hydrus-1D

0

2
0.6
4

E 0 (cm ·día -1 )P (cm ·día -1 )0

Observado
Simul. van Genuchten

P
E0

6

Simul. Vogel y Císlerová
1.2? (kP a)0-50

Profundidad -20 cm-100? (kP a)0-50

Profundidad -40 cm-100? (kP a)0-50

Profundidad -60 cm

22-Dec17-Dec12-Dec

7-Dec

2-Dec

27-N ov

22-N ov

17-N ov

12-N ov

7-N ov

2-N ov

28-Oct23-Oct18-Oct13-Oct8-O ct3-O ct28-S ep

-100

Días

Fig. 1. Precipitación, E0 y potenciales mátricos observados (perfil C) y simulados para el periodo de validación.

Agradecimientos. Este trabajo ha sido posible gracias a una beca deinvestigación y Postgrado F.P.U. del Ministerio de Educación Ciencia y a
los proyectos PIRIHEROS (REN2003-08678/HID) y CANOA (CGL200404919-C02-01/HID), así como al acuerdo CSIC-DGCONA (RESEL). Los
autores agradecen la colaboración de: Dr. Francesc Gallart, Dr. Jérôme
Latron, Rafael Poyatos y Núria Martínez-Carreras.

Referencias
Abbasi, F., Feyen, J. y van Genuchten, M.Th. 2004. Two-dimensional
simulation of water flow and solute transport below furrows: model
calibration and validation. J. Hydrol., 290, 63-79 pp.

260

Anderton, S., Latron, J. y Gallart, F. 2002a. Sensitivity analysis and multiresponse, multi-criteria evaluation of a physically-based distributedmodel. Hydrological Process., 16, 333-353 pp.
Anderton, S., Latron, J., White, S., Llorens, P., Salvany, M.C. y Gallart, F.
2002b. Internal validation of a physically-based distributed modelusing data from a Mediterranean mountain catchment. Hydrol. Earth
Syst. Sci., 6, 67-83 pp.
Christiansen, J.S., Thorsen, M., Clausen, T., Hansen, S. y Refsgaard, J.C.
2004. Modelling of macropore flow and transport processes atcatchment scale. J. Hydrol., 299, 136-158 pp.
Flühler, H., Durner, W. y Flury, M. 1996. Lateral solute mixing processes
a key for understanding field-scale transport of water and solutes.
Geoderma, 70, 165-183 pp.
Gallart, F., Latron, J. y Llorens, P. 1993. Analysis of the hydrological roleof old agrocultural hillslope terraces using TOPMODEL concepts. In:

Rubio et al. Modelización del flujo transitorio en suelos franco-limosos utilizando Hydrus-1D

Robinson, M. (Ed.), Methods of hydrologic basin comparison, report
120, Institute of Hydrology, Wallingford, 131-136 pp.
Gallart, F., Latron, J., Llorens, P. y Rabadà, D. 1997. Hydrological
functioning of Mediterranean mountain basins in Vallcebre, Catalonia:
Some challenges for hydrological modeling. Hydrol. Process., 11,1263-1272 pp.
Gardner, W.R. 1967. Water uptake and salt distribution patterns in saline
soils. In: Isotope ans radiation techniques in soil physics and irrigation
studies. Int. Atomic Energy Agency Proc. Ser., 335-341 pp.
Gerke, H.H. y van Genuchten, M. Th. 1993. A dual-porosity model for
simulating the preferential movement of water and solutes instructured porous media. Water Resour. Res., 29, 305-319 pp.
Gwo, J.P., Jardine, P.M., Wilson, G.V. y Yeh, G.T. 1995. A multiple-poreregion concept to modelling mass transfer in subsurface media. J.
Hydrol., 164, 217-237 pp.
Jarvis, N.J. 1994. The MACRO model (Version 3.1). Technnicaldescription and sample simulations. Reports and Dissert, 19, Dept.
Soil Sci., Swedish Univ. Agric. Sci., Uppsala, Sweden, 51pp.
Jarvis, N.J. y Larson 1998. Modeling the impact of preferential flow on
nonpoint source pollution. In: Selim, H.M.; Ma, L. (Eds), Physicalnonequilibrium in soils: Modeling and application, Ann Arbor Press,
Chelsea, MI, 195-221 pp.
Latron, J. 2003. Estudio del funcionamiento hidrológico de una cuenca
Mediterránea de montaña (Vallcebre, Pirineo catalán). Tesis Doctoral.
Facultat de Geologia. Universitat de Barcelona. 269 pp.
Luckner, L., van Genuchten, M.Th. y Nielsen, D.R. 1989 A consistent setof parametric models for the two-phase flow of immisceible fluidsin thesubsurface. Water Resour. Res., 25, 2187-2193 pp.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivityof unsaturated porous media. Water Resour. Res., 12, 513-522 pp.
Pruess, K. y Wang, J.S.Y. 1987. Numerical modelling of isothermal andnon-isothermal flow in unsaturated fractured rock. A review. In:
Evans, D.D. y Nicholson, T.J. (Eds), Flow and transport throughunsaturated fractured rock, Geophysics monograph, 42. American
Geophysical Union, Washington, DC, 11-22 pp.

Richards, L.A. 1931. Capillary conduction of liquids through porous
mediums. Physics, 1, 318-333 pp.
Ross, P.J. y Smettem, K.R. 2000. A simple treatment of physicalnonequilibrium water flow in soils. Soil Sci. Soc. Am. J., 64, 19261930 pp.
Simunek, J., Huang, K., Sejna, M. y van Genuchten, M.Th. 1998a. The
HYDRUS-1D software package for simulating the one-dimensional
movement of water, heat and multiple solutes in variable-saturated
media. Version 1.0, IGWMC-TPS-70, International Ground Water
Modeling Center, Colorado School of Mines, Golden, Colorado, 162
pp.
Simunek, J., Sejna, M. y van Genuchten, M.Th. 1998b. The HYDRUS-2D
software package for simulating two-dimensional movement of water,
heat and multiple solutes in variably saturated media. Version 2.0,
IGWMC-TPS-53, International Ground Water Modeling Center,
Colorado School of Mines, Golden, Colorado, 251 pp.
Skopp, J. 1981. Comment on Micro-meso and macroporosity of soil. Soil
Sci. Sco. Am. J., 45, 1246 pp.
Smith, M., Allen, R.G., Monteith, J.L, Perrier, A., Pereira, L. y Segeren,
A. 1992. Report of the expert consultation on procedures for revisionof FAO guidelines for prediction of crop water requirements. United
Nations-Food and Agriculture Organization, Rome, 54 pp.
Sole, A., Josa, R., Pardini, G., Aringhieri, R., Plana, F. y Gallart, F. 1992 .
How mudrock and soil physical properties influence badlandformation at Vallcebre (Pre-Pyrennes, NE Spain). Catena, 19: 287-300
pp.
Van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J., 44,892-898 pp.
Van Genuchten, M.Th. 1987. A numerical model for water and solute
movement in and below the root zone. Reseach report nº 121, U.S.
Salinity Laboratory, USDA, ARS, Riverside, California.
Van Genuchten, M. Th., Leij, F.J y Yates, S.R.. 1991.
The RETC Code for Quantifying the Hydraulic Functions of
Unsaturated Soils, Version 1.0. EPA Report 600/2-91/065, U.S.
Salinity Laboratory, USDA, ARS, Riverside, California.
Vogel, T. y Císlerová, M. 1988. On the reliability of unsaturated hydraulicconductivity calculated from the moisture retention curve. Transport
in Porous Media, 3, 1-15 pp.

261