HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (13) Citing Articles via Google Scholar Google Scholar Articles by Bitterlich, S. Articles by Knabner, P. Search for Related Content PubMed Articles by Bitterlich, S. Articles by Knabner, P. GeoRef GeoRef Citation Agricola Articles by Bitterlich, S. Articles by Knabner, P. Related Collections Vadose Zone Processes and Chemical Transport

ORIGINAL RESEARCH

Inverse Estimation of the Unsaturated Soil Hydraulic Properties from Column Outflow Experiments Using Free-Form Parameterizations

^a Institute of Applied Mathematics, Univ. of Erlangen-Nuremberg, Martensstr. 3, 91058 Erlangen, Germany
^b Institute of Geoecology, Braunschweig Technical Univ., Langer Kamp 19c, 38106 Braunschweig, Germany
^* Corresponding author (w.durner{at}tu-bs.de)

Received 27 November 2003.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY APPLICATION RESULTS DISCUSSION CONCLUSIONS REFERENCES

 
Inverse methods are increasingly used to estimate the hydraulic properties of unsaturated soils. The method generally uses a weighted least-squares approach in which numerically simulated data are fitted to measured data. In this study we used inverse methods to estimate the unsaturated soil hydraulic properties from continuous and multistep column outflow experiments. The method employs piecewise polynomial functions to obtain a free-form parameterization of the hydraulic properties, rather than fixed functional forms typical of the van Genuchten–Mualem and Brooks–Corey–Burdine models. For the polynomial functions we used quadratic B-splines and piecewise cubic Hermite interpolation. The method leads to local parameterizations that can also be hierarchic, depending on the invoked number of degrees of freedom. Since a suitable number of degrees of freedom cannot be defined a priori, we embedded the estimation method into a multilevel procedure, which also included a stability analysis, based on singular value decomposition of the sensitivity matrix. The optimization procedure was made more stable by imposing monotonicity constraints on the hydraulic functions. Tests with synthetic and measured data from column outflow experiments show the validity and robustness of the method.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY APPLICATION RESULTS DISCUSSION CONCLUSIONS REFERENCES

 
KNOWLEDGE of the hydraulic properties of unsaturated soils (i.e., the retention function relating the water content with the pressure head , and the hydraulic conductivity function K as a function of or ) is essential for most or all studies involving water flow and solute transport in the vadose zone. Most traditional methods to determine these properties require relatively restrictive initial and boundary conditions, and thus can be time-consuming, laborious, and expensive. Moreover, these methods generally aim at identifying either the retention or conductivity properties. Because of major advances in computational techniques and computer power, estimation of the unsaturated soil hydraulic properties by means of inverse modeling has now become an attractive alternative to traditional methods (Hopmans et al., 2002).

Popular laboratory approaches for the inverse estimation of the soil hydraulic properties have been one-step or multistep outflow methods (Kool et al., 1987; Durner et al., 1999a). Whereas the general feasibility of outflow experiments for identifying hydraulic functions has been demonstrated in several studies (e.g., Hopmans and Simunek, 1999), the applicability and success of this method is often still limited by the ill-posedness of the problem being considered. Ill-posedness in general is affected by nonoptimal experimental design, poor measurement quality, errors in the assumed process model, or unsuitable parameterizations of the hydraulic functions (Hornung, 1983, 1990; Kool et al., 1987; Yeh and Simunek, 2002; Vrugt et al., 2003).

Parameterizations of the soil hydraulic functions serve to constrain the inverse problem. The number of degrees of freedom determines the flexibility of the parameterized hydraulic functions. Most commonly used are the functions of van Genuchten–Mualem (van Genuchten, 1980) and Brooks–Corey–Burdine (Brooks and Corey, 1966), which contain a limited number of adjustable coefficients. Since the individual coefficients are related to properties of the pore-size distribution, they are generally referred to as "parameters" (e.g., parameter n in the van Genuchten equation refers to the width of the pore-size distribution). The low number of parameters may be achieved in part by coupling the retention function with the conductivity function through common parameters (again, parameter n can be considered as an example). Furthermore, in these low-parameterized functions, the degrees of freedom have a global character in that parameter changes will affect the hydraulic functions over the entire pressure head domain (–,0), even if only one degree of freedom is modified.

The use of low-parameterized functions can cause problems in cases where the selected hydraulic functions are not flexible enough to represent the actual shape of the hydraulic properties, such as for structured soils or for soils having pore-size distributions that do not have an approximately lognormal shape. In such cases, systematic errors in the fit of the retention function may propagate into the estimated conductivity function in a disproportionate manner (Durner et al., 1999a). To overcome this problem, piecewise continuous and composite (bimodal) functions have been introduced (Othmer et al., 1991; Durner, 1992; Wilson et al., 1992; Ross and Smettem, 1993; Mohanty et al., 1997). Unfortunately, when such functions are used in models that couple the conductivity function with the retention curve (e.g., Burdine, 1953; Mualem, 1976), the approach may lead to instability in the conductivity prediction near saturation (Durner, 1994). Furthermore, the identifiability of individual parameters in these composite functions is generally much worse than for the low-parameterized functions. Consequently, their use in inverse procedures for estimating the hydraulic properties has been limited to feasibility studies so far (e.g., Zurmühl and Durner, 1998). One major difference between the use of unimodal and bimodal descriptions is that some of the parameters in the bimodal approach no longer have an effect on the global shape of the hydraulic functions.

