Inverse modeling of a radial multistep outflow experiment for determining unsaturated hydraulic properties

Advances in Water Resources Vol. 22, No. 5, pp. 431–444, 1999 䉷 1999 Elsevier Science Ltd Printed in Great Britain. All rights reserved PII: S 0 3 0 9 - 1 7 0 8 ( 9 8 ) 0 0 0 3 0 - X 0309-1708/99/$ - see front matter

Inverse modeling of a radial multistep outflow experiment for determining unsaturated hydraulic properties Stefan Finsterle* & Boris Faybishenko Lawrence Berkeley National Laboratory, Earth Sciences Division, One Cyclotron Road, Mail Stop 90-1116, CA 94720, Berkeley, USA (Received 4 January 1998; accepted 31 July 1998)

Modeling flow and solute transport in the unsaturated zone on the basis of the Richards equation requires specifying values for unsaturated hydraulic conductivity and water potential as a function of saturation. The objectives of the paper are to evaluate the design of a transient, radial, multi-step outflow experiment, and to determine unsaturated hydraulic parameters using inverse modeling. We conducted numerical simulations, sensitivity analyses, and synthetic data inversions to assess the suitability of the proposed experiment for concurrently estimating the parameters of interest. We calibrated different conceptual models against transient flow and pressure data from a multi-step, radial desaturation experiment to obtain estimates of absolute permeability, as well as the parameters of the relative permeability and capillary pressure functions. We discuss the differences in the estimated parameter values and illustrate the impact of the underlying model on the estimates. We demonstrate that a small error in absolute permeability, if determined in an independent experiment, leads to biased estimates of unsaturated hydraulic properties. Therefore, we perform a joint inversion of pressure and flow rate data for the simultaneous determination of permeability and retention parameters, and analyze the correlations between these parameters. We conclude that the proposed combination of a radial desaturation experiment and inverse modeling is suitable for simultaneously determining the unsaturated hydraulic properties of a single soil sample, and that the inverse modeling technique provides the opportunity to analyze data from nonstandard experimental designs. 䉷 1999 Elsevier Science Ltd. All rights reserved. Keywords: multi-step outflow, unsaturated hydraulic parameters, parameter estimation.

with the methodology used for the subsequent data analysis. We discuss this aspect of parameter estimation in this paper, and examine an experimental design suitable for the determination of unsaturated hydraulic properties, taking advantage of the flexibility provided by inverse modeling techniques. The Tough233 simulator is used for forward modeling, and Itough215 for solving the inverse problem. We focus on the determination of the absolute permeability (k abs, m 2), and the dependence of capillary pressure (p c, Pa) and relative permeability (k r) on saturation (S), believed to be key characteristics affecting model predictions. Note, however, that the relative importance of these and other parameters depends on the intended use of the model, and should be studied prior to performing laboratory or field experiments. The sensitivity of the model

1 INTRODUCTION Modeling transient water flow and contaminant transport in the vadose zone on the basis of the Richards equation requires knowledge of the unsaturated hydraulic conductivity and the water retention characteristics of the soil. The determination of unsaturated hydraulic properties requires two steps: (1) an experimental procedure to obtain relevant data; and (2) a method applied to derive the parameter values from these data. Since both steps are equally important in parameter estimation, the appropriateness of a proposed experimental procedure must be evaluated together *Corresponding author. Tel: 001 510 486 5205; Fax: 001 510 486 5686; e-mail: [email protected] 431

432

S. Finsterle, B. Faybishenko

predictions with respect to the parameters determines the required level of accuracy with which the parameters are to be estimated. Once these requirements have been identified, an appropriate laboratory or fleld experiment can be designed. We first review some of the methods proposed for the estimation of unsaturated hydraulic properties. It is worth mentioning up front that all methods are ‘indirect’ to a greater or lesser extent. For example, even though the determination of the saturated permeability based on flow rate data is considered to be a ‘direct’ measurement, the value is actually obtained by inverting Darcy’s Law, assuming onedimensional, steady-state flow with a constant pressure gradient (i.e. using a specific model of the flow system). The experiment has to be carefully designed to make sure that the assumptions underlying this simple model are satisfied. Any discrepancies between the experiment and the model lead to an error in the permeability estimate. This simple example illustrates that the design of the experiment has to be adapted to the model used for data analysis, or that the model must be modified if the experimental conditions deviate from the simplifying model assumptions. Flexible parameter estimation techniques such as inverse modeling provide opportunities for new experimental designs with higher accuracy and efficiency in the determination of unsaturated hydraulic properties. The following review is loosely arranged from ‘direct’ to more indirect methods, with models of increasing complexity. Probably the most direct method involves pointwise construction of capillary pressure curves by measuring S and p c under equilibrium conditions. For use in numerical models, a functional form is selected and matched to the experimental data. These parametric models can be considered mere fitting functions. However, they contain parameters considered representative of the pore structure. An example is the l parameter of the Brooks–Corey model,2 which can be related to the pore size distribution or fractal dimension of the pore space.32 Consequently, there are attempts to estimate these model parameters from soil texture and other properties that are relatively easy to measure.34,42 A similar procedure is involved when predicting unsaturated hydraulic conductivity from water retention curves. Pore connectivity models relating capillary pressure and relative permeability functions were introduced by Burdine3 and Mualem,29 and are used in the expressions of Brooks and Corey,2 van Genuchten49 and Russo.39 These conventional methods are time-consuming and expensive. Moreover, they do not permit the determination of both relative permeability and capillary pressure curves during a single experiment, which often leads to inconsistent parameters and thus unreliable model predictions.5,27 The disadvantages of the traditional methods were discussed by Salehzadeh and Demond.41 They also reviewed recent measurement techniques and designed a laboratory apparatus for rapid and simultaneous measurement of capillary pressure and relative permeability functions.

