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 (3) Citing Articles via Google Scholar Google Scholar Articles by Kelleners, T. J. Articles by Skaggs, T. H. Search for Related Content PubMed Articles by Kelleners, T. J. Articles by Skaggs, T. H. GeoRef GeoRef Citation Agricola Articles by Kelleners, T. J. Articles by Skaggs, T. H. Related Collections Lysimeter/Rhizosphere Studies Inverse Procedures/Parameter Estimation Hydraulic Conductivity

ORIGINAL RESEARCH

Inverse Analysis of Upward Water Flow in a Groundwater Table Lysimeter

T. J. Kelleners^a,b,c,*, R. W. O. Soppe^b,d, J. E. Ayars^b, J. imnek^a,e and T. H. Skaggs^a

^a USDA-ARS, George E. Brown, Jr. Salinity Lab., 450 W. Big Springs Road, Riverside, CA 92507
^b USDA-ARS, Water Management Research Lab., 9611 S. Riverbend Ave., Parlier, CA 93648
^c Presently at Dep. of Plants, Soils, and Biometeorology, Utah State Univ., Logan, UT 84322
^d Alterra-ILRI, P.O. Box 47, 6700 AA Wageningen, The Netherlands
^e Presently at Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92521
^* Corresponding author (tkelleners{at}cc.usu.edu)

¹ The mention of trade or manufacturer names is made for information only and does not imply an endorsement, recommendation, or exclusion by the USDA-ARS.

Received 9 August 2004.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY HYPOTHETICAL SOIL MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 
The accuracy of numerical water flow models for the vadose zone depends on the estimation of the soil hydraulic properties. In this study, the hydraulic parameters for a silty clay soil in a large lysimeter were determined through inverse modeling of a fallow period with upward water flow from a shallow groundwater table. Parameter uniqueness was studied by simulating a hypothetical soil with known hydraulic properties under comparable conditions. Sensitivity analysis showed that the pressure head h(z,t), the volumetric water content (z,t), and the cumulative bottom flux Q(t) were least sensitive to the residual volumetric water content _r and the pore-connectivity parameter in the van Genuchten–Mualem (VGM) model. Parameter response surfaces showed that least squares fitting with (z,t) data is more likely to result in a unique hydraulic parameter set than least squares fitting with h(z,t) or Q(t) data. With only (z,t) in the objective function, the least squares minimization algorithm was capable of finding the correct soil hydraulic parameters, provided that _r and were fixed and that multiple initial parameter estimates were used. The protocol that was developed for the hypothetical soil was subsequently applied to the actual groundwater table lysimeter. The soil hydraulic parameters for the lysimeter for two (x,y) locations were determined using (z,t) data as measured by capacitance sensors. The variability in the optimized inverse of the air-entry value and the saturated hydraulic conductivity K_s in the VGM model was relatively high because of the high parameter correlation between these parameters. The optimized soil hydraulic properties can be used to study capillary rise from the groundwater table.

Abbreviations: CV, coefficient of variability • EC, electrical conductivity • LM, Levenberg–Marquardt • NRMSE, normalized root mean square error • VGM, van Genuchten–Mualem

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY HYPOTHETICAL SOIL MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 
ENVIRONMENTAL IMPACT and irrigation water management studies often rely on numerical models to predict water flow and solute transport in the vadose zone. The predictions of these models are sensitive to the input of the soil hydraulic properties. The soil hydraulic properties are usually described using empirical or semitheoretical water retention and hydraulic conductivity functions (e.g., Gardner, 1958; Brooks and Corey, 1964; Mualem, 1976; van Genuchten, 1980; Kosugi, 1999).

Traditionally, hydraulic properties have been measured in the laboratory using soil samples taken from the field. Applying the soil hydraulic properties derived from these small-scale experiments to the larger field scale often leads to disappointing results. Soil disturbance during sampling and spatial variability in the soil hydraulic properties play major roles in the failure to obtain suitable field-scale parameter values for the hydraulic functions (Kool et al., 1987).

Alternatively, soil hydraulic properties may be determined in situ. This avoids soil disturbance and provides an estimate of the hydraulic properties that integrates soil heterogeneity. Also, if the soil hydraulic functions are ultimately to be used to analyze field-scale processes, in situ estimation intuitively seems more appropriate than laboratory analysis.

Parameter optimization has emerged as an important technique for estimating soil hydraulic parameters. Both laboratory and field experiments can now be run and analyzed under a large variety of flow conditions. In the past, these experiments were restricted to certain well-defined conditions that allow the calculation of the hydraulic parameters by solving analytical solutions of the relevant flow equations. Today, the hydraulic parameters can be obtained from flow experiments by combining a numerical flow model with a parameter estimation code. Recent advances in the computational power of computers have made this method even more attractive.