This study focuses on local rather than global parameterization (i.e., perturbations of single parameters affect just segments of the hydraulic functions). We investigate the use of free-form parameterizations of the hydraulic properties using piecewise polynomial functions in which each degree of freedom influences the hydraulic functions only locally within a small range of water potentials. The term free-form is used to indicate that no a priori assumption is made about the specific shape of the hydraulic functions and the number of parameters (degrees of freedoms). Splines and other interpolation methods have been used previously to describe the soil hydraulic properties (Erh, 1972; Kastanek and Nielsen, 2001; Prunty and Casey, 2002). However, to our knowledge, our study is the first where this type of description is used for inverse modeling of transient flow experiments.

The free-form parameterization inverse approach will be used here to analyze several transient outflow experiments in which an initially saturated soil column is drained by decreasing the pressure head at the lower outflow boundary. Water flow across the upper boundary is forced to be zero. The suitability of outflow experiments to identify hydraulic parameters has been investigated since the early 1980s (e.g., see review by Kool et al., 1987). The current standard in applying the technique was described by Hopmans et al. (2002), and Durner et al. (1999a) published a review on problems and current limitations of the methodology.

	THEORY

TOP ABSTRACT INTRODUCTION THEORY APPLICATION RESULTS DISCUSSION CONCLUSIONS REFERENCES

 
Water flow through the soil column is described with the Richards equation in one spatial dimension

[1]

where t (T) is time, z (L) is depth, positive downwards, (L) is the pressure head, () represents the water retention curve, and K() (L T^–1) is the hydraulic conductivity function. Inherent to this description is the assumption that water flow takes place as isothermal one-phase flow in a rigid porous medium. Equation [1] is subject to the following initial and boundary conditions:

[2]

where L denotes the length of the soil column and T denotes the end time.

Free-Form Parameterizations
DuChateau (1997) showed that values of the pressure head during an outflow experiment satisfy the inequality

[3]

if the function _LB(t) is continuously differentiable with ^'_LB < 0 for 0 < t < T and _LB(0) = L. Consequently, an outflow experiment does not provide information for estimating the hydraulic functions outside of the pressure interval [, 0], where = _LB(T) – L. Furthermore, the hydraulic functions are assumed to be constant when the soil is saturated ( > 0), such that it is sufficient to describe the hydraulic functions on the interval [, 0]. For physical reasons the hydraulic functions are additionally assumed to be nonnegative and monotonously increasing. The problem of monotonicity will be discussed in more detail below.

For the interval [, 0] we use a free-form cascaded parameterization. For this purpose we partition the interval [, 0] into:

[4]

For simplicity we will use equidistant partitioning with h = _j – _j–1 for j = 2, ..., r. We can now proceed in several slightly different ways. Below we describe two possible free-form parameterizations: quadratic B-spline parameterization and piecewise cubic Hermite interpolation.

Quadratic B-Spline Parameterization
Quadratic B-splines consistent with the partitioning of Eq. [4] are defined as (e.g., de Boor, 1978; Prenter, 1975)

[5]

where ₁, ₂ and _r, _r+1 are defined within the boundaries at and 0, respectively. Outside the interval [, 0] the functions are kept constant. These quadratic B-splines serve as shape functions such that the hydraulic functions are defined by

[6]

where = r + 1 is the number of degrees of freedom for each function () and K() assuming partitioning (Eq. [4]) with r nodes, and p_j (j = 1, ..., 2) denotes the degrees of freedom of the parameterization (which must be determined). Because of the local support of the shape functions, each degree of freedom p_j or p_j+ influences or K, respectively, only locally on the subinterval (_j–2, _j+1). It is possible to choose several alternative piecewise polynomial shape functions instead of quadratic B-splines. However, linear B-splines are unsuitable as shape functions because they lead to a nondifferentiable function () with a corresponding discontinuous water capacity function.

Since parameterization (Eq. [6]) does not automatically ensure monotonicity in the soil hydraulic functions, we require the p_i to satisfy additional constraints of the form

[7]

where = ₁ < ... < = 0 is a second partitioning of the interval [, 0] with a finer subdivision. For the parameterization using quadratic B-splines we need at least one additional subdivision of partitioning (Eq. [4]). Even then, however, monotonicity cannot be guaranteed completely. In some cases, numerical tests showed minor violations of the monotonicity requirement of the hydraulic functions. Having nonmonotonous hydraulic functions is physically unrealistic and could hamper convergence of numerical solutions of the Richards equation. For these reasons we also considered piecewise cubic Hermite interpolation since in this method monotonicity in the degrees of freedom implies monotonicity in the parameterized hydraulic functions.

Piecewise Cubic Hermite Interpolation
We again use Eq. [4] to partition the interval [, 0]. The degrees of freedom are now defined as the values of the hydraulic functions at the nodes of the partitioned interval; that is,

[8]

On each subinterval (_j, _j+1) the functions () and K() are determined by cubic Hermite interpolation between p_j and p_j+1 or p_r+j and p_r+j+1, respectively (Prenter, 1975). Hermite interpolation requires derivatives of the and K functions at each node in addition to the values of and K. From Fritsch and Carlson (1980) we know that the derivatives can be computed from the values p_j, j = l, ..., 2r, such that the resulting interpolates are monotonously increasing if the degrees of freedom satisfy

[9]

where r (=) is the number of degrees of freedom for each unknown function.

Solution of the Forward Problem
The Richards equation was numerically solved using the hybrid mixed finite element method of Schneid and Knabner (1997). This method uses the Darcian flow rate

[10]

