Evaluation of various infiltration models

Evaluation of various infiltration models

Many approaches have been presented to solve the problem of infiltration of water into a soil. Youngs (1995) reviews the historic development of infiltration theory including the classic solutions based on the Richards equation. For a comprehensive review of analytical and empirical solutions see Kutilek & Nielsen (1994). A typical use of these parameterized solutions is as fitting models for measured data sets. The optimized parameters serve as a convenient, condensed description of the data and are also used for predictive purposes.

Uncertainty is inherent in model predictions and is caused by several factors including imperfectly known model parameters and potential model error, e.g., errors due to violation of assumptions made in the development of the model. It would seem desirable to report, as a rule, an estimate of uncertainty together with any model prediction. While the model error will in general be unknown, the parameter uncertainty can be estimated in form of confidence intervals if the parameters are obtained by fitting measured data. This estimate is derived from the shape of the objective function in the immediate vicinity of the optimized parameter set and will be precise for a linear fitting model with zero model error and measurement errors that are uncorrelated and normally distributed. Using the Monte-Carlo method, it is possible to obtain the actual distributions of optimized parameters that can be used to evaluate the accuracy of estimated confidence limits for nonlinear fitting models.

We compared infiltration equations in terms of their fitting ability and accuracy of estimated parameter confidence intervals. Specifically, we tested six models for three different levels of simulated measurement errors. The objective was to identify the models that provide the best parameter and uncertainty estimates while converging fast and reliably.

Selection of infiltration equations was necessarily somewhat arbitrary given their great number, however, we included those that we believe to be widely used. Empirical models tested were those proposed by Mezencev (1948) and by Horton (1940). Four models derived from physical considerations were tested: Green & Ampt (1911), Parlange et al. (1982), Philip (1957), and Swartzendruber (1987). Within the framework of assumptions made in the development of these models, their parameters can be physically interpreted. However, when applied to field-obtained infiltration data, the physically based models are used in a manner effectively similar to the empirical models because a number of the assumptions are typically violated (e.g., uniformity of medium properties and of initial water content). Thus for the purposes of this study, each model's parameter vector was considered unknown and subjected to the fitting procedure.

Reference data free of measurement errors are prerequisite to generating data sets with controlled 'measurement' errors. For this purpose, cumulative infiltration I as a function of time t was predicted using a one-dimensional numerical scheme solving the Richards equation. Assumptions inherent in the simulations include a rigid homogeneous medium, isothermal conditions, and a continuous air phase at atmospheric pressure at all times.

The hydraulic properties of a given soil determine the shape of the respective infiltration curve, possible shapes ranging from significantly curved over the time frame of interest to quickly approaching steady state. To cover two points in the soil hydraulic spectrum, reference data were obtained for a sandy loam and for a clay. For each, two I(f) reference data sets over 5 hr and 20 hr were obtained and represented by 70 and 140 (I, t) points, respectively. Points were selected so as to be equal-spaced on a t^0.5 -axis, thus providing the highest density at the origin.

Fitting data sets over two different periods (5 hr and 20 hr) was done to identify the effect, if any, of the later part of the curve on goodness of fit and on the length of confidence intervals for optimized parameters. In addition, it allowed comparison of optimized parameters between long and short sets to identify time-dependence of the parameters.

Soil water retention and conductivity curves were described by van Genuchten's (1980) relationships, Equations (1) through (3). Using data given by van Genuchten (1980), q_r, q_s, a, n, K_s, and l, were set equal to 0.0, 0.446, 0.152 m^-1, 1.17, 3.417 10^-5 m/hr, and 0.5, respectively, for the clay. In case of the sandy loam, the hydraulic parameters were set to 0.1346, 0.3213, 1.74 m^-1, 1.8646, 3.5125 10^-3 m/hr, and -0.4509, respectively (Wendroth, 1993). Initial condition was equilibrium with a negative soil-water pressure head of 5 m and 15 m at the soil surface for the sandy loam and clay, respectively. Constant-head boundary conditions were applied at the top and bottom of the modeled spatial domain. Surface ponding depth was set to 5 cm for both soil types. Pressure head at the bottom was fixed at zero, simulating a stationary groundwater table. For consistency with the equilibrium initial condition, bottom depth was 15 m for the clay and 5 m for the sandy loam.

The solution h(z, t) was obtained from a finite-element, Picard time-iterative model using linear elements and a Galerkin formulation at each time step. An adaptive-grid algorithm was developed that periodically recreates the finite-element grid, depending on the position of the wetting front, to keep the highest node density at the critical portion of the h(z)-curve. The number of grid nodes may change over the course of a simulation. For the cases presented here, it was on the order of 1000 at all times. Total mass balance errors were 1.62% and 0.69% for the sand and clay simulations, respectively. The adaptive-grid finite-element Richards-Equation Solver RES-ID, written in C for UNIX, will be documented and released into the public domain.

Six infiltration equations with either two or three parameters were tested. The model equations are listed in TABLE 1, together with respective two-letter abbreviations.