Numerous studies have been published on the combined use of flow experiments and parameter optimization to estimate the soil hydraulic parameters. Reviews of these studies can be found in Kool et al. (1987), Hopmans and imnek (1997), and Hopmans et al. (2002). A recurring issue is the problem of parameter uniqueness. If the objective function to be minimized (usually the sum of squared differences between measured and calculated flow variables) does not have a clear global minimum, the inversion will be nonunique. The occurrence of local minima in the objective function may also cause nonuniqueness. The occurrence of uniqueness problems depends on the soil type, the boundary conditions, the type of data used in the objective function, and the parameter estimation algorithm.

Several investigators have tried to infer soil hydraulic properties by applying the inverse modeling technique to upward infiltration experiments. The experiments provide information on the wetting branch of the soil hydraulic properties without the need to account for the macropore flow that might occur during downward (gravity) infiltration. So far, several laboratory experiments on small soil samples have been conducted using either a flux condition (Hudson et al., 1996) or a pressure condition (imnek et al., 2001; Young et al., 2002) at the bottom boundary.

In this work, the hydraulic properties of a silty clay soil in a weighing lysimeter are determined by inverse modeling of a summer fallow period where capillary rise from the groundwater table replenishes the depleted root zone. The lysimeter is 4 m long, 2 m wide, and 3 m deep, and as such, constitutes an intermediate scale between the sample scale and the field scale. The Hydrus-1D model (imnek et al., 1998a), which includes a parameter estimation routine, is used to solve the flow problem. The objectives of the study were (i) to investigate the issue of parameter uniqueness for the upward infiltration problem by testing the inverse procedure on a hypothetical soil with known soil hydraulic properties and (ii) to determine the soil hydraulic properties of the silty clay soil in the lysimeter.

	THEORY

TOP ABSTRACT INTRODUCTION THEORY HYPOTHETICAL SOIL MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 
Governing Flow Equation
The water flow calculations are conducted with the Hydrus-1D model, which numerically solves the Richards equation:

[1]

where is the volumetric water content (L³ L^–3), t is the time (T), z is the vertical coordinate (L) (positive upward), h is the pressure head (L), and K is the hydraulic conductivity (L T^–1).

Initial and boundary conditions need to be specified to solve Eq. [1]. The initial condition is specified in terms of pressure head or water content:

[2]

[3]

The soil surface boundary is described by an atmospheric boundary condition that switches between a prescribed flux condition and a prescribed head condition, depending on the prevailing transient pressure head at the soil surface. Whenever the surface pressure head is in the range h_A < h(0,t) < h_S, the atmospheric boundary is a prescribed flux condition:

[4]

where R is the potential rate of infiltration or evaporation (L T^–1) under the prevailing atmospheric conditions, and h_A (L) and h_S (L) are minimum and maximum allowed values of the surface pressure head, respectively. If the surface pressure reaches h_A or h_S, the surface boundary switches to a prescribed pressure head condition: h(0,t) = h_A or h(0,t) = h_S. In this study, we use h_A = –10⁵ cm (hypothetical soil), h_A = –10⁶ cm (lysimeter soil), and h_S = 0 cm.

A pressure head boundary condition is specified at the bottom of the lysimeter:

[5]

where h₀ is the prescribed value of the pressure head (L) and L is the depth (L) of the lysimeter.

Soil Hydraulic Properties
The soil hydraulic properties are described with the VGM model (van Genuchten, 1980; Mualem, 1976):

[6]

[7]

[8]

where S_e is the effective saturation, _r is the residual volumetric water content, _s is the saturated volumetric water content, is the inverse of the air-entry value (L^–1), n is a pore-size distribution index, K_s is the saturated hydraulic conductivity (L T^–1), and is a pore-connectivity parameter.

Inverse Procedure
The objective function that is minimized during the parameter estimation process is (e.g., imnek et al., 1998a):

[9]

where b is the vector of parameters to be optimized (e.g., b = [_r, _s, , n, K_s, ]), q_k^* is the measured value for the kth measurement set at depth z_j and at time t_i, q_k is the corresponding predicted value for parameter vector b, m_q is the total number of measurement sets, n_qk is the number of depths z for the kth measurement set, and o_qjk is the number of observation times t for depth z and measurement set k.

The different measurement sets are weighted using weighting coefficients w_k:

[10]

where _k² is the variance of the data in the kth measurement set. The use of w_k minimizes unwanted differences in weighting among data types caused by differences in magnitude and the number of data points.

Minimization of the objective function is accomplished by using the Levenberg–Marquardt (LM) nonlinear minimization method (Marquardt, 1963; imnek et al., 1998a). The LM method uses the steepest descent method when the objective function is far from its minimum and switches to the inverse Hessian method as the minimum is approached (imnek and Hopmans, 2002). The LM method is a local gradient-type search algorithm, as opposed to global search algorithms that search the entire parameter space. Local search algorithms are generally sensitive to the initial parameter estimates. Consequently, different initial estimates need to be examined.

	HYPOTHETICAL SOIL