While the approaches discussed previously rely on direct observations and make use of geometric models of the pore space, parameter estimation by inverse modeling involves process simulation, i.e. it uses the equations governing flow in unsaturated porous media and determines the parameters by minimizing the differences between model predictions and observed data. An advantage of inverse modeling is that any type of data can be used for parameter estimation, provided that the calculated system response is sensitive to the parameters of interest. Further, numerical simulation and inversion techniques impose few restrictions on the experimental layout and flow processes considered. We focus here on applications of inverse modeling to data obtained under transient, unsaturated flow conditions. Zachmann et al.51 were the first to analyze drainage data using inverse methods. Kool et al.25 and Parker et al.31 applied inverse modeling techniques to one-step outflow experiments. Kool and Parker24 estimated parameters of a hysteretic water retention model based on synthetically generated pressure and saturation data from ponded infiltration and redistribution events, and provided a detailed description of the inverse modeling framework. Problems of non-uniqueness and parameter identiflability were analyzed by Russo et al.40 using infiltration data, by Toorman et al.46 for synthetic data from one-step outflow experiments, and by Simunek and van Genuchten43 for three-dimensional disc permeameter infiltration experiments. These theoretical studies point out the need to include additional data such as tensiometer measurements or prior parameter information to constrain the solution of the inverse problem. Van Dam et al.48 studied data from a one-step axial outflow experiment and suggested adding independent measurements of p c and k r to overcome non-uniqueness problems. Eching and Hopmans11 analyzed one-step and multi-step experiments and found that a simultaneous inversion of cumulative outflow and soil water pressure data gave unique and stable solutions. Applications of inverse modeling techniques to actual field data were described by Dane and Hruska,6 Kool et al.,26 Finsterle and Pruess17 and Zijlstra and Dane.52 The important issue of how unsaturated soil hydraulic functions determined in the laboratory relate to field conditions was discussed by Eching et al.12 The main objectives of this paper are to evaluate the design of a transient, radial, multi-step outflow experiment and to determine unsaturated hydraulic parameters using inverse modeling. Measurements were carried out on undisturbed loamy soil cores10 and analyzed by inverse modeling using Itough2.15 The paper is structured as follows. We first discuss theoretical aspects of the estimation of unsaturated hydraulic parameters, formulate the forward and inverse problems, and introduce the design criteria. Next, we describe the materials and methods used in the experiment. The experimental design is then evaluated based on synthetic data. Finally, we present the analysis of the multi-step drainage experiment and discuss the resulting estimates.

Inverse modeling of multistep outflow experiment function:2

2 THEORETICAL CONSIDERATIONS In this section, we first discuss the model used to solve the forward problem, i.e. the governing equations describing unsaturated flow. Then we formulate the inverse problem, and finally introduce a number of criteria suitable for evaluating experiments designed for parameter estimation purposes. 2.1 Formulation of the forward problem

The accumulation term M represents mass per unit volume (2)

where f is porosity and r is water density. The mass flux term is given by Darcy’s law F ¼ ¹ kabs

kr r (ⵜp ¹ rg) m

1 l pc ¼ ¹ pe ·Se ¹

(3)

Here, k abs (m 2) is the absolute permeability, k r is the relative permeability, a dimensionless number between zero and one as a function of saturation S, m (Pa·s) is the dynamic viscosity, p (Pa) is the soil water pressure, and g (m 2 s ¹1) is the gravitational acceleration vector. This form of the governing equation makes it explicit that the commonly used hydraulic conductivity (m s ¹1) is not only a function of the porous material, but also of fluid density and viscosity which are slightly affected by pressure and temperature changes. Estimating the absolute permeability instead of hydraulic conductivity ensures that these dependencies are properly taken into account. While the effect is minor for water flow under ambient conditions, it may become significant when inverting gas flow data (for an example, see Finsterle and Persoff16), or when performing predictions under strongly nonisothermal conditions such as encountered during remediation of contaminated sites by means of steam flooding. A crucial part of the conceptual model is the choice of the characteristic curves, expressing capillary pressure p c and relative permeability as a function of S. Many consistent parametric models have been proposed in the literature (for a comprehensive review, see Durner8,9). We restrict our discussion to the most widely used models. Brooks and Corey (BC) introduced the following capillary pressure

(4)

Here, p e and l are fitting parameters sometimes referred to as air entry pressure and pore size distribution index, respectively. The effective liquid saturation S e is defined as Se ¼

The general-purpose numerical simulator Tough233 is used to solve the forward problem. While Tough2 is able to handle nonisothermal flow of multiple components in up to three phases, we discuss here only the terms involved in unsaturated flow of water according to the Richards equation. The integral finite difference method employed in Tough2 solves the following mass balance equation for an arbitrary subdomain V n bounded by the surface G n and an inward normal vector n, where t is time and q is a local sink or source term: Z Z Z d M dV ¼ F·ndG þ q dV (1) Gn Vn dt Vn M ¼ fr

433

Sl ¹ Slr 1 ¹ Slr

(5)

where S lr is the residual liquid saturation. Introducing eqn (4) in the model suggested by Burdine3 yields the following expression for relative permeability: 2 þ 3l kr ¼ Se l

(6)

Applying Mualem’s model yields39 kr ¼ Set þ 2 þ 2=l

(7)

where t is an additional fitting parameter which accounts for pore connectivity and tortuosity effects. Note that with t¼1, eqn (7) is identical to eqn (6). We will refer to eqns (5) and (6) as the BCB model, and eqn (5) in combination with eqn (7) as the BCM model. Alternative expressions are given by van Genuchten (VG):49 1 pc ¼ (Se¹ 1=m ¹ 1)1=n a

(8)

With m ¼ 1 ¹ 2/n and applying Burdine’s model, the relative permeability becomes   m (9) kr ¼ S2e 1 ¹ (1 ¹ S1=m e ) Assuming m ¼ 1 ¹ 1/n and introducing eqn (8) into Mualem’s model yields   m 2 (10) kr ¼ Ste 1 ¹ (1 ¹ S1=m e ) We will refer to eqns (8) and (9) as the VGB model, and eqn (8) in combination with eqn (10) as the VGM model. For convenience and based on a weak analogy to the BC model described by Morel–Seytoux et al.,28 we will refer to parameter n as the pore size distribution index (PSDI), 1/a as the air entry pressure (AEP), and t as the tortuosity factor. While the two models, BC and VG, exhibit only minor differences in the capillary pressure function for intermediate and low liquid saturations, the system behavior may differ significantly near full saturation.8,9 Also note that the pore connectivity models by Burdine and Mualem require a single value for the residual liquid saturation S lr. However, the physical meaning of this parameter is ambiguous.30 In the water retention curve, it indicates the asymptotic saturation at which the capillary pressure approaches infinity. In the relative permeability function, it indicates the saturation at which the liquid phase becomes discontinuous, i.e. where saturation cannot be further reduced by

434

S. Finsterle, B. Faybishenko

hydromechanical processes. It is obvious that the values for S lr do not have to coincide in the two different curves. It is important to realize that the parameters estimated by inverse modeling are not intrinsic properties of the soil, but are model parameters strictly related to the specific formulation as outlined in this section. Transferring parameters from one model to another may lead to conceptual and thus numerical prediction inconsistencies. The same problem arises when using ‘directly’ measured parameters values in an analytical or numerical prediction model. On the other hand, the parameters estimated by inverse modeling can be considered optimal for the specific forward model used to simulate the experiment. 2.2 Formulation of the inverse problem Solving the inverse problem is usually defined as the estimation of parameters by calibrating a model against the observed data. In the broader sense of model identification, however, inverse modeling also requires identifying the most suitable conceptual model, which includes the functional form of the characteristic curves. In this section we focus on the parameter estimation procedure. We follow the standard procedure and minimize some measure of the differences between the observed and predicted system responses, which are assembled in the residual vector r with elements ri ¼ yⴱi ¹ yi (p)