Note that both the GA and the PA model are nonexplicit expressions requiring iterative solutions in t for each / to be predicted. Since cumulative infiltration grows monotonically with time, a bisection method was employed in both cases to find the value of I whose argument t matches the time for which the respective prediction is to be made. The iterative nature of the bisection process causes an increase in computational cost by one or several orders of magnitude, depending on the desired accuracy.

The six-parameter infiltration model by Haverkamp et al. (1990, implicit equation; recently, an explicit approximation has been published by Barry et al., 1995) was also tested but found to be poorly suited for fitting. The model was excluded from this study because the nonuniqueness problems resulting from the large number of parameters prevented effective parameter estimation.

Parameter optimization was performed by squared-residual minimization using the program LM-OPT (Clausnitzer & Hopmans, 1995), an implementation of the Levenberg-Marquardt (LM) algorithm. The LM algorithm combines the standard steepest-descent and quadratic-extrapolation minimization methods using a unitless parameter I that is updated after each iteration. The l-updating procedures in LM-OPT are extensions of the standard LM procedure, aimed at improving efficiency by reducing the number of iterations needed to reach convergence at an optimum. LM-OPT was used to minimize an objective function q defined as the sum of normalized residuals according to

where t_i, I_i, r_i, and s²_i are the time and cumulative infiltration, residual, and measurement variance, respectively, at data point T. The total number of data points is N. The predicted cumulative infiltration is obtained from the model function 7(t_i, p) utilizing the parameter vector p, which is different for each of the infiltration equations of TABLE 1. Each optimization was terminated once there was less than 0.001 relative improvement in the objective function q.

'Measured' data sets were generated from each of the reference sets by imposing normally distributed, independent 'measurement' errors in I on the I(t) reference data. To identify the effect of the magnitude of the measurement errors, three values for the standard deviation of the imposed errors, s_E, were chosen for each soil. These error levels S_E were defined so that for any two adjacent data points I(t_i) and I(t_i+1), I(t_i) + ks_E would be approximately equal to I(t_i+1) - ks_E, with k equal to three, two, or one. The corresponding s_E-values are 0.000007 m, 0.000014 m, and 0.000021 m for the clay, and 0.0001 m, 0.0002 m, and 0.0003 m for the sandy loam. Five thousand data sets were generated for each reference set and error level. Setting variance values s²i equal to the square of the respective s_E-value, each model equation was fitted to each of the error-imposed data sets. Thus, 5000 optimized parameter sets were obtained for each model equation and given combination of reference data and error level. Initial parameter estimates for fitting the error-imposed data sets were obtained by fitting each model equation to the respective reference curve first.

REFERENCES

BARRY D.A.; PARLANCE, J.Y.; HAVERKAMP, R.; ROSS, P.J. Infiltration under ponded conditions: 4. An explicit predictive infiltration formula, Soil Science, v. 160, p.8-17,1995.

CLAUSNITZER, V.; HOPMANS, J.W. Nonlinear parameter estimation: LM-OPT, General purpose optimization code based on the LcvenDerg-Marcjuardt algorithm, Version 1.0 for UNIX, Land, Air and Water Resources Paper 100032, University of California, Davis, 1995.

GREEN, W.A.; AMPT, G.A. Studies on soils physics: 1. The flow of air and water through soils, Journal of Agricultural Science, v.4, p. 1-24,1911.

HAVERKAMP, R., J.Y. PA, J.L. STARR, G. SCHMITZ, AND C. FUENTES. Infiltration under ponded conditions: 3. A predictive equation based on physical parameters, Soil Science, v.149, p.292-300, 1990.

HORTON, R.E. An approach towards a physical interpretation of infiltration capacity, Soil Science Society of America Proceedings, v.5, p.399-417, 1940.

KUTILEK, M.; NIELSEN, D.R. Soil Hydrology, p. 140-159, Catena-Verlag, 1994.

MEZENCEV, V.J. Theory of formation of the surface runoff [Russian], Meteorologia e Hidrologia, v.3, p.33-40, 1948.

PARLANCE, J.Y.; LISLE, L; BRADDOCK, R.D.; SMITH, R.E. The three-parameter infiltration equation, Soil Science, v.133, p.337-341, 1982.

PHILIP, J.R, The theory of infiltration: 4. Sorptivity and algebraic infiltration equations, Soil Science, v.84, p.257-264, 1957.

SWARTZENDRUBER, D. A quasi-solution of Richards' equation for the downward infiltration of water into soil, Water Resources Research, v.23,p.809-817, 1987.

VAN GENUCHTEN, M.Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Science Society of America Journal, v.44, p.892-898, 1980.

WENDROTH, O.; EHLERS, W.; HOPMANS, J. W.; KAGE, H.; HALBERTSMA, J.; WOSTEN, J.H.M. Reevaluation of the evaporation method for determining hydraulic functions in unsaturated soils, Soil Science Society of America Journal, v.57, p. 1436-1443,1993.

YOUNGS, E. The physics of infiltration, Soil Science Society of America Journal, v.59, p 307:313, 1995.

Recebido para publicação em 16.04.97

Aceito para publicação em 06.05.97

Publication Dates

• Publication in this collection
  31 May 2005
• Date of issue
  June 1997