as a second variable which, like the pressure head, is approximated explicitly. This mixed discretization leads to local mass conservation of water. Discretization in time was accomplished using the backward Euler method, while the cumulative water outflow rate was computed from the volumetric flow rate using a trapezoidal quadrature rule.

Treatment of the Inverse Problem
In the following, the relative permeability and capillary pressure functions are both partitioned into the same number of intervals (alternative approaches are possible). Consequently, let p represent the vector of the 2 adjustable degrees of freedom of the invoked parameterization. The inverse problem is solved by means of a least-squares approach in which the objective function

[11]

is minimized subject to the monotonicity constraints given by Eq. [7] or [9]. In Eq. [11], Y_m contains the measured values, Y_s are the numerically calculated values as a function of p, N is the number of measurements, and W is the diagonal weighting matrix.

According to DuChateau (1997), measurements of the pressure head at the upper boundary of the soil column (z = 0) and the cumulative outflow rate at the lower boundary (z = L) permit a unique solution of the inverse problem. We hence use these two data types in our objective function. Since only the derivative of the water content appears in Eq. [1], the saturated water content, _s, and the residual water content, _r, cannot be determined simultaneously by inverse modeling of the outflow data. For this reason we fixed the value of (0) = _s in the minimization. To improve estimation of the hydraulic conductivity near saturation, the saturated conductivity K(0) = K_s was also fixed. Minimization of the objective function O(p), with the monotonicity restrictions, was accomplished using a sequential quadratic programming method. Details of this procedure are given by Schittkowski (1988).

Sensitivity Analysis and Multilevel Procedure
Inverse modeling results using free-form parameterization depend very much on the number of degrees of freedom in each of the hydraulic functions. The number of degrees of freedom determines the flexibility of the parameterized hydraulic functions. This flexibility and the effects of measurement errors in the data both increase with increasing . Hence, an appropriate value of cannot be established a priori but must be determined during the solution of the inverse problem. The aim is to achieve a good balance between having flexibility in the hydraulic functions and limiting the effect of errors in the data.

We performed a linear sensitivity analysis similar to that used by Chavent et al. (1994). Assume that the minimum of Eq. [11] for data set Y_m occurs at p, and the minimum for the data Y_m + Y_m at p + p, where Y_m denotes the measurement errors and p the associated parameter deviations. A first-order expansion leads to the approximate problem of minimizing

[12]

The minimum of Eq. [12] in the least squares sense occurs at

[13]

This equation shows that the measurement errors and parameter changes are coupled via the sensitivity matrix [Y_s(p)]/p containing the derivatives of the simulated data with respect to the parameters. We assume here that 2 < N and that the sensitivity matrix has a full rank 2.

Singular value decomposition (Golub and van Loan, 1996) of the sensitivity matrix leads to a relationship between the identification error and the measurement error as follows

[14]

in which the right-hand side is the solution of Eq. [12] and [13], and where Y_m and the sensitivity matrix are appropriately modified to incorporate the weights in W. Furthermore, _j (j = l, ..., 2) are the singular values of the sensitivity matrix, whereas {u_j} and {v_j} form orthogonal bases of the range of the sensitivity matrix and the range of the transposed sensitivity matrix, respectively. Since the singular values determine the influence of data errors on the identification errors, we define two characteristic numbers, the spectral condition number of the sensitivity matrix.

[15]

and the maximum error amplification

[16]

which may be used to quantify the ill-posedness of the inverse problem. The sensitivity matrix may also be used to compute the gradient of the objective function, which is required for the invoked sequential quadratic programming method. Efficient methods for calculating the sensitivity matrix, which are based on applying the chain rule and do not use finite difference approximations, are described by Bitterlich and Knabner (2002b).

The minimization approach for the objective function and the sensitivity analysis were embedded into a multilevel algorithm. Minimization starts with the least possible number of degrees of freedom. The number of degrees of freedom is subsequently increased step by step until a certain stopping criterion is satisfied, or the maximum number of degrees of freedom is achieved. At each level the initial parameter values for the minimization are set by interpolating the minimization result of the previous level to the new parameterization. By making use of the sensitivity analysis, a stopping criterion is formulated that terminates the multilevel identification process if the characteristic numbers exceed given tolerance parameters:

[17]

where indicates the presumed measurement error of the adopted experimental setup. Suitable tolerance parameters can be found from numerical experiments by comparing the development of singular values with the fitting errors. Most problems will show a slow increase in the characteristic numbers. When the characteristic numbers increase dramatically, the multilevel estimation process should be terminated. We refer the reader to Bitterlich and Knabner (2001)(2002a, 2002b, 2003) and Igler and Knabner (2000) for additional details about the identification method.

	APPLICATION

TOP ABSTRACT INTRODUCTION THEORY APPLICATION RESULTS DISCUSSION CONCLUSIONS REFERENCES

 
The free-form parameterization approach was used to estimate the soil hydraulic properties from a set of synthetic outflow data, with and without adding noise, and from several measured data sets. Below we show results for the synthetic outflow data with added noise, and for data from two different outflow experiments (continuous and multistep outflow).

The synthetic data were generated by simulating an outflow experiment for a soil column of Guelph loam (drying curve) with a length of 15 cm. The hydraulic functions used for this simulation were given by the van Genuchten–Mualem model with _s = 0.52, _r = 0.218, = 0.0115 cm^–1, n = 2.03, m = l – l/n, K_s = 1.3167 cm h^–1, and = 0.5 (van Genuchten, 1980). The pressure applied at the lower boundary was decreased continuously from the initial equilibrium condition with _LB(0) = +15 cm to _LB(40) = –200 cm after 40 h. The simulated cumulative outflow and the pressure head at the upper boundary were used in the optimization, first without noise and then with 5% Gaussian noise added.