TOP ABSTRACT INTRODUCTION THEORY HYPOTHETICAL SOIL MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 
Model Setup
Inverse estimation procedures for the lysimeter were investigated by considering a hypothetical soil with known soil hydraulic properties. Upward flow from a shallow groundwater table was simulated for conditions closely resembling the summer fallow period in the lysimeter. A homogeneous silty clay soil was assumed with a constant groundwater table at 1 m below the soil surface. The vertical profile was divided into 200 elements, each with a thickness of 0.5 cm. The parameters describing the hydraulic properties of the soil were obtained from the Rosetta pedotransfer program of Schaap et al. (2001). Soil texture data (2.6% sand, 43.9% silt, and 53.5% clay) and dry bulk density data (1.38 g cm^–3) were used as input for Rosetta. These soil data were obtained during a field experiment at the location where the lysimeter soil was taken (see Kelleners et al., 2004). The data were averages for the B1 horizon of this soil (20–75 cm below soil surface). The Rosetta-predicted soil hydraulic parameter values were: _r = 0.101, _s = 0.492, = 0.015 cm^–1, n = 1.321, K_s = 3.47 cm d^–1, and = –1.055. Note that the physical meaning of as a tortuosity factor is lost by allowing < 0. With = –1.055 the VGM model should be interpreted as a semitheoretical fitting equation (Schaap and Leij, 2000).

A 100-d fallow period was simulated with an initially dry soil profile. The initial pressure head values were assumed to increase linearly over five depth intervals with h_i(z = 0 cm) = h_A = –100000 cm, h_i(–10) = –30000 cm, h_i(–50) = –15000 cm, h_i(–70) = –1000 cm, h_i(–90) = –30 cm, and h_i(–100) = 0 cm. This initial pressure head distribution was intended to represent a soil profile that is depleted by root water uptake from a preceding crop. The top 10 cm of the soil was considered further depleted by continued evaporation from the soil surface after irrigation was halted toward the end of the growing season. An atmospheric boundary condition was specified at the top boundary with a constant potential evaporation rate of 0.5 cm d^–1 and zero rainfall and irrigation. At the bottom boundary a groundwater table condition was specified (h₀ = 0 cm).

The bottom boundary during the inverse analysis was specified as a "time-variable boundary condition" in Hydrus, albeit with constant zero pressure head. This prevented negative pressure heads at the bottom boundary when the LM algorithm tested _s values that were higher than the true _s value of 0.492. This problem only occurred because _i was used to specify the initial condition during the inverse analysis (with _i, each Hydrus run starts with the conversion of _i into h_i using the current estimate of the soil hydraulic parameters). The problem is a consequence of the incompatibility between _i(z) and h₀(t) at t = 0 and z = –L. Specifying the bottom boundary as a time-variable boundary condition may still result in some artificial redistribution near the bottom during the first few time steps, but the effect of the redistribution on the calculated water flow is negligible.

Data Generation
Simulated pressure heads and volumetric water contents at the end of each day at depths of 5, 15, 25, 35, 45, 55, 65, 75, 85, and 95 cm below the soil surface were recorded for use in the inverse analysis. This resulted in 1000 data points for (z,t) and h(z,t) each. As an example, the resulting volumetric water contents as a function of time and depth are shown in Fig. 1 . Note that during the 100-d calculation period the upward moving wetting front passes all measurement depths, thus providing a maximum amount of information for the inverse analysis. In addition, the cumulative upward flux through the bottom of the soil profile (22.1 cm) was subdivided into 1-d intervals, resulting in 100 Q(t) data points.

View larger version (22K): [in this window] [in a new window]   Fig. 1. Simulated volumetric water content for the 100-d fallow period with hypothetical soil hydraulic properties. The numbers in the figure refer to the depth below the soil surface.

Sensitivity Analysis
A sensitivity analysis was performed to study the relationship between the measurement data and the model parameters for the upward flow problem. The higher the sensitivity of a model parameter to the data, the higher the chance that the parameter is identifiable during the inverse process. Sensitivity coefficients for the hypothetical soil were calculated from (imnek et al., 1998b):

[11]

where s_kl is the change in variable q_k corresponding to a 1% change in parameter b_l, e_l is the lth unit vector, and b = 0.01 b. The 0.01b_l term in Eq. [11] is included to allow comparison between different parameters, independent of their invoked unit or absolute value.