(11)

yⴱi

is an observation at a given point in space and Here, time, and y i is the corresponding model prediction, which depends on the vector p of the unknown model parameters. In inverse modeling, the distribution of the final residuals is expected to be consistent with the distribution of the measurement errors, provided that the true system response is correctly identified by the model. If the error structure is assumed to be Gaussian, the objective function to be minimized can be infened from maximum-likelihood considerations to be the sum of the squared residuals weighted by the inverse of the covariance matrix C yy:20 ¹1 Z(p) ¼ rT Cyy r

(12)

An iterative procedure is needed to minimize eqn (12). The Levenberg–Marquardt modification of the Gauss–Newton algorithm1 was found to be suitable for our purposes. Under the assumption of normality and linearity, a detailed error analysis of the final residuals and the estimated parameters can be conducted.17 For example, the covariance matrix of the estimated parameter set is given by: ¹1 ¹1 J) Cpp ¼ s20 (JT Cyy

(13)

where J is the Jacobian matrix at the solution. Its elements are the sensitivity coefficients of the calculated system response with respect to the parameters: Jij ¼ ¹

⳵ri ⳵yi ¼ ⳵pj ⳵pj

(14)

In eqn (13), s20 is the estimated error variance, a goodnessof-fit measure given by s20 ¼

¹1 rT Cyy r M¹N

(15)

where M is the number of observations, and N is the number of parameters. The inverse modeling formulation outlined above is implemented in a computer program named Itough2.15 The selection of the most appropriate model given a set of candidate models will be discussed in Section 2.3. 2.3 Design criteria Prior to testing, design calculations should be conducted by means of synthetic data inversions to evaluate the ability of the proposed experiment to estimate the parameters of interest. Performing synthetic inversions reduces the risk of producing an ill-posed inverse problem when analyzing the experimental data. Design calculations provide insight into the sensitivity of each potential observation with respect to the parameters. Synthetic model calibrations can be performed to assess whether the inverse problem is well-posed. A number of criteria have been proposed to measure the performance of an experimental design.22,23,44,45 In this section we discuss the design criteria used to evaluate the suitability of multi-step laboratory outflow experiments for the determination of unsaturated hydraulic properties. As a general rule, the experiment should provide sufficiently sensitive data so that the parameters can be determined with an acceptably low estimation uncertainty. Furthermore, the inverse problem should be well-posed, leading to a unique and stable solution. It is obvious that highly sensitive data provide the most valuable information about the model parameters. We introduce a dimensionless sensitivity measure for parameter j that consists of a sum of the absolute values for the sensitivity coefficients ⳵y i/⳵p j, scaled by the measurement error jyi and the expected parameter variation jpj : M X ⳵yi jpj (16) Qj ¼ · i ¼ 1 ⳵pj jyi This aggregate measure of sensitivity is similar to the one proposed by Kool and Parker.24 However, we use the sum rather than an integral for both space and time, and scale the sensitivity coefficients by jpj : Note that jyi can be viewed as an actual or required measurement accuracy, and jpj as the expected or attainable parameter uncertainty. We believe that using the discretized form for Q j better reflects the situation where discrete data are available from individual sensors at different locations, and that calibration occurs at selected points in time. In other words, we use the same information that will be available for estimation, which allows us to study how repositioning sampling points affects the overall sensitivity. The scaling of the sensitivity coefficients is necessary to compare the

Inverse modeling of multistep outflow experiment contribution of different data types to the solution of the inverse problem, and to evaluate the sensitivity of parameters that vary over different scales depending on their measurement units and physical nature. It is important to realize that a sensitivity analysis does not address the question whether the parameters of interest can be determined independently and with sufficient accuracy. High sensitivity is only a necessary, but not sufficient condition for a well-posed inverse problem (for a detailed discussion of this point see Finsterle and Persoff16). It is essential to perform synthetic inversions to address issues such as uniqueness, stability, and estimation uncertainty.4,24 Here, we focus on the information provided by the covariance matrix C pp, which contains the variances and correlation structure of the estimated parameter set. A measure of overall estimation uncertainty is given by the trace of the covariance matrix. Since strong parameter correlations usually lead to an illposed inverse problem with large estimation uncertainties, an experiment should be designed such that the key parameters can be identified as independently as possible. To analyze the overall parameter correlation, we define jⴱp as the conditional standard deviation of a single parameter, i.e. the uncertainty of a parameter assuming that all other parameters are fixed, and j p as the joint standard deviation, i.e. the square-root of the diagonal element of C pp. The conditional standard deviation of parameter j is the inverse of the jth diagonal element ¹1 J). We of the Fisher information matrix F ¼ s0¹ 2 (JT Cyy propose to interpret the ratio   jp ⴱ (17) kj ¼ jp j as an aggregate measure of parameter correlation, i.e. of how independently parameter j can be estimated. A value close to one signifies an independent estimate, whereas small values indicate a loss of parameter identiflability as a result of its correlation to other uncertain parameters. As mentioned earlier, inverse modeling provides the optimum parameter set for the given conceptual model without indicating whether the model adequately describes the salient features of the flow system. If competing models have been developed and matched to the data, a criterion is needed to decide which of the alternatives is preferable. A number of tests for model discrimination have been described in the literature.4,39 The simplest test is based on the goodness-of-fit measure given by eqn (15). However, since the match can always be improved by adding more fitting parameters, the criterion should contain the number of parameters to guard against overparameterization. We will use the Akaike information criterion (AIC) for model discrimination tests.4 For normally distributed residuals, AIC can be written as AIC ¼ (M ¹ N)s20 þ loglCyy l þ M·log(2p) þ 2N:

(18)

435

A test should be designed so it produces data that allow discriminating among a set of model alternatives.

3 MATERIALS AND METHODS Various outflow experiment designs for determining unsaturated hydraulic properties have been proposed. Experiments with flow along the core axis were first described by Richards,35 Richards and Moore36 and Gardner.18 An experiment with radial flow geometry was proposed by Richards et al.,38 and was further developed by several investigators.19,21,37 Gardner19 developed three methods for the determination of the unsaturated hydraulic conductivity using radial flow experiments. They include: (1) the constant flux method, in which pressure is recorded in time; (2) the constant pressure method, in which the flow rate is measured; and (3) a method that requires instantaneous injection or removal of a small quantity of water and monitoring the associated pressure change. For these methods, Gardner derived analytical solutions, assuming a ratio of a sample radius to a cup radius of at least 10. Richards and Richards37 and Klute et al.21 developed analytical solutions for radial flow experiments, taking into account the actual ratio between the core radius and porous cylinder radius, as well as the membrane impedance. These solutions, which are analogous to Gardner’s solution for the constant pressure method, assumed a linear relationship between the moisture content and pressure for each pressure step, thus limiting the step size that can be taken without loss of accuracy. Doering7 investigated both a one-step procedure and small increments of pressure changes. Valiantzas47 proposed correction terms to Gardner’s solution for a onestep experiment for several forms of the diffusivity function. It should be noted that gravity was neglected in all the proposed analytical solutions. Elrick and Bowman13 found that a small amount of air could diffuse through the porous membrane when the pressure exceeds 10 kPa. If the extraction cylinder is water filled, the appearance of air changes the flow rate, and this air must be removed by means of a flushing procedure, potentially affecting the boundary conditions. Klute et al.21 discussed the advantages of using a radial instead of the more common axial flow geometry. They pointed out that soil shrinkage during drying is significantly reduced in a design with a central porous cylinder, thus avoiding loss of contact between the sample and the extracting cylinder. Furthermore, the air trapped in the porous cylinder can be easily removed with mininimal disturbance of the boundary condition. As a result of the shorter flow distance, a larger sample volume can be tested in shorter times.19 Fig. 1 shows a schematic of a flow cell apparatus for radial, single- and multi-step desaturation experiments on soil samples. Faybishenko14 and Dzekunov et al.10 have tested the design and performed a series of radial flow experiments using the central porous cylinder for both

436

S. Finsterle, B. Faybishenko and inverse models, correctly representing the conditions encountered by the vertical tensiometer. The samples investigated were saturated under vacuum in order to minimize the effect of entrapped air.50

4 RESULTS AND DISCUSSION 4.1 Evaluation of experimental design

Fig. 1. Schematic of apparatus for radial flow experiment.

injection and extraction of water. A second porous cylinder is used as a monitoring tensiometer. They conducted wetting and drying experiments by applying one-step, multistep, and continuously changing boundary pressures under isothermal and nonisothermal conditions using cores of different sizes. The experimental procedure can be described as follows. Soil cores (22 cm long and 15–18 cm in diameter) are cut in the field using a special cylindrical knife which minimizes the disturbance of the sample.10,14 The cores are conserved in a solid metal or plastic cylinder, and the annulus is filled in with a paraffin–tar mixture. In the laboratory, a hole with a diameter of 2.3 cm is drilled co-axially at the center of the core, in which a ceramic cylinder with an air entry pressure of about 1 bar is inserted. The top and bottom surfaces of the core are covered with end caps. The central porous cylinder is connected to a burette to measure the cumulative water discharge. The burette is connected to a vacuum pump. In order to inhibit air accumulation in the cylinder, which would affect the outflow measurements,13,21 the opening of the porous cylinder is directed downward, and the inner void space is kept in an air-filled rather than water-filled state. Such a design allows the extracted water to freely flow into the measuring burette. A port allows air exchange between the core and the atmosphere. A tensiometer for water potential measurements is inserted near the outer wall of the flow cell. The radial flow geometry allows vertical installation of a tensiometer. As confirmed by our numerical simulations, the pressure measured by the tensiometer is very close to the vertically averaged pressure at any radial distance, i.e. the system behavior can be accurately described by a one-dimensional, radial model. Also note that the flexibility offered by inverse modeling does not require that the average pressure be measured by the tensiometer. For example, for soil cores with higher permeability and weaker capillarity, where gravity effects may become significant, a two-dimensional (r–z) representation of the core can be used in the forward

In this section, we discuss numerical simulations of a multistep, radial outflow experiment and assess the suitability of the design for the estimation of unsaturated hydraulic properties. We model a synthetic multi-step desaturation experiment where the pressure at the center of a cylindrical core is reduced in discrete steps from ¹2 to ¹10, ¹20, ¹30, ¹60 and ¹90 kPa after 1, 2, 3, 5, and 10 days, respectively. It is assumed that the cumulative outflow is recorded as a function of time, and that a tensiometer is installed for capillary pressure measurements. The standard errors of the outflow and pressure measurements are assumed to be 2% of the measured value. Design calculations have to rely on preliminary information about the system. If prior knowledge is poor, the conclusions regarding the optimum design have to be assessed for a number of alternative conceptual models, and a robust design has to be chosen that comprises a wide range of conditions possibly encountered during the experiment. As an example, we present here only the results obtained with the BCB model, which turned out to be a likely candidate model for the soils investigated. The performance of the proposed radial flow experiment is analyzed assuming that three parameters are to be estimated based on capillary pressure and cumulative outflow data. The three parameters are the logarithm of the absolute permeability, log(k abs), the pore size distribution index, l, and the logarithm of the air entry pressure, log(p e). Fig. 2 shows the simulated system behavior calculated with the base-case parameter set given in Table 1. The outflow and capillary pressure curves represent potential data to be used in an inversion for determining the three parameters of interest. While the pressure prescribed at the central extraction cylinder is reduced in discrete steps, the capillary pressure at the outer wall of the flow cell decreases rather smoothly. The cumulative outflow through the extraction cylinder reaches 900 ml at the end of the experiment, draining the sample to about 50% of its initial water content. We first perform a sensitivity analysis to determine the optimal location of the tensiometer. The location within the core sample yielding the maximum total sensitivity provides the most information about the model parameters to be estimated. Simulations were performed using the parameter set in Table 1, and the sensitivity measure (eqn (16)) was evaluated for different potential locations of the monitoring tensiometer along the radius of the core. The dependence of sensitivity on sensor location results from the lack of capillary equilibrium during the transient phase of the

Inverse modeling of multistep outflow experiment

437

Fig. 2. Simulated cumulative outflow through central extraction cylinder and capillary pressure near wall of flow cell. The total pore volume is 1820 ml.

experiment. Fig. 3 shows the sensitivity measure as a function of radial distance of the tensiometer from the core axis. Recall that a pressure boundary condition is prescribed at the ceramic cylinder in the center of the core. Consequently, the sensitivity is zero at the extraction cylinder. Data sensitivity increases and reaches a maximum if the tensiometer is placed at a radial distance of about 0.028 m from the core axis, before it starts to decrease towards the wall of the flow cell. Note, however, that the decrease of sensitivity with radial distance is very minor. Because of the non-linearity of unsaturated flow problems, the curves shown in Fig. 3 depend on the base-case parameter set, i.e. the point of maximum sensitivity may shift if soil properties vary. For example, if the permeability of the sample is higher than expected, the zone of low sensitivity around the core axis is larger. While an optimum design requires the tensiometer to be located relatively close to the center of the core, a robust design suggests its installation near the outer wall of the flow cell to avoid the low sensitivity zone in the case of higher perrneability. If permeability happens to be lower than expected, the sensitivity at the outer boundary decreases slightly, but remains at an acceptable level. Installation of the tensiometer at the cell wall has the additional advantage of minimizing its impact on the water flowing towards the center of the core. As mentioned earlier, strong correlations among the parameters may severely affect the quality of the inverse Table 1. Base case parameter set and expected parameter variation Parameter

Base case value

jp

log(k abs) (m 2) pore size distribution index l log(p e) (Pa)

¹13.00 0.25 3.00

1.00 0.25 0.50

modeling results if random or systematic errors are present. Because parameter correlations are not addressed by a standard sensitivity analysis, design calculations should also include synthetic data inversions to examine the potential estimation uncertainty and correlation structure. The design should then be modified to minimize the trace of the resulting covariance matrix or some other uncertainty measures.44 For problems with few unknown parameters, contouring the objective function eqn (12) is a means to visualize the well- or ill-posedness of the inverse problem. Points of equal objective function lie on continuous surfaces in the parameter space. Fig. 4 shows contour plots in the three parameter planes: (a) log(k abs) ¹ log(p e); (b) log(k abs)¹l; and (c) l ¹ log(p e). The top and middle row of panels show, respectively, the objective function obtained when only flow rate or only pressure measurements are available. The bottom row results from combining the two types of observations. The shape, size, orientation, and convexity of the minimum provides information about the uniqueness and stability of the inversion, and represents the uncertainty and correlation structure of the estimated parameter set. Furthermore, the presence of local minima can readily be detected. The planes shown in Fig. 4 intersect the global minimum. In the absence of measurement errors, the objective function visualized in Fig. 4 is devoid of local minima. The central panel in Fig. 4 reveals that the joint estimation of log(k abs) and l is likely to be unstable if only pressure measurements were available. The combination of pressure and flow rate data yields a well-defined global minimum. Note that the orientation of the minima from the flow rate data tend to be orthogonal to those from the pressure data, i.e. when combined, a well-developed minimum results. In this example, the topology shown in the bottom row in Fig. 4 is similar to that shown in the middle row, because pressure

438

S. Finsterle, B. Faybishenko

Fig. 3. Dimensionless sensitivity of capillary pressure summed over time with respect to log(k abs), l, and log(p e) as a function of tensiometer location. A time-varying suction pressure is prescribed at the extraction cylinder.

data have a higher relative sensitivity given the assumed measurement accuracy of 2%. Synthetic inversions demonstrate that the global minimum is accurately identified by the Levenberg–Marquardt minimization algorithm. Recall that the value of the objective function is usually evaluated only at a few points in the parameter space along the path taken by the minimization algorithm. However, information about the structure of the objective function can be obtained from an analysis of the covariance matrix, which relies on a local examination of the objective function curvature at the minimum. In this case, the covariance matrix calculated at the minimum indicates that it should be feasible to obtain accurate parameter estimates using the multi-step desaturation experiment. 4.2 Analysis of multi-step experiment We analyze a radial, multi-step outflow experiment performed on a soil of loam with a bulk density of 1.44 g cm ¹3, and a porosity of 0.46. The sequence of applied suction pressures was chosen such that the data allow for a pointwise construction of the water retention curve for comparison with inverse modeling results (see discussion of Fig. 8 later). In this experiment, the applied pressure was temporarily increased as shown in Fig. 5, allowing water to redistribute within the core. However, since no water was supplied through the central cylinder, the saturation changes induced by the temporary pressure increase were minor, and hysteretic effects are expected to be too small to be identified. The cumulative outflow data are shown in Fig. 6. A numerical model (see Section 2.1) of the radial outflow experiment was developed, taking into account the flow resistance of the ceramic cylinder as well as gravity effects that turned out to be very small for the soil studied. Inverse modeling was then used to analyze the data for the

parametric models described in Section 2.1. Following the proposed test design (see Section 4.1), the cumulative water discharge through the ceramic extraction cylinder and the capillary pressure near the outer wall of the flow cell were used to estimate absolute permeability and unsaturated hydraulic properties. The number of parameters subjected to the estimation process were varied to study the impact of correlations and the effects of overparametrization. Inversions using the BCM, BCB, VGM and VGB models were performed, estimating between N ¼ 2 and 5 parameters as shown in Table 2. The first inversion (Case A) involved the estimation of five parameters. From this case, we determined the parameters with the lowest sensitivity and largest estimation uncertainty. Residual liquid saturation and tortuosity factor were then removed from the parameter vector by fixing them at the base-case values shown in Table 3. This elimination process based on sensitivity was continued until only two parameters (PSDI and AEP) were estimated in Case E. First, we evaluate the performance of each inversion using the criteria summarized in Table 4. They include the objective function Z (eqn (12)) as a goodness-of-fit measure, the trace of the covariance matrix C pp as an aggregate measure of estimation uncertainty, and the Akaike information criterion (eqn (18)). Because no tortuosity factor t is present in the BCB and VGB model, these models do not appear in Cases A and B. Also note that BCM and BCB become identical in Cases C, D and E by setting t to the base-case value of one. The Z-values in Table 4 confirm the obvious fact that for a given model, the match can always be improved by increasing the number of fitting parameters. Taking the goodnessof-fit as the only selection criterion is thus inappropriate and usually leads to over-parameterization. Note that increasing the number of parameters also increases the correlations among the parameters, which results in higher estimation

Inverse modeling of multistep outflow experiment

439

Fig. 4. Contours of the objective function in the three parameter planes: (a) log(k abs)-log(p e), (b) log(k abs) ¹ l, and (c) l ¹ log(p e). The first row shows the objective function when only cumulative outflow data are used, the second row includes only pressure data, and the third row comprises both the cumulative outflow and pressure data. The planes intersect the parameter space at the global minimum, i.e. they contain the best estimate parameter set.

uncertainties if the match is not significantly improved (see Table 5). The trace of the estimation covariance matrix is an aggregate measure of overall estimation uncertainty, and is, therefore, used as one of the key measures for evaluating inversion results. The models with the lowest AIC values were selected for further discussion. Recall that the AIC takes into account both goodness-of-fit and parsimony of the model. Based on the criteria summarized in Table 4, VGB can be discarded as a candidate model. It fails to match the data, rendering the estimated parameters meaningless. The strong rejection of the VGB model is somewhat surprising given its apparent closeness to the VGM model. However, the two models differ significantly in the predicted relative

permeability near full saturation, making them distinguishable when trying to concurrently match flow and pressure data. On the other hand, both the VGM and BCM models were equally capable of matching the pressure and flow rate data as shown in Figs. 5 and 6. In order for the VGM model to match the data, it is necessary to estimate tortuosity as an additional parameter. If tortuosity is fixed at a value of tVGM ¼ 0:5, then the criteria discussed previously all favor the BCB over the VGM model. In the remainder of this section we discuss the estimated parameters obtained with the BCB model, Case D, and the VGM model, Case B, which realized the lowest AIC. We first note (Table 5) that the estimation uncertainties are reasonable, and that they are lower for the BCB model

440

S. Finsterle, B. Faybishenko

Fig. 5. Comparison between observed (symbols) and simulated (lines) pressure using the BCB (solid), VGM (dash–dotted) and VGB (dashed) model. The prescribed pressure at the extraction cylinder is shown as a dashed line.

than for the VGM model. This is mainly because estimating t as an additional parameter of the VGM model leads to relatively strong correlations among all the parameters, as indicated by the low k-values determined from eqn (17). Comparing the estimates themselves, it seems that the approximate relations for AEP, 1/a ⬇ p e, and PSDI, n ⬇ l þ 1, hold. However, there is a large discrepancy in the estimated k abs values. The VGM model requires a significantly larger absolute permeability to match the outflow data. A high k abs value seems necessary to compensate for the sharp decline of the relative permeability curve near full

saturation. Fig. 7 shows the BCB and VGM relative permeability functions for the respective best-estimate parameter set given in Table 5. Despite the large gap between the two curves, the resulting effective permeability, which is the physical property determining unsaturated flow, is very similar. To illustrate this, we have multiplied the VGM relative permeabilities with the ratio of the estimated absolute permeabilities, yielding almost identical effective permeabilities over the range of saturations encountered during the experiment. Predictions made for unsaturated flow are, therefore, expected to be almost identical for

Fig. 6. Comparison between observed (symbols) and simulated (lines) cumulative outflow using the BCB (solid), VGM (dash–dotted) and VGB (dashed) model.

Inverse modeling of multistep outflow experiment Table 2. Case definition for inversions Case

N

Parameters

A B C D E

5 4 4 3 2

SDI, log(AEP), log(k abs), S lr, t PSDI, log(AEP), log(k abs), t PSDI, log(AEP), log(k abs), S lr PSDI, log(AEP), log(k abs) PSDI, log(AEP)

PSDI, pore size distribution index: l for BC models, n for VG models. AEP, air entry pressure: p e for BC models, 1/a for VG models.

both models, despite the differences in absolute permeability. The two models yield different results only for saturations greater than 0.98 (see Fig. 7). The fact that vastly different k abs values were obtained illustrates the model character of the parameters estimated by inverse modeling. In this case, k abs turns out to be a mere fitting parameter. Saturated permeability depends almost exclusively on the distribution of the largest pores, which are drained immediately without having a significant contribution to the observed cumulative outflow. Since calibration occurs against data obtained under partially saturated conditions, little or almost no information about the presence or absence of macropores is contained in the data. Any deviation from the uni-modal pore size distribution underlying both the BC and VG model will result in a significantly different k abs value. This important point was also discussed by Durner8,9 and prompted Luckner et al.27 to match the predicted and experimental hydraulic conductivities at fluid phase contents less than full saturation. If k abs and water retention characteristics were determined independently, i.e. by a Darcy experiment and a conventional desaturation experiment, any error in the model that predicts relative permeabilities from capillary pressure data would result in a gross error in the effective permeability of unsaturated soils. Such an error is introduced, for example, by macropores not accounted for in a uni-modal pore size distribution model. If estimating k abs along with the unsaturated hydraulic properties, such errors can be partly compensated for. Furthermore, calibrating against both pressure and flow rate data reduces the risk that an error in the pressure data or the hydraulic models used to derive eqns (6), (7), (9) and (10) is forced into an error in the relative permeability curves. A disadvantage of inverse modeling is that the estimated parameters can only be used for predictions of similar Table 3. Base case parameter sets Parameter PSDI log(AEP) (Pa) log(k abs) (m 2) t S lr

BCB/BCM

VGB/VGM

0.25 3.00 ¹13.00 1.00 0.05

1.25 3.00 ¹13.00 0.50 0.05

PSDI, pore size distribution index: l for BC models, n for VG models. AEP, air entry pressure: p e for BC models, 1/a for VG models.

441

Table 4. Goodness-of-fit (Z), trace of estimation covariance matrix (C pp), and Akaike information criterion (AIC) for all inversions Case(N) A(5) B(4) C(4) D(3) E(2)

Model BCM VGM BCM VGM BCB/BCM VGM VGB BCB/BCM VGM VGB BCB/BCM VGM VGB

Z 59.9 59.7 62.1 59.9 62.9 70.4 1823.0 63.4 70.6 1865.3 91.8 339.0 1932.1

Trace (C pp)

AIC

⫻ 10

¹9.8 ¹10.1 ¹9.7 ¹11.8 ¹8.9 ¹1.4 1751.2 ¹10.3 ¹3.1 1791.5 16.9 264.2 1856.3

1.8 1.2 1.0 1.9 3.9 2.9 4.5 5.1 1.3 9.2 2.3 6.5 4.7

¹1

⫻ 10 ¹2 ⫻ 10 ¹2 ⫻ ⫻ ⫻ ⫻ ⫻ ⫻

10 ¹3 10 ¹2 10 ¹1 10 ¹4 10 ¹4 10 ¹3

processes and saturation ranges to those encountered during calibration. It is obvious that for saturated zone modeling, k abs should be determined from an experiment performed under fully saturated conditions. Ideally, we would like to have a single parameter set that can be applied to both the saturated and unsaturated zone. Provided that an independent measurement of k abs is available, this information can be used as an argument to discriminate between the VGM and BCB model, which both perform equally well under unsaturated conditions. For the loamy soils investigated here, saturated permeability is measured to be about 10 ¹13 m 2 or less, i.e. more than an order of magnitude lower than the estimate obtained with the VGM model. Given this information, it is reasonable to favor the BCB model with an estimated k of 5.6 ⫻ 10 ¹14 m 2 over the VGM model, thus expanding the applicability of the parameter set to studies of both the unsaturated and saturated zone. We note in passing that if k abs is fixed at 10 ¹13 m 2, the VGM Model fails to match the data, even though parameters n and 1/a are optimized, whereas an acceptable match can be obtained with the BCB model after optimizing for l and p e (see Case E). Fig. 8 shows that the capillary pressure curves obtained with the two models and the respective best estimates are very similar. The symbols shown in Fig. 8 represent a conventionally determined capillary pressure curve, assuming that equilibrium was reached at the end of each pressure step. The curve was constructed by taking the average saturation calculated from the cumulative outflow data and the corresponding capillary pressure measured with the tensiometer. The equilibrium assumption is not quite fulfilled during this transient experiment, i.e. the pressure at the outer cell wall is in fact not identical to the pressure at the extraction cylinder. This leads to a shift in the capillary pressure curve to the left, and the capillary pressure may be underestimated by up to 15 kPa. Since the actual flow problem is solved in inverse modeling, transient effects are automatically taken into account, yielding a more accurate estimate of the capillary pressure curve.

442

S. Finsterle, B. Faybishenko Table 5. Best estimates and their uncertainty BCB, Case D

VGM, Case B

Parameter

Estimate

Std. Dev.

k

Estimate

Std. Dev.

k

PSDI log(AEP, Pa) log(k abs, m 2) t

0.105 3.05 ¹13.25 1.00

0.003 0.04 0.06 –

0.35 0.29 0.32 –

1.122 3.26 ¹11.90 ¹0.39

0.006 0.06 0.11 1.01

0.16 0.17 0.17 0.12

PSDI – Pore Size Distribution Index: g for BC models, n for VG models. AEP – Air Entry Pressure: pe for BC models, 1/a for VG models.

Fig. 7. Relative permeability function derived by inverse modeling for the Brooks–Corey–Burdine (BCB) and the van Genuchten–Mualem (VGM) model.

5 CONCLUSIONS In this paper, we evaluated the design of a radial, multi-step desaturation experiment on a single core, and analyzed the transient flow and pressure data using inverse modeling

techniques. We determined the absolute and relative permeabilities as well as the capillary pressure function for a loam soil. The inverse code Itough2 used in this study is based on the general-purpose multiphase flow simulator Tough2, which provides the necessary flexibility to

Fig. 8. Capillary pressure functions directly inferred from the data assuming equilibrium, and by inverse modeling using the BCB and VGM models. The direct measurements (symbols) are obtained by volume averaging of saturation using cumulative outflow data and the pressure measured by the tensiometer. The results obtained by inverse modeling (lines) take into account the non-equilibrium and transient effects.

Inverse modeling of multistep outflow experiment analyze data from non-standard, specially designed experiments. These experiments can be optimized to produce data that contain significantly more information about the system to be modeled. We have shown that the radial flow experiment provides sensitive data for the simultaneous determination of absolute permeability and the parameters of the relative permeability and capillary pressure functions. Installation of the tensiometer near the outer wall of the flow cell is a slightly sub-optimal, but robust configuration. ColIecting transient flow data from a multi-step experiment provides the information needed to constrain the effective permeability governing unsaturated flow. In many cases, applying conventional curve fitting procedures to capillary pressure data collected under equilibrium conditions does not allow one to distinguish between alternative conceptual models such as the Brooks–Corey–Burdine or van Genuchten–Mualem model. As a consequence, predictions made with the resulting parameter set may be erroneous if the wrong model is chosen, and if absolute permeability and unsaturated hydraulic properties are determined independently. It is therefore important to numerically simulate a transient experiment, capturing the relevant processes governing unsaturated flow, as opposed to inferring effective permeability from geometric pore size distribution models. We have used three criteria to evaluate inversions that use different conceptual models and have different numbers of adjustable parameters. We have demonstrated that the goodness-of-fit criterion is insufficient and misleading. It has to be complemented by an aggregate measure for estimation uncertainty, and a penalty term to guard against overparameterization. We have pointed out in this paper that the estimated parameters are not intrinsic properties of the porous medium; they are related to the functional model being used as illustrated by the dependence of the absolute permeability estimate on the hydraulic model. If an independent measurement of absolute permeability (or any value of effective permeability) were made, the inverse solution can be further constrained, making it possible to select the model that is more likely to be true. On the other hand, if no such measurement is available, the parameter value concurrently estimated by inverse modeling partly compensates for the error in the model, making the subsequent predictions more accurate. The proposed experimental design and analysis procedure will be used in the future to investigate additional effects caused by temperature changes, entrapped air, anisotropy, and hysteresis.

ACKNOWLEDGEMENTS This work was partially supported by the Environmental Management Science Program under a grant from EM-52, Office of Science and Technology, and Office of Energy

443

Research, of the US Department of Energy under Contract no. DE-AC03-76SF00098. We thank K. Pruess and E. Sonnenthal (LBNL) for their reviews of an earlier draft of this paper. The valuable commments and suggestions of three anonymous reviewers are gratefully acknowledged.

REFERENCES 1. Beck, J. V., and K. J. Arnold, Parameter Estimation in Engineering and Science. Wiley, New York, 1977. 2. Brooks, R. H. and Corey, A. T., Hydraulic properties of porous media, Hydrology Paper, Vol. 3. Colorado State University, Fort Collins, 1964, pp. 1–27. 3. Burdine, N. T. Relative permeability calculations from pore size distribution data. Petr. Trans. Am. Inst. Mining Eng, 1953, 198, 71–78. 4. Carrera, J. and Neuman, S. P. Estimation of aquifer parameters under transient and steady state conditions: 1. Maximum likelihood method incorporating prior information. Water Resour. Res., 1986, 22(2), 199–210. 5. Clausnitzer, V., Hopmans, J. W. and Nielsen, D. R. Simultaneous scaling of soil water retention and hydraulic conductivity curves. Water Resour. Res., 1992, 28(1), 19–31. 6. Dane, J. H. and Hruska, S. in situ determination of soil hydraulic properties during drainage. Soil Sci. Soc. Am. J., 1983, 47, 619–624. 7. Doering, E. J. Soil-water diffusivity by one-step method. Soil Sci., 1965, 99(5), 322–326. 8. Durner, W., Vorhersage der hydraulischen Leitfa¨higkeit strukturierter Bo¨den, Bayreuther Bodenkundliche Berichte, Vol. 20, Bayreuth, Germany, 1991. 9. Durner, W. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res., 1994, 30(2), 211–223. 10. Dzekunov, N. E., Zhemov, I. E. and Faybishenko, B. A., Thermodynamic Methods of Investigating the Water Regime in the Vadose Zone. Nedra, Moscow, 1987 (in Russian). 11. Eching, S. O. and Hopmans, J. W. Optimization of hydraulic functions from transient outflow and soil water pressure head data. Soil Sci. Soc. Am. J., 1993, 47, 619–624. 12. Eching, S. O., Hopmans, J. W. and Wallender, W. W. Estimation of in situ unsaturated soil hydraulic functions from scaled cumulative drainage data. Water Resour. Res., 1994, 30(8), 2387–2394. 13. Elrick, D. L. and Bowman, D. H. Note on an improved apparatus for soil moisture flow measurements. Soil Sci. Soc. Am. Proc., 1964, 28, 450–453. 14. Faybishenko, B. A., Water–salt Regime of Soils Under Irrigation. Agropromizdat, Moscow, 1986 (in Russian). 15. Finsterle, S., Itough2 command reference, Version 3.1, Report LBNL-40041. Lawrence Berkeley National Laboratory, Berkeley, CA, 1997. 16. Finsterle, S. and Persoff, P. Determining permeability of tight rock samples using inverse modeling. Water Resour. Res., 1997, 33(8), 1803–1811. 17. Finsterle, S. and Pruess, K. Solving the estimation– identification problem in twophase flow modeling. Water Resour. Res., 1995, 31(4), 913–924. 18. Gardner, W. R. Calculation of capillary conductivity from pressure plate outflow data. Soil Sci. Soc. Am. Proc., 1956, 20, 317–320. 19. Gardner, W. R., Measurement of capillary conductivity and diffusivity with a tensiometer, Trans. 7th International

444

20. 21. 22.

23.

24. 25.

26. 27.

28.

29. 30.

31.

32. 33.

34.

S. Finsterle, B. Faybishenko Congress of Soil Science, Vol. 1. Madison, WI, Elsevier, Amsterdam, 1960, pp. 300–305. Gill, P. E., Murray, W. and Wright, M. H., Practical Optimization, Academic, San Diego, CA, 1981. Klute, A., Whisler, F. D. and Scott, E. J. Soil water diffusivity and hysteresis data from radial flow pressure cells. Soil Sci. Soc. Am. Proc., 1964, 28, 160–163. Knopman, D. S. and Voss, C. I. Behavior of sensitivities in the advection–dispersion equation: implications for parameter estimation and sampling design. Water Resour. Res., 1988, 24(2), 225–238. Knopman, D. S. and Voss, C. I. Multiobjective sampling design for parameter estimation and model discrimination in groundwater solute transport. Water Resour. Res., 1989, 25(10), 2245–2258. Kool, J. B. and Parker, J. C. Analysis of the inverse problem for transient unsaturated flow. Water Resour. Res., 1988, 24(6), 817–830. Kool, J. B., Parker, J. C. and van Genuchten, M. Th. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation: I. Theory and numerical studies. Soil Sci. Soc. Am., 1985, 49, 1348–1354. Kool, J. B., Parker, J. C. and van Genuchten, M. Th. Parameter estimation for unsaturated flow and transport models — a review. J. Hydrol., 1987, 91, 255–293. Luckner, L., van Genuchten, M. Th. and Nielsen, D. R. A consistent set of parametric models for the two-phase flow of immiscible fluids in the subsurface. Water Resour. Res., 1989, 25(10), 2187–2193. Morel-Seytoux, H. J., Meyer, P. D., Nachabe, M., Touma, J., van Genuchten, M. Th. and Lenhard, R. J. Parameter equivalence for the Brooks–Corey and van Genuchten soil characteristics: preserving the effective capillary drive. Water Resour. Res., 1996, 32(5), 1251–1258. Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res., 1976, 12(3), 513–522. Nimmo, J. R. 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., 1991, 27(4), 661– 662. Parker, J. C., Kool, J. B. and van Genuchten, M. T. Determining soil hydraulic properties from one-step outflow experiments by parameter identification, II. Experimental studies. Soil Sci. Soc. Am. J., 1985, 49, 1354–1359. Perfect, E., McLaughlin, N. B., Kay, B. D. and Topp, G. C. An improved fractal equation for the soil water retention curve. Water Resour. Res., 1996, 32(2), 281–287. Pruess, K., Tough2 — a General-purpose Numerical Simulator for Multiphase Fluid and Heat Flow, Report LBL-29400. Lawrence Berkeley National Laboratory, Berkeley, CA, 1991. Rawls, W. J. and Brakensiek, D., Estimation of soil water retention and hydraulic properties. In Unsaturated Flow in Hydrologic Modeling, Theory and Practice, ed. H. J. MorelSeytoux. Kluwer, Dordrecht, 1989, pp. 275–300.

35. Richards, L. A. Capillary conductance of liquids through porous mediums. Physics, 1931, 1, 318–333. 36. Richards, L. A. and Moore, D. C. Influence of capillary conductivity and depth if wetting on moisture retention in soils. Trans. Am. Geophys. Union, 1952, 33, 531–540. 37. Richards, L. A. and Richards, P. L. Radial-flow cell for soilwater measurements. Soil Sci. Soc. Am. Proc., 1962, 26(6), 515–518. 38. Richards, L. A., Russel, M. B. and Neal, O. R. Further developments on apparatus for field moisture studies. Soil Sci. Soc. Am. Proc., 1937, 2, 55–64. 39. Russo, D. Determining soil hydraulic properties by parameter estimation: on the selection of a model for the hydraulic properties. Water Resour. Res., 1988, 24(3), 453–459. 40. Russo, D., Bresler, E., Shani, U. and Parker, J. C. Analyses of infiltration events in relation to determining soil hydraulic properties by inverse problem methodology. Water Resour. Res., 1991, 27(6), 1361–1373. 41. Salehzadeh, A. and Demond, A. H. Apparatus for the rapid automated measurement of unsaturated soil transport properties. Water Resour. Res., 1994, 30(10), 2679–2690. 42. Saxton, K. E., Rawls, W. J., Romberger, J. S. and Papendick, R.I. Estimating generalized soil-water characteristics from texture. Soil Sci. Soc. Am. J., 1986, 50, 1031–1036. 43. Simunek, J. and van Genuchten, M. Th. Estimating unsaturated soil hydraulic properties from tension disc infiltrometer data by numerical inversion. Water Resour. Res., 1996, 32(9), 2683–2696. 44. Steinberg, D. M. and Hunter, W. G. Experimental design: review and comment. Technometrics, 1984, 28, 71–97. 45. Sun, N. and Yeh, W. Coupled inverse problems in groundwater modeling, 2. Identifiability and experimental design. Water Resour. Res., 1990, 26(2), 2527–2540. 46. Toorman, A. F., Wierenga, P. J. and Hills, R. G. Parameter estimation of hydraulic properties from one-step outflow data. Water Resour. Res., 1992, 28(11), 3021–3028. 47. Valiantzas, J. D. A simple approximation equation to calculate diffusivities from onestep outflow experiments. SSSAJ, 1989, 53, 342–349. 48. van Dam, J. C., Strickler, J. N. M. and Droogers, P., From One-step to Multistep Determination of Soil Hydraulic Functions by Outflow Experiments, Report no. 7. Department of Hydrology, Soil Physics and Hydraulics, Agricultural University, Wageningen, The Netherlands, 1990. 49. van Genuchten, M. Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J., 1980, 44, 892–898. 50. Voorhees, W. B. A vacuum chamber for nondestructive wetting of secondary soil aggregates. Soil Sci. Soc. Am. Proc., 1964, 28, 291–292. 51. Zachmann, D. W., DuChateau, P. C. and Klute, A. The calibration of the Richards flow equation for a draining column by parameter identification. Soil Sci. Soc. Am. J., 1981, 45, 1012–1015. 52. Zijlstra, J. and Dane, J. H. Identification of hydraulic parameters in layered soils based on a quasi-Newton method. J. Hydrol., 1996, 181, 233–250