The measured data were obtained by means of both a continuous and a multistep outflow experiment on an undisturbed core taken from a field site near Bayreuth, Germany. The sample (diameter 10 cm, vertical length 15 cm) was taken at the 10- to 25-cm depth from the Ap (0–35 cm) horizon of a Stagnic Cambisol. The soil is a very dense loamy sand (78% sand, 18% silt, 4% clay, 1% organic matter content) with a volumetric saturated water content of 0.31 and a saturated hydraulic conductivity of 100 cm d^–1. The samples were mounted on a sintered glass plate (height 0.7 cm, air entry = –150 cm, K_s = 120 cm d^–1) and slowly saturated from the bottom with de-aired water. The saturated conductivity was measured by establishing a constant-gradient flux from the bottom to the top of the column. The outflow experiments started by applying a continuously or stepwise increasing suction to the lower boundary. A fully automated computer controlled device was used to control the experiment and monitor outflow at the bottom, while a tensiometer at a depth of 3.5 cm from the top of the soil core provided pressure head data. A detailed description of the laboratory device is given by Zurmühl (1998).

Two types of outflow experiments were performed sequentially on the same sample. In the "Continuous" experiment, the pressure at the bottom of the column was decreased continuously (and linearly in time) from +15.2 cm at t = 0 to –60 cm at t = 127 h. In the "Multistep" experiment the pressure at the bottom was decreased stepwise from +15.2 cm to –60 cm at t = 100 h. The sample was resaturated between the two experiments as described above.

As mentioned, values of _s and K_s were fixed during the optimization. Since the objective function consisted of two types of data with different units and magnitudes, we weighted the data with the inverse of the squares of their mean values (Kool et al., 1987). Also, the multilevel algorithm requires the parameterization to be refined by equidistant bisectioning of the pressure interval [, 0]. This means that we used a sequence of equidistant partitions with r = 2, 3, 5, 9, 17, and 33 nodes in the multilevel algorithm.

To compare the free-form optimization method with a more classic approach, additional optimizations were performed with the standard van Genuchten–Mualem functions. In this approach the van Genuchten coefficients , n, and _r were optimized using the Levenberg–Marquardt method. A modified version of the HYDRUS-1D code (Simunek et al., 1998) was used to solve the direct problem. Both _s and K_s were again fixed at their measured values as in the free-form parameterization method, the tortuosity parameter was assumed to be 0.5, and m was constrained by the relation m = 1 – l/n. The objective function was taken to be the same as for the free-form method.

	RESULTS

TOP ABSTRACT INTRODUCTION THEORY APPLICATION RESULTS DISCUSSION CONCLUSIONS REFERENCES

 
Synthetic Data
Figure 1 shows results of the inverse analysis using a quadratic B-spline parameterization for fitting the noisy synthetic data, with r = 2, 5, and 33 nodes. We note that for quadratic B-spline parameterization each of the hydraulic functions has one degree of freedom more than the number of nodes. For piecewise cubic Hermite interpolation, the number of nodes is the same as the number of degrees of freedom. The "true" hydraulic functions used to generate the synthetic data are also shown for comparison. Deviations between the synthetic data and the simulations, and between the true and optimized functions, are shown below the figures (amplified by a factor of two for better visual inspection).

View larger version (35K): [in this window] [in a new window]   Fig. 1. Quadratic B-spline parameterization results. Estimated (solid line) and original (dashed line) hydraulic functions and corresponding simulated (black line) and measured (points) data with 5% Gaussian noise. From left to right: 2, 5, and 33 nodes. Corresponding residuals are magnified by a factor of two (area plots).

View larger version (14K): [in this window] [in a new window]   Fig. 2. Evaluation of the objective function in the multilevel algorithm for the considered examples and values of the condition number of the sensitivity matrix for different parameterization levels.

For the case without added noise (results not further shown here), we obtained almost identical results for r = 2, 3, and 5, with nearly perfect fits of the data, while values of the objective function were essentially zero. Optimization of the data without noise produced about the same low fitting errors in the hydraulic conductivity as optimization of the data with 5% noise. Not fixing K_s produced much larger errors in both cases.

Cubic Hermitian parameterization was found to produce very similar results as for the B-splines (not further shown here) for both the noise-free and noisy data. Initially, for low r values (r < 5), values of the objective function and the fitting error were greater than those for quadratic parameterization because of less flexibility (one less degree of freedom). For the higher parameterization levels, we obtained almost identical results for the objective functions and hydraulic properties. Hence, except for the monotonicity issues discussed above, no reasons exist for preferring one method over the other. Rather, one may want to decide a posteriori whether one of the two methods is superior for a specific case.

Fit to Experimental Data
The free-form parameterization approach appears particularly suited for the hydraulic functions of porous media having bimodal pore-size distributions, as are commonly found for fine-textured and aggregated soils, or fractured rock (Othmer et al., 1991; Coppola, 2000; Tuller and Or, 2002; Peters and Klavetter, 1988). While the Bayreuth soil used in this study had a relatively unimodal pore-size distribution, attempts to describe its flow behavior with the classic van Genuchten–Mualem model were not very successful (Durner et al., 1999b). Neither did the results improve by increasing the number of degrees of freedom in the classic functions (e.g., by allowing the saturated conductivity, K_s, and tortuosity coefficient, , to be adjustable parameters). Moreover, the classic functions obtained for this soil using inverse modeling differed for different experimental procedures (continuous vs. multistep outflow) and for different methods of weighting the data in the objective function (outflow vs. outflow and pressure head data). These difficulties led us to believe that for this problem either the measurements were incorrect, the Richards equation was invalid, or the invoked analytical van Genuchten–Mualem functions for the hydraulic properties were too restrictive.