Sensitivities s_kl of the pressure head h(z,t), the volumetric water content (z,t), and the cumulative bottom flux Q(t) were calculated for all six soil hydraulic parameters. All three measurement sets showed a decrease in sensitivity according to n > > _s > K_s > > _r, with n being the most influential parameter and _r being the least influential parameter. As an example, Fig. 2 shows the sensitivity of the pressure head at the 5-, 15-, and 35-cm depths to the six hydraulic parameters. Maximum sensitivity is observed at times when the upward moving wetting front passes a certain depth and large changes in the pressure head occur. Sensitivities are highest for the shallowest depth; this is due primarily to the nonlinear nature of the retention curve (small changes in the water content of relatively dry soil result in large changes in the pressure head). By comparison, the differences are less pronounced for the sensitivity of the water contents at different depths (data not shown). Figure 2 also shows that the duration of elevated sensitivity levels is longest for the shallowest depth. This is attributed to the more diffuse moisture front at this depth, a consequence of the increased distance from the groundwater table.

View larger version (17K): [in this window] [in a new window]   Fig. 2. Sensitivity of the pressure head in the hypothetical soil to a 1% change in the soil hydraulic parameters in the van Genuchten–Mualem model at depths of 5, 15, and 35 cm below the soil surface.

imnek and van Genuchten, 1996

View larger version (15K): [in this window] [in a new window]   Fig. 3. Sensitivity of the cumulative upward bottom flux in the hypothetical soil to a 1% change in the soil hydraulic parameters in the van Genuchten–Mualem model.

We need to stress that two-dimensional response surfaces provide only cross sections of the full six-dimensional parameter space. As such, these two-dimensional surfaces do not provide full proof about the uniqueness of the inverse problem. Nevertheless, the response surfaces are useful to study the behavior of the objective function in the parameter space. The inverse parameter estimation technique is expected to be unsuccessful if response surfaces do not display a clearly defined global minimum in the two-dimensional parameter planes (imnek et al., 1998b).

First we studied parameter uniqueness with only (z,t), h(z,t), or Q(t) in the objective function. Contour plots of the response surfaces for –n, –K_s, and n–K_s are shown in Fig. 4 . Single minima are evident for (,n;) (Fig. 4a) and (n,K_s;) (Fig. 4g) but not for (,K_s;) (Fig. 4d). In contrast, identifiable minima are absent in all response surfaces for h(z,t) and Q(t). From Fig. 4 it appears that (z,t) data contain more useful information for uniquely identifying , n, and K_s than h(z,t) or Q(t) data. Similarly, response surfaces also showed that _r and _s can only be identified with (z,t) data in the objective function (comparison not shown). With solely h(z,t) and Q(t) data in the objective function, identifying water retention parameters is not possible because there is no information about the volumetric water content of the soil. Note that some of the local minima in Fig. 4 are artifacts resulting from the selected discretization of 50 by 50 points. The lower the resolution, the more artificial local minima will appear, especially near parameter combinations where changes only gradually, such as near the global minimum.

View larger version (67K): [in this window] [in a new window]   Fig. 4. Contour lines of the objective function for the parameter combinations –n, –Ks, and n–Ks using either (z,t), h(z,t), or Q(t) data in the objective function.

View larger version (55K): [in this window] [in a new window]   Fig. 5. Contour lines of the objective function for nine parameter combinations using (z,t) data in the objective function.

View larger version (25K): [in this window] [in a new window]   Fig. 6. Contour lines of the objective function for the parameter combinations –n, –Ks, and n–Ks using (z,t) and Q(t) data in the objective function.

Model grid and boundary conditions for the inverse analysis were the same as for the forward simulation. However, the initial conditions were specified in terms of water content _i(z) instead of pressure head h_i(z), consistent with h(z,t) data not being available for the lysimeter. The LM minimization method was repeated 50 times for each case, using 50 different initial estimates of the soil hydraulic parameters. Each initial parameter set was generated by random selection from predetermined parameter intervals ([0–0.13] for _r, [0.3–0.6] for _s, [0.0005–0.05 cm^–1] for , [1.05–2.0] for n, [0–25 cm d^–1] for K_s, and [–5–10] for ). These intervals were also used as bounds during the parameter optimizations. The choice of 50 realizations was a compromise between statistical considerations (more realizations give more reliable means) and the computational effort. We omitted LM minimizations that experienced numerical errors during the Hydrus model runs or that did not converge within 20 iterations.

The results of the parameter estimations for the hypothetical soil are summarized in Table 1 [(z,t) data], Table 2 [Q(t) data], and Table 3 [(z,t) and Q(t) data]. Each table shows four different cases. In the first case all six hydraulic parameters are optimized simultaneously. In the other three cases _r, , and both _r and are fixed to their true values while the remaining parameters are optimized. Fixing parameters simplifies the minimization process by reducing the degrees of freedom in the inverse process. We chose to fix _r and because of the relatively low sensitivity of (z,t) and Q(t) to these two parameters.

View this table: [in this window] [in a new window]   Table 1. Results of the hydraulic parameter estimation for the hypothetical soil with volumetric water content (z,t) data in the objective function.

View this table: [in this window] [in a new window]   Table 2. Results of the hydraulic parameter estimation for the hypothetical soil with cumulative bottom flux Q(t) data in the objective function.

View this table: [in this window] [in a new window]   Table 3. Results of the hydraulic parameter estimation for the hypothetical soil with volumetric water content (z,t) and cumulative bottom flux Q(t) data in the objective function.

[12]

where b_l is the optimized value of the lth soil hydraulic parameter, _l is the average optimized value of the lth soil hydraulic parameter, and N is the number of converged LM minimizations. The normalized root mean square error (NRMSE, %) in Tables 1, 2, and 3 is a measure of the variability of the optimized parameters around the true parameter value:

[13]

where b_l^* is the true value of the lth soil hydraulic parameter. Both CV and NRMSE allow comparisons among different parameters, irrespective of their invoked unit or absolute value. The success rate in the last column of Tables 1, 2, and 3 refers to the number of minimizations for which all optimized parameter values were within 5% of their true values (number before the slash) as a fraction of the total number of converged minimizations N (number behind the slash).

Table 1 shows that fixing _r and is essential if only (z,t) data are used in the objective function. Even then, only 17 of 35 converged minimizations arrive at the correct soil hydraulic parameters. Also, considerable variability remains in the optimized values of (CV = 31.0%) and K_s (CV = 103.0%) with _r and fixed. The uncertainty in the and K_s values is reflected in the high parameter correlation of between 0.95 and 0.98 in the converged minimizations (numbers not shown). Note that –1 stands for perfect negative correlation, 0 for no correlation, and +1 for perfect positive correlation. The difficulty of identifying and K_s separately was anticipated because the response surface (,K_s;) (Fig. 4d) lacked a clear minimum. Table 1 also shows that failure to fix during the inverse process will lead to severe errors in the optimized value of K_s (NRMSE of 463.1 and 477.5%).

None of the LM minimizations with solely Q(t) data in the objective function converged to the true soil hydraulic parameters (Table 2). The absence of a single minimum in the response surfaces for –n, –K_s, and n–K_s (Fig. 4c, 4f, and 4i) already suggested that Q(t) data alone were insufficient to identify uniquely the hydraulic parameters. It is interesting to note that the absence of water content data in the objective function did not result in large errors in the optimized _s values (NRMSE 2.3–15.2%). It appears that the use of water content data for the initial condition helped in narrowing the range of possible values for this parameter. The error in the optimized _r value was larger (NRMSE 45.7–59.5%) because the lowest _i value of 0.138 at z = 0 was higher than the true _r value of 0.101. Estimation of _r therefore required extrapolation beyond the measurement range. With pressure head instead of water content as the initial condition in the inverse analysis it would have been impossible to approximate even _s because only the difference _s – _r could have been optimized with only Q(t) data in the objective function.

Combining (z,t) and Q(t) data in the objective function for the upward flow problem (Table 3) reduces the deviations between the optimized and true values of K_s and for all cases (judging from the NRMSE values). Comparison between Tables 1 and 3 also shows that the SSQ is lowered by combining (z,t) and Q(t) in the objective function, signaling a better fit between the measured and simulated water contents. However, _r and still need to be fixed during the inverse process to obtain the true soil hydraulic parameter values (17 of 20 converged minimizations are successful). With _r and fixed, the variability in the optimized values of and K_s is now small (NRMSE < 3%), despite the fact that the parameter correlation between these parameters remains high (0.97; numbers not shown).

The above results agree with findings of imnek and van Genuchten (1996), who, after studying downward infiltration from a tension disc infiltrometer, concluded that cumulative infiltration alone will not provide a unique solution for the inverse problem. In subsequent studies these authors suggested augmenting the cumulative infiltration data with measured final water contents (imnek and van Genuchten, 1997) or with (z,t) data (imnek et al., 1999) to obtain unique solutions. imnek et al. (1999) also noted that augmenting cumulative flux data with h(z,t) data did not improve the results. It appears that the upward flux problem (this study) and the downward infiltration study (imnek and coworkers) pose similar challenges for inverse parameter estimation. This is not a complete surprise since both flow problems involve infiltration into a dry soil from a pressure head boundary condition.

Stability of the Inverse Solutions
The stability of the soil hydraulic parameter estimates is examined by altering the setup of the inverse process in four ways. The computational effort is limited to the case where only (z,t) data are used for fitting, and where _r and are being fixed (success rate of 17/35 in Table 1). This case is most relevant for the studied groundwater table lysimeter, with no useful Q(t) or h(z,t) data available. First, the adequacy of using 50 replicates is demonstrated by increasing the number of LM minimizations from 50 to 100. Second, the effect of both underestimating and overestimating the _r value during the inverse analysis is quantified by fixing _r to 0.075 and 0.125, respectively (true _r = 0.101). Third, we assessed the consequences of using = 0.5, a number suggested by Mualem (1976) for most natural soils, instead of = –1.055. Fourth, we tested the impact of measurement errors in the water content data by adding normally distributed random errors with zero mean and 0.0025 cm³ cm^–3 standard deviation to the (z,t) data.

The results of the stability analysis are summarized in Table 4. Increasing the number of LM minimizations from 50 to 100 did not result in significant changes in the soil hydraulic parameter estimates (compare Tables 1 and 4). The biggest changes occurred in the mean K_s value (increases from 4.64 to 4.98 cm d^–1) and in the NRMSE value for K_s (increases from 141.9 to 152.7%). These increases are insignificant for a parameter like K_s, which is often found to vary one or two orders of magnitude, even in homogeneous materials. This confirms that 50 replicates suffice to obtain reliable statistics for the inverse solutions.

View this table: [in this window] [in a new window]   Table 4. Stability of the hydraulic parameter estimates for the hypothetical soil with volumetric water content (z,t) data in the objective function.

The incorporation of random measurement errors in the water content data does not significantly change the outcome of the inverse solutions. The mean parameter values remain about the same, the CV and NRMSE values increase slightly, and the SSQ values increase to 1.25 x 10^–2. Furthermore, the success rate is still 5/37, despite the erroneous water contents. Note that a standard deviation for the water content measurement error of 0.0025 cm³ cm^–3 is representative of certain electromagnetic techniques like time domain reflectometry (e.g., Heimovaara and Bouten, 1990; Lambot et al., 2002). The use of other, less consistent sampling techniques might increase the standard deviation and thereby the uncertainty in the estimated soil hydraulic parameters. The occurrence of persistent instrument errors might increase the uncertainty even more.

Implications
The response surfaces already indicated that the upward flow problem might suffer from nonuniqueness problems, even when (z,t) and Q(t) data are combined. Different combinations of parameter values may lead to equally small values of the objective function. This is especially true for the and K_s parameters. The LM minimization method will not always find the exact global minimum under these circumstances. It is therefore not surprising that the inverse solutions yield a collection of parameter sets that perform about equally well.

We would like to stress that the failure to consistently recover the true parameters during the inverse solutions is not due to the LM algorithm itself. There is simply no distinct global minimum in the objective function. The use of a global search algorithm (e.g., Abbaspour et al., 1997; Lambot et al., 2002), which examines the complete parameter space, would probably not help under these circumstances. A global search algorithm would only help if the LM minimization was frequently ending up in a local minimum that is distinctly different from the global minimum. Inspection of the individual minimization results showed that this is seldom the case.

The nonuniqueness problem implies that we cannot expect to find a unique set of soil hydraulic parameters for the groundwater table lysimeter. We can only expect to find a collection of parameter sets that describe the data equally well. Analysis of the CV values will provide important information about the reliability of the calculated mean parameter values. The stability analysis indicated that 50 LM minimization runs will suffice for this purpose.

The failure to consistently find the true soil hydraulic parameters for the upward flow problem introduces uncertainty into the modeled soil water fluxes for the groundwater table lysimeter. Calculated state variables like water content and pressure head will be a function of the chosen parameter set, which may be nonunique. However, the use of nonunique parameter sets may still yield valuable information on spatial and temporal trends in the state variables and in the water balance components.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION THEORY HYPOTHETICAL SOIL MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 
Experimental Setup
The groundwater table lysimeter (2 m width, 4 m long, and 3 m deep) was installed in 1995 at the USDA-ARS facility in Parlier, CA. The top 1.7 m of the lysimeter consisted of an undisturbed soil monolith from the west side of the San Joaquin Valley, CA. The soil was a saline silty clay (fine, smectitic, thermic Sodic Haploxerert). The bottom 1.3 m of the lysimeter consisted of disturbed soil from the same location, hand packed to a dry bulk density of 1.3 to 1.35 g cm^–3 (Schneider et al., 1996). A truck scale measured changes in soil water content with time. A mariotte bottle was used to maintain the groundwater table at about 1.0 m below soil surface. Over the years, the lysimeter was planted with cotton (Gossypium hirsutum L.), safflower (Carthamus tinctorius L.), and alfalfa (Medicago sativa L.), and was irrigated with good quality water (electrical conductivity, EC = 0.4 dS m^–1) by subsurface drip, surface drip, and sprinkler systems. The EC of the groundwater in the lysimeter was about 14 dS m^–1.

All components of the lysimeter water balance were measured directly except evapotranspiration, which was calculated by difference on an hourly basis. The depth of the groundwater table was measured with an observation well in the center of the lysimeter using a pressure transducer. The distribution of soil water content with depth was monitored at two locations in the lysimeter with multisensor EnviroSCAN capacitance probes (Sentek Pty Ltd., Kent Town, South Australia).¹ The capacitance probes each held 11 sensors at 10-cm intervals starting at 5.5 cm (Location L1) and 7.0 cm (Location L2) below the soil surface. The frequency readings of the capacitance sensors were converted into volumetric water contents by the procedure of Kelleners et al. (2004). In this procedure, the frequency was related to the soil permittivity by an equation that described the electromagnetic properties of the sensor–plastic access tube–soil system. The permittivity was then related to the volumetric water content using the empirical model of Malicki et al. (1996).

Model Setup
A summer fallow period (10 July 2001–26 Sept. 2001) was selected for the inverse analysis. During this period the root zone, which was depleted of water by the preceding safflower crop, was gradually replenished by capillary rise from the groundwater table. Daily rainfall (practically none) and reference evapotranspiration for the simulation period were taken from a California Irrigation Management Information System weather station located about 500 m from the lysimeter. The reference evapotranspiration was converted to potential evaporation for bare soil using the dual crop coefficient procedure of Allen et al. (1998). For two 4-d periods (3–6 Aug. 2001 and 14–17 Sept. 2001), the lysimeter was covered with a plastic sheet during sprinkler irrigation of the surrounding (fallow) field. Potential evaporation during these times was set to zero. The depth of the groundwater table as measured in the observation well varied between 88 and 100 cm below the soil surface during the calculation period.

The safflower crop preceding the fallow period was drip irrigated and planted in rows. This resulted in a horizontally and vertically heterogeneous soil moisture pattern. The heterogeneous pattern persisted after the crop was harvested and the irrigation was halted. During the inverse analysis of the fallow period only one-dimensional vertical flow was assumed for a single (x,y) location. Thus, the water balance fluxes measured with the lysimeter could not be used in the inverse analysis because they represent an areal average for the complete lysimeter. In contrast, the capacitance probes, which represent single locations, could be used, provided that the horizontal soil water fluxes at these locations were negligible compared with the vertical (upward) fluxes. Alternatively, the lysimeter could be described with a complete three-dimensional model. This was not attempted because of the lack of spatial data (e.g., initial water contents for only two [x,y] locations) and because of the large computational requirements.

Thus, only water content data derived from the capacitance sensors were used in the inverse analysis. However, not all capacitance data were included in the objective function. The top two L1 sensors (at the 5.5- and 15.5-cm depths) and the top three L2 sensors (at depths of 7 to 27 cm) were excluded because the upward moving wetting front never reached these sensors. In fact, these sensors recorded a decrease in the water content during the 79-d calculation period as a result of continued evaporation to the atmosphere. Elimination of the "drying" sensors enabled a pure estimate of the soil hydraulic properties during wetting. These wetting soil hydraulic properties might be different from the drying soil hydraulic properties because of hysteresis (e.g., Dane and Wierenga, 1975). The problem of "drying" sensors did not occur for the hypothetical soil because of the extremely dry initial conditions in that case. The bottom L1 sensor (at the 105.5-cm depth) and the bottom L2 sensor (at the 107-cm depth) were excluded because they were always below the groundwater table and therefore added little information for the parameter estimation process. The water content data used for fitting were restricted to eight L1 sensors (25.5–95.5 cm depth) and seven L2 sensors (37–97 cm depth).

Only the top 100 cm of the lysimeter was modeled with Hydrus-1D. The soil profile was divided into 109 elements that varied in size from 0.1 cm at the soil surface to 1.0 cm at depth. An atmospheric boundary condition was used at the top, and a variable pressure head boundary was used at the bottom. The capacitance probe data from 9 July 2001 were used to specify the initial _i(z) condition. In the objective function, simulated water contents were averaged over 6-cm depth intervals to be consistent with the vertical zone of influence of the capacitance sensors (e.g., Paltineanu and Starr, 1997). The value of h_A was decreased from –10⁵ to –10⁶ cm after trial calculations showed that maintaining h_A = –10⁵ sometimes caused unrealistically high water contents in the topsoil.

Note that a homogeneous soil profile was assumed in the model while the lysimeter soil actually contained three horizons (silty clay Ap, 0–20 cm; silty clay B1, 20–75 cm; and clay B2, 75–120 cm). This means that the optimized soil hydraulic parameters constitute a compromise between the hydraulic properties of the B1 horizon and the B2 horizon (water content data from shallow depths were not included in the objective function). The hydraulic properties should therefore be viewed as effective properties for the subsoil. In practice, the parameters will mainly represent the B1 horizon, since the hydraulic parameters are most sensitive to the water content changes at shallow depths where the soil is initially relatively dry. Simultaneous optimization of the soil hydraulic properties for the individual soil horizons was not attempted because of the nonuniqueness problems discussed above, and because previous studies indicated that this requires not only (z,t) data, but also h(z,t) and Q(t) data (e.g., Jacques et al., 2002; Ritter et al., 2003).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY HYPOTHETICAL SOIL MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 
Inverse Solutions
The earlier results showed that the parameters _r and need to be fixed to obtain unique soil hydraulic parameter sets. Even then, not all optimizations for the hypothetical soil arrived at the true parameter values; with (z,t) data in the objective function, only 17 out of 35 converged optimizations resulted in parameter sets that were <5% from the true parameter values. The problem now is to estimate _r and for the lysimeter soil. The value of _r was fixed to 0.078, which is an estimate of the hygroscopic water content based on thermogravimetric measurements on undisturbed samples from a soil similar to that of the lysimeter (Kelleners et al., 2004). The stability analysis for the hypothetical soil showed that the estimation of is more critical, so two different values were selected. The value of was fixed to 0.5, as suggested by Mualem (1976), and to –1.055, as predicted by Rosetta. Each of the LM minimizations was repeated 50 times with 50 initial parameter estimates selected randomly from the same predetermined parameter intervals as for the hypothetical soil. These parameter intervals were also used as bounds during the minimization process.

The results of the parameter estimation for L1 and L2 in the lysimeter soil are shown in Table 5. The optimized _s and n values of 0.483 to 0.487 and 1.178 to 1.320, respectively, fall within the expected range for a silty clay soil. The low CV values for _s and n indicate that these parameters were reliably identified. The optimized values for and K_s were 0.002 to 0.004 cm^–1 and 0.1 to 0.32 cm d^–1, respectively. These values are relatively low but not unrealistic. The elevated CV values for and K_s show that there is more uncertainty in these parameters. Also, there is strong correlation between and K_s (correlation coefficient > 0.99). Low CV values for _s and n, high CV values for and K_s, and high parameter correlation between and K_s were also found for the hypothetical soil when (z,t) data were used for fitting (Table 1, _r and fixed).

View this table: [in this window] [in a new window]   Table 5. Results of the soil hydraulic parameter estimation for the water content monitoring locations L1 and L2 in the groundwater table lysimeter.

It is important to note that the converged LM minimizations generally end up at approximately the same values for and K_s, despite the high parameter correlation between these parameters. The sole exception is one solution for L2 with = –1.055 that yielded K_s = 17.52 cm d^–1 and = 0.022 cm^–1. However, the SSQ value of 1.90 for this solution was clearly out of line with all other solutions, and was therefore not included in Table 5. Apparently, the LM minimization ended up in a local minimum in this case. It appears that high parameter correlation does not prevent the LM algorithm from finding the approximate solution, provided multiple initial parameter estimates are used. However, parameter correlation slows down the minimization process in the sense that more iterations are required to arrive at a solution. This is probably the reason why several of the attempted minimizations did not converge within 20 iterations, both for the hypothetical soil and the lysimeter soil.

Soil Hydraulic Properties
The soil hydraulic parameter estimates for = 0.5 were selected for further analysis. We preferred the = 0.5 estimates over the = –1.055 estimates because of the higher number of converged minimizations (Table 5). The calculated soil hydraulic properties for L1 and L2 are depicted in Fig. 7 , together with the Rosetta prediction based on soil texture and dry bulk density data. Two water retention curves and two hydraulic conductivity curves are calculated for each location, representing the upper and lower bounds according to the hydraulic parameter estimates.

View larger version (16K): [in this window] [in a new window]   Fig. 7. Optimized (a) soil water retention and (b) hydraulic conductivity functions for the groundwater table lysimeter for Locations 1 and 2 compared with the Rosetta predicted soil hydraulic functions for a silty clay soil.

Measured vs. Calculated Water Content
Measured vs. calculated water contents with time for all depths included in the objective function are shown in Fig. 8 . The effect of uncertainty on the soil hydraulic parameter estimates is assessed by showing both the lowest and highest calculated water contents at each depth. The speed of the upward moving wetting front is described reasonable well at both locations and does not differ much for the range of parameter estimates. However, the rate at which the water content increases at individual depths is generally overestimated. This may be related to the measurement volume of the capacitance sensor, which was approximated by a vertical zone of influence of 6 cm. Increasing or decreasing the zone of influence in Hydrus did not result in significant improvements (results not shown). Since the measurement volume of the sensor in the presence of a sharp wetting front remains unknown, no further conclusions are warranted.

View larger version (41K): [in this window] [in a new window]   Fig. 8. Measured (symbols) and calculated (lines) volumetric water content for the 79-d fallow period in the groundwater table lysimeter for (a) Location 1 and (b) Location 2. The numbers in the figures refer to the depth below the soil surface of the center of the sensors.

View larger version (15K): [in this window] [in a new window]   Fig. 9. Initial and final measured (symbols) and calculated (lines) volumetric water content with depth for the 79-d fallow period in the groundwater lysimeter for (a) Location 1 and (b) Location 2.