We now discuss the two consecutive outflow experiments (continuous and multistep) performed on the same sample. We first discuss the continuous outflow experiment and compare results obtained with the free-form model with those using the classic van Genuchten–Mualem approach. Figure 3 shows the measured cumulative outflow (top row) and pressure head (second row) data. Also shown are the simulated data using the free-form model with r = 2 (left column), 9 (central column), and 33 (right column), as well as amplifications of the fitting errors (black area plots). The third and fourth rows of Fig. 3 depict the corresponding hydraulic functions. We selected r = 9 for the central row based on an analysis of the error functional and the condition number for increasing levels of parameterization (Fig. 4) . The fitting error decreased significantly until r = 9, whereas the condition number increased continually with increasing parameterization level. From this we conclude that r = 9 provides an appropriate degree of flexibility in fitting the data. To compare this fit with the optimal fit that can be obtained by using the classic van Genuchten parameterization, we include the latter one in the central column of Fig. 3, together with an enlargement of the deviations (gray area).

View larger version (35K): [in this window] [in a new window]   Fig. 3. Continuous outflow results for the Bayreuth soil. From top to bottom: Observed (dots, every second data point shown) and fitted cumulative outflow data, observed and fitted pressure head data at the 3.5-cm depth, and the corresponding retention and conductivity functions for quadratic B-spline parameterization. From left to right: 2, 9, and 33 nodes (solid black lines). The central plots also show results obtained by fitting the van Genuchten–Mualem model to the data (dashed lines). Residuals are magnified by a factor of two for cumulative outflow and by a factor of four for the pressure head data, respectively (area plots).

View larger version (15K): [in this window] [in a new window]   Fig. 4. Values of the objective function in the multilevel algorithm for the Bayreuth continuous outflow experiment and values of the condition number of the sensitivity matrix for different parameterization levels.

Inspection of the hydraulic functions (Fig. 3, central column) shows that the estimated free-form retention function with r = 9 and the van Genuchten–Mualem retention function are almost identical. This suggests that the differences in simulated flow behavior were caused primarily or only by differences in the estimated hydraulic conductivity functions.

Figures 3 and 4 showed the results for the continuous outflow experiment. Similar results for the multistep outflow experiment are given in Fig. 5 and 6 . Analysis of the error functional and condition number (Fig. 6) shows that r = 9 again provides an appropriate degree of flexibility since the condition number increased substantially whereas the fitting error had reached its minimum. Similarly as for the continuous outflow experiment, the classic van Genuchten–Mualem model fitted the first part of the experimental data very well (Fig. 5), whereas the second part again showed systematic deviations. The free-form model nearly perfectly reproduced both the outflow and the pressure head data for the entire experiment, except for some dynamic effects near the air entry pressure (to be discussed below). Even the amplified fitting errors did not show significant differences between the measurements and the simulations.

View larger version (35K): [in this window] [in a new window]   Fig. 5. Multistep outflow results for the Bayreuth soil. From top to bottom: Observed (red dots) and fitted cumulative outflow data, observed and fitted pressure head data at the 3.5-cm depth, and the corresponding retention and conductivity functions. From left to right: 2, 9, and 33 nodes (solid black lines). The central plot also shows results obtained by fitting the standard van Genuchten–Mualem model to the data (dashed blue lines). All residuals shown are magnified by a factor of two (area plots).

View larger version (16K): [in this window] [in a new window]   Fig. 6. Values of the objective function in the multilevel algorithm for the Bayreuth multistep outflow experiment and values of the condition number of the sensitivity matrix for different parameterization levels.

Comparison of Fig. 3 and 5 shows that the free-form retention functions identified from the continuous and multistep outflow experiments are very similar. Some differences are apparent between the conductivity functions, especially near saturation where the sensitivity is too low. We conclude from the almost perfect fit of the observed outflow and pressure head data, and from the similarity of the resulting hydraulic functions for the two outflow experiments, that the outflow processes were described well with the governing Richards equation. The only notable deviation during multistep outflow occurred at water contents close to the air entry pressure. Probably, the air-phase permeability in this pressure range was not high enough to allow for an unrestricted flow of air. This dynamic equilibrium phenomenon is often observed during multistep outflow experiments (e.g., Wildenschild et al., 1999) and may require application of a two-phase (air, water) flow model (Luckner et al., 1989; Schultze et al., 1999). Since this effect is generally restricted to only a relatively narrow pressure range, and its magnitude is often small, it should be negligible in most practical applications.

	DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY APPLICATION RESULTS DISCUSSION CONCLUSIONS REFERENCES

 
Probably the most important finding of our study is the remarkable robustness of the free-form hydraulic function parameter estimation process. Convergence toward the final hydraulic property estimates was not affected by the initial values for the free-form parameterization. Also, using a seemingly excessive number of degrees of freedom did not negatively affect the results, except possibly for some local deviations in the retention function because of overfitting. Several previous studies have raised questions of how many coefficients could be reasonably well optimized simultaneously from a flow experiment (e.g., Mous, 1993; Zurmühl and Durner, 1998). We believe that the data should decide the optimum number of degrees of freedom. When the data are of high quality and the experimental setup is appropriate with respect to the desired information content, then simultaneous optimization of 20 or more spline coefficients should not be a problem, as was shown in this study. We found that the number of degrees of freedom in the free-form approach is actually of secondary importance since in the overparameterized case small perturbations will only reduce the smoothness of the hydraulic functions, but not negatively affect the overall results. While the uncertainty as such of the individual coefficients of the free-form parameterization will increase with increasing number of nodes, these coefficients will not be used in the same way as individual soil hydraulic parameters (e.g., such as those used in solute transport models such as dispersivity).

The proposed method also avoids potential problems related to the shape of the hydraulic conductivity function near saturation. Inverse modeling with low-parameterized global hydraulic functions is unable to address the frequently observed inconsistency between the optimized and measured saturated hydraulic conductivity, especially for structured media (Mohanty et al., 1997; Schaap and Leij, 2000). This has also lead to the recommendation not to use K_s as a matching factor in hydraulic conductivity prediction models, but to prefer some measured value away from saturation (van Genuchten and Nielsen, 1985; Hopmans et al., 2002). Free-form estimation of the conductivity function requires fixing the conductivity at saturation because of the extremely low sensitivity of the function on outflow information near saturation. However, the resulting interpolation between the estimated unsaturated conductivity and the measured value at saturation will always be consistent with respect to the observed flow process.

Potential drawbacks of the proposed free-form parameterization method are (i) difficulties to implement free-form descriptions of the hydraulic properties in numerical models for variably saturated flow, (ii) the empirical (nonphysical) nature of the free-form parameterization coefficients, (iii) problems of extrapolating the resulting hydraulic functions beyond the moisture range covered by the experiment, and (iv) increased risk of obtaining hydraulic functions that are physically meaningless when the experimental data are incomplete or erroneous.

The first drawback should be relatively unimportant since many numerical models now read hydraulic properties from tables and/or operate internally with interpolations between tabulated values (e.g., Simunek et al., 1998). Since Hermitian and B-spline interpolations lead to smooth functions with monotonous and continuous derivatives for both the retention and conductivity function, their use in simulation models should not pose undue problems.

The second drawback (the empirical nature of the free-form parameters) could actually be an advantage in that it immediately reflects the nonphysical nature of most or all currently available analytical hydraulic property descriptions. Assigning physical meaning to soil hydraulic parameters (e.g., the tortuosity coefficient and other parameters in the van Genuchten–Mualem functions) can be misleading (Schaap and Leij, 2000) since this ignores the fact that the functions are essentially empirical expressions. This is further reflected by the debate on the physical meaning of the residual water content (e.g., Nimmo, 1991), and the saturated water content and saturated hydraulic conductivity (van Genuchten and Nielsen, 1985). This is even more the case for descriptions of the hydraulic parameters of soils having bimodal pore-size distributions. Such descriptions were initially presented within a purely empirical context to improve the ability to fit observed data (Durner, 1994), but were later redefined in terms of parameters reflecting the statistical properties of textural and structural pore systems (Coppola, 2000; Miyamoto et al., 2003).

A seemingly more serious problem of the proposed free-form method is the extrapolation of the functional relationships to moisture conditions outside the range of calibration. Since the free-form parameterization assumes local estimation, any extrapolation must be done by users with or without additional information of soil texture and related data. While in some cases such extrapolation may be justified (e.g., with certain pore-size distribution models and/or using independent estimates of the residual water content using pedotransfer function [Schaap et al., 1998]), measurements in the wet range generally provide little information about the shape of the hydraulic functions in the dry range. Overall, the free-form functions reflect the real information content that can be gained from the calibration experiment.

The fourth potential drawback addresses the independent parameterization of the retention and the conductivity functions. Decoupling enhances the risk of obtaining physically unrealistic results (e.g., a constant water content across a certain range of pressure heads and a simultaneously decreasing hydraulic conductivity function), especially when the measurements are incomplete or erroneous. In the classic approaches, the conductivity function is often coupled with the retention function by means of a statistical pore-size distribution model. One major reason for sticking to this coupling in the inverse analysis is to limit the number of unknown coefficients, thus avoiding nonuniqueness problems. However, our study indicates that the optimal number of parameters will depend on both the design of the selected flow experiment and the quality of the measured data. Coupling the retention and conductivity functions using shared parameters is not always necessary and may actually reduce flexibility in the hydraulic description. Still, we recognize that coupling serves as a stabilizing feature in hydraulic parameter estimation analyses since it forces the constitutive functions into shapes that are physically consistent. It may also allow identification of the hydraulic functions from outflow data alone in cases where the decoupled case could fail. Since coupling the conductivity curve with free-form parameterized retention curve using the model of Mualem (1976), or any other predictive pore-size distribution model poses no inherent problems, we will address this approach in a future study.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION THEORY APPLICATION RESULTS DISCUSSION CONCLUSIONS REFERENCES

 
We presented a free-form parameterization approach for inverse estimation of the unsaturated soil hydraulic properties. The invoked parameterizations using quadratic B-splines and piecewise cubic Hermite interpolation avoid the often experienced problems of ill-posedness of inverse problems, provided an appropriate number of degrees of freedom is selected. The free-form estimation approach was embedded into a multilevel procedure that also involved a sensitivity analysis. An important feature of the multilevel approach is that the number of degrees of freedom determining the flexibility of the resulting hydraulic functions is not fixed a priori, but is selected as a function of the quality of the available data. The approach provides great flexibility in the hydraulic property description depending on soil type and the quality and amount of experimental information. Tests with synthetic data and real data showed the versatility and robustness of the method for a variety of experimental conditions.

	REFERENCES

TOP ABSTRACT INTRODUCTION THEORY APPLICATION RESULTS DISCUSSION CONCLUSIONS REFERENCES

 

• Bitterlich, S., and P. Knabner. 2001. Form-free and cascadic identification of material laws in (un)saturated fluid flow from column experiments. Z. Angew. Math. Mech. 81(Suppl. 1):13–16.
• Bitterlich, S., and P. Knabner. 2002a. Adaptive and form-free identification of nonlinearities in fluid flow from column experiments. p. 63–74. In Z. Chen and R.E. Ewing (ed.) Fluid flow and transport in porous media: Mathematical and numerical treatment. AMS, Providence, RI.
• Bitterlich, S., and P. Knabner. 2002b. An efficient method for solving an inverse problem for the Richards equation. J. Comput. Appl. Math. 147:153–173.
• Bitterlich, S., and P. Knabner. 2003. Experimental design for outflow experiments based on a multi-level identification method for material laws. Inverse Prob. 19:1011–1030.
• Brooks, R.H., and A.T. Corey. 1966. Properties of porous media affecting fluid flow. J. Irrig. Drain. Div. Proc. ASCE 92(IR2):61–88.
• Burdine, N.T. 1953. Relative permeability calculations from pore size distribution data. Pet. Trans. AIME 198:71–78.
• Chavent, G., J. Zhang, and C. Chardaire-Riviere. 1994. Estimation of mobilities and capillary pressure from centrifuge experiments. p. 265–272. In H.D. Bui et al. (ed.) Inverse problems in engineering mechanics. Balkema, Rotterdam, The Netherlands.
• Coppola, A. 2000. Unimodal and bimodal descriptions of hydraulic properties for aggregated soils. Soil Sci. Soc. Am. J. 64:1252–1262.[Abstract/Free Full Text]
• de Boor, C. 1978. A practical guide to splines. Springer, Berlin.
• DuChateau, P. 1997. An inverse problem for the hydraulic properties of porous media. SIAM J. Math. Anal. 28:611–632.
• Durner, W. 1992. Predicting the unsaturated hydraulic conductivity using multi-porosity water retention curves. p. 185–202. In M.Th. van Genuchten et al. (ed.) Proc. of the International Workshop, Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils. Univ. of California, Riverside.
• Durner, W. 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res. 30:211–233.
• Durner, W., E. Priesack, H.J. Vogel, and T. Zurmühl. 1999b. Determination of parameters for flexible hydraulic functions by inverse modeling. p. 817–829. In M.Th. van Genuchten et al. (ed.) Proc. Int. Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. Univ. of California, Riverside.
• Durner, W., B. Schultze, and T. Zurmühl. 1999a. State-of-the-art in inverse modeling of inflow/outflow Experiments. p. 661–681. In M.Th. van Genuchten et al. (ed.) Proc. Int. Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. Univ. of California, Riverside.
• Erh, K.T. 1972. Application of spline functions to soil science. Soil Sci. 114:333–338.
• Finsterle, S., and J. Najita. 1998. Robust estimation of hydrogeologic model parameters. Water Resour. Res. 34:2939–2947.
• Fritsch, F.N., and R.E. Carlson. 1980. Monotone piecewise cubic interpolation. SIAM J. Numer. Anal. 17:238–246.
• Golub, G.H., and C.F. van Loan. 1996. Matrix computations. 3rd ed. John Hopkins University Press, Baltimore.
• Hopmans, J., and J. Simunek. 1999. Review of inverse estimation of soil hydraulic properties. p. 643–659. In M.Th. van Genuchten et al. (ed.) Proc. Int. Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. Univ. of California, Riverside.
• Hopmans, J.W., J. Simunek, N. Romano, and W. Durner. 2002. Simultaneous determination of water transmission and retention properties. Inverse methods. p. 963–1008. In J.H. Dane and G.C. Topp (ed.) Methods of soil analysis. Part 4. SSSA Book Ser. 5. SSSA, Madison, WI.
• Hornung, U. 1983. Identification of nonlinear soil physical parameters from an input-output experiment. p. 227–237. In P. Deuflhard and E. Hairer (ed.) Finite Elements in Water Resour. 18. Springer, Berlin.
• Hornung, U. 1990. Parameter identification. p. 755–764. In D.G. DeCoursey (ed.) Proc. Int. Symp. on Water Quality Modeling of Agricultural Non-Point Sources. Part 2. 19–23 June 1988. Utah State Univ., Logan.
• Igler, B., and P. Knabner. 2000. Structural identification of nonlinear coefficient functions in transport processes through porous media. p. 157–178. In H.J. Bungartz et al. (ed.) Lectures on Applied Mathematics. Springer, Berlin.
• Kastanek, F.J., and D.R. Nielsen. 2001. Description of soil water characteristics using cubic spline interpolation. Soil Sci. Soc. Am. J. 65:279–283.[Abstract/Free Full Text]
• Kool, J.B., J.C. Parker, and M.Th. van Genuchten. 1987. Parameter estimation for unsaturated flow and transport models—A review. J. Hydrol. (Amsterdam) 91:255–293.
• Luckner, L., M.Th. van Genuchten, and D.R. Nielsen. 1989. A consistent set of parametric models for the two-phase flow of immiscible fluids in the subsurface. Water Resour. Res. 25:2187–2193.
• Miyamoto, T., T. Annaka, and J. Chikushi. 2003. Soil aggregate structure effects on dielectric permittivity of an Andisol measured by time domain reflectometry. Available at www.vadosezonejournal.org. Vadose Zone J. 2:90–97.[Abstract/Free Full Text]
• Mohanty, B.P., R.S. Bowman, J.M.H. Hendrickx, and M.Th. van Genuchten. 1997. New piecewise-continuous hydraulic functions for modeling preferential flow in an intermittent-flood-irrigated field. Water Resour. Res. 33:2049–2063.
• Mous, S.L.J. 1993. Identification of the movement of water in unsaturated soils: The problem of identifiability of the model. J. Hydrol. (Amsterdam) 143:153–167.
• Mualem, Y. 1976. A new model for predicting the hydraulic properties of unsaturated porous media. Water Resour. Res. 12:1248–1254.
• Nimmo, J.R. 1991. Comment on the treatment of residual water content in "A consistent set of parametric models for the two-phase flow of immiscible fluids in the subsurface" by L. Luckner et al. Water Resour. Res. 27:661–662.
• Othmer, H., B. Diekkrüger, and M. Kutilek. 1991. Bimodal porosity and unsaturated hydraulic conductivity. Soil Sci. 152:139–149.
• Peters, R.R., and E.A. Klavetter. 1988. A continuum model for water movement in an unsaturated fractured rock mass. Water Resour. Res. 24:416–430.
• Prenter, P.M. 1975. Splines and variational methods. John Wiley and Sons, New York.
• Prunty, L., and F.X.M. Casey. 2002. Soil water retention curve description using a flexible smooth function. Available at www.vadosezonejournal.org. Vadose Zone J. 1:179–185.[Abstract/Free Full Text]
• Ross, P.J., and K.J.R. Smettem. 1993. Describing soil hydraulic properties with sums of simple functions. Soil Sci. Soc. Am. J. 57:26–29.
• Schaap, M.G., and F.J. Leij. 2000. Improved prediction of unsaturated hydraulic conductivity with the Mualem–van Genuchten model. Soil Sci. Soc. Am. J. 64:843–851.[Abstract/Free Full Text]
• Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 1998. Neural network analysis for hierarchical prediction of soil hydraulic properties. Soil Sci. Soc. Am. J. 62:847–855.[Abstract/Free Full Text]
• Schittkowski, K. 1988. Solving constrained nonlinear least squares problems by a general purpose SQP-method. In K.-H. Hoffmann et al. (ed.) Trends in mathematical optimization. Birkhäuser, Basel, Switzerland.
• Schneid, E., and P. Knabner. 1997. numerical solution of unsteady saturated/unsaturated flow through porous media. p. 387–343. In M. Feistauer et al. (ed.) Numerical modelling in continuum mechanics, Part II.
• Schultze, B., O. Ippisch, B. Huwe, and W. Durner. 1999. Dynamic nonequilibrium in unsaturated water flow. p. 877–892. In M.Th. van Genuchten et al. (ed.) Proc. Int. Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. Univ. of California, Riverside.
• Simunek, J., and M. Sejna, and M.Th. van Genuchten, 1998. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat and multiple solutes in variably-saturated media. Version 2.0. IGWMC-TPS-70. Int. Groundwater Modeling Center, Colorado School of Mines, Golden, CO.
• Tuller, M., and D. Or. 2002. Unsaturated hydraulic conductivity of structured porous media: A review of liquid configuration–based models. Available at www.vadosezonejournal.org. Vadose Zone J. 1:14–37.[Abstract/Free Full Text]
• van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Abstract/Free Full Text]
• van Genuchten, M.Th., and D.R. Nielsen. 1985. On describing and predicting the hydraulic properties of unsaturated soils. Ann. Geophys. 3:614–628.
• Vrugt, J.A., W. Bouten, H.V. Gupta, and J.W. Hopmans. 2003. Toward improved identifiability of soil hydraulic parameters: On the selection of a suitable parametric model. Available at www.vadosezonejournal.org. Vadose Zone J. 2:98–113.[Abstract/Free Full Text]
• Wildenschild, D., J.W. Hopmans, and J. Simunek. 1999. Flow rate dependence of soil hydraulic characteristics. Soil Sci. Soc. Am. J. 65:35–48.
• Wilson, G.P., P. Jardine, and J. Gwo. 1992. Modeling the hydraulic properties of a multi-region soil. Soil Sci. Soc. Am. J. 56:1731–1737.[Abstract/Free Full Text]
• Yeh, T.-C., and J. Simunek. 2002. Stochastic fusion of information for characterizing and monitoring the vadose zone. Available at www.vadosezonejournal.org. Vadose Zone J. 2002:207–221.
• Zurmühl, T. 1998. Capability of convection dispersion transport models to predict transient water and solute movement in undisturbed soil columns. J. Contam. Hydrol. 30:99–126.
• Zurmühl, T., and W. Durner. 1998. Determination of parameters for bimodal hydraulic functions by inverse modeling. Soil Sci. Soc. Am. J. 62:874–880.[Abstract/Free Full Text]

This article has been cited by other articles:

				J. A. Vrugt, P. H. Stauffer, Th. Wohling, B. A. Robinson, and V. V. Vesselinov Inverse Modeling of Subsurface Flow and Transport Properties: A Review with New Developments Vadose Zone J., May 27, 2008; 7(2): 843 - 864. [Abstract] [Full Text] [PDF]

				G. Crescimanno and P. Garofalo Application and Evaluation of the SWAP Model for Simulating Water and Solute Transport in a Cracking Clay Soil Soil Sci. Soc. Am. J., October 27, 2005; 69(6): 1943 - 1954. [Abstract] [Full Text] [PDF]

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire