Parameter estimation of two-fluid capillary pressure–saturation and permeability functions

Advances in Water Resources Vol. 22, No. 5, pp. 479–493, 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 2 5 - 6 0309-1708/99/$ - see front matter

T

C

Parameter estimation of two-fluid capillary pressure–saturation and permeability functions J. Chen, J.W. Hopmans* & M.E. Grismer University of California, Department of Land, Air & Water Resources, Hydrology Program, 123 Veihmeyer, Davis, CA 95616, USA (Received 17 December 1997; revised 17 June 1998; accepted 27 June 1998)

Accurate characterization and modeling of subsurface flow in multi-fluid soil systems require development of a rapid, reproducible experimental method that yields the information necessary to determine the parameters of the capillary pressure–saturation and permeability functions. Previous work has demonstrated that parameter optimization using inverse modeling is a powerful approach to indirectly determine these constitutive relationships for air–water systems. We consider extension of the inverse parameter estimation method to the modified multi-step outflow method for two-fluid (air–water, air–oil and oil–water) flow systems. The wellposedness of the proposed parameter estimation problem is evaluated by its accuracy, uniqueness and parameter uncertainty. Seven different parametric models for the capillary pressure– saturation and permeability functions were tested in their ability to fit the multi-step outflow experimental data; these included the van Genuchten–Mualem (VGM) model, van Genuchten–Burdine (VGB), Brooks–Corey–Mualem (BCM), Brooks–Corey– Burdine (BCB), Brutsaert–Burdine (BRB), Gardner–Mualem (GDM) and Lognormal Distribution–Mualem (LNM) models. The VGB, BCM and BCB models fitted the multi-step outflow data poorly, and resulted in relatively large values for the rootmean-squared residuals. It was concluded that the remaining VGM, LNM, BRB and GDM models successfully characterized the multi-step experimental data for two-fluid flow systems. Although having one parameter less than the other models, the GDM model’s performance was excellent. Finally, we conclude that optimized capillary pressure–saturation functions for the oil–water and air–oil could be predicted from the respective air–water function using interfacial tension ratios. 䉷 1999 Elsevier Science Ltd. All rights reserved.

time-consuming, and can often not be replicated1–3. Prediction methods4,5 are typically based on a simplified conceptual soil pore model for predicting permeability from poresize distribution information obtained from capillary pressure–saturation data. Alternatively, the similarity assumption is applied to predict unknown capillary pressure or permeability functions from known values of easy-to-measure soil or fluid properties6,7. In either case, independent data are required to test the validity of the prediction models. In recent years, the parameter estimation approach in combination with inverse flow modeling has been developed to complement existing measurement and predictive methods for estimation of the constitutive relationships for air–water systems. Advantages of the parameter estimation method include: (1) the experiment for which the modeling is to be conducted may be transient, and thereby faster than

1 INTRODUCTION Successful environmental protection and remediation strategies associated with hydrocarbon contamination of soil and ground water requires modeling of multi-fluid flow and transport in subsurface soil systems. However, the implementation of such models is often hampered by lack of sufficient information regarding the capillary pressure– saturation and permeability functions for the different soil materials. Methods to determine soil capillary pressure and permeability functions are either based on measurements or are partially derived from prediction models. Equilibrium and steady-state experimental methods are usually highly restrictive by constrained initial and boundary conditions, *Corresponding author. Email: [email protected] 479

P

a

480

J. Chen et al.

T

C

P

a

Fig. 1. Experimental setup for multi-step outflow. Pressure transducers T 2 and T 3 measured the wetting (w) and non-wetting (nw) fluid pressure, respectively. Pressure transducer T 1 was used to measure outflow rate.

steady-state methods; (2) boundary and initial conditions of the experiment can be freely chosen; and (3) both the capillary pressure–saturation and permeability function parameters are obtained simultaneously for the same soil sample. Estimation of the constitutive functions for air–water systems using the inverse method has become increasingly attractive8–17. The fundamental assumption of inverse modeling parameter estimation is that the constitutive functions can be described by a parametric model for which the unknown parameters can be estimated by minimization of deviations between observed and predicted state variables such as flux density or capillary pressure. To determine whether the inverse problem can be solved, it must be ‘correctly posed’18. Ill-posed inverse problems may result in no solution, non-uniqueness (more than one solution), or instability (solution is sensitive to small changes in input data). Previous studies have shown that the experimental method plays an important role in determining the wellposedness of the parameter estimation problem and indicated that a minimum amount of data must be collected to characterize the simulated flow process. For example, Gardner’s19 one-step transient outflow method may be an ill-posed parameter estimation problem, yielding a nonunique set of parameters for the constitutive relationships. Van Dam et al.20 suggested that cumulative outflow be measured during multiple outflow steps, so as to yield sufficient information for determination of a unique parameter set using the inverse method. Including capillary pressure measurements during the outflow experiment further resulted in improved parameter sensitivity15. Also, Eching and Hopmans9 concluded that measurement of capillary pressure in addition to outflow measurements provided

adequate information for unique solutions of the inverse parameter estimation problem. Although inverse parameter estimation approaches have been developed and increasingly applied in recent years, their application has been limited to air–water flow soil systems. In the modeling of these soil systems, the traditional Richards’ assumption of neglecting the influence of the air phase on water flow is applied. Here, we expand application of the inverse parameter estimation method to other two-fluid flow systems, including air–oil and oil– water. We examine use of the inverse parameter estimation method, using two-fluid multi-step outflow experiments in conjunction with numerical solution of the two-fluid flow equations, to determine the parameters characterizing the constitutive relationships for multi-fluid flow systems. We first use the van Genuchten–Mualem functions because of their proven fitting ability and widespread use, and their successful application by Eching et al.10 to the results of multi-step outflow experiments for air–water systems. We also evaluate the application of the van Genuchten–Burdine (VGB), Brooks–Corey–Mualem (BCM), Brooks–Corey– Burdine (BCB), Brutsaert–Burdine (BRB), Gardner– Mualem (GDM) and lognormal distribution–Mualem (LNM) relationships in the inverse method. Russo and associates21,14 studied the applicability of various constitutive hydraulic relationships for air–water systems, whereas Demond and Roberts2 compared various parametric models for general two-fluid flow systems. However, to date no information is available on the applicability of different constitutive relationships for parameter estimation of multi-fluid flow systems. Consequently, in addition to evaluating the use of the multi-step outflow method in conjunction with the inverse parameter estimation technique

Pressure–saturation and permeability functions

481

Table 1. Physical properties of fluids at 20⬚C Interfacial tension (N/m) Soltrol 130 (Columbia) Soltrol 220 (Lincoln) Viscosity (Ns/m 2 ⫻ 10 ¹3)

Air–oil 0.0239 a 0.0259 a Oil

Oil–water

Air–water 0.0681 a

0.0259 a 0.036423

C

Air

Water

0.0181 Soltrol 130 Soltrol 220

T

1.00

P

1.44 3.92

a

Density (kg/m 3) 1.20 Soltrol 130 Soltrol 220 a

998

762 803

Independently measured.

for estimating constitutive relationship parameters of multifluid flow systems, we also consider the applicability of several different capillary pressure–saturation and permeability relationships to this approach.

2 MATERIALS AND METHODS 2.1 Experimental apparatus Multi-step transient outflow experiments were modified from Eching and Hopmans9,10 to include measurement of the non-wetting fluid pressure (Fig. 1), and were conducted in a constant-temperature (20⬚C) laboratory22. Soils used in the transient outflow experiments were air-dried, sieved to less than 2 mm in diameter, and uniformly packed in a 6.0 cm inside-diameter and 7.6 cm high brass soil core that was placed in a standard Tempe pressure cell. Two soils, a Columbia fine sandy loam (63.2% sand, 27.5% silt and 9.3% clay) and a Lincoln sand (88.6% sand, 9.4% silt and 2.0% clay) were used in the multi-step experiments. These experiments were conducted using air or oil as the non-wetting (nw) fluid, displacing either water, or oil as the wetting (w) fluid, assuming one-dimensional vertical flow in the soil core. Soltrol 130 and Soltrol 220 were used as the oil phase for the Columbia and Lincoln soils, respectively. The density, viscosity, and interfacial tension values for the experimental fluids are listed in Table 1. A 0.74 cm thick porous ceramic plate at the base of the soil core allowed drainage of the wetting fluid, but was impermeable to the non-wetting fluid. Applied air pressures (expressed in cm water height) for the air–water system were 60, 80, 120, 200, 400, and 700 cm above atmospheric pressure for the Columbia soil. Selected equivalent air pressure steps for air–water of the Lincoln soil were 40, 60, 80, 100, 150, 200 and 400 cm above atmospheric pressure. These steps were smaller than those for Columbia soil, because of Lincoln soil’s coarser texture. For both soil types, applied air pressure

increments were chosen such that drainage volumes were approximately equal for each pressure step. The pressure steps for the other fluid pairs were reduced, to account for differences in interfacial tension values between the tested fluid pair and air–water, thereby creating approximately equal drainage volumes for each wetting fluid during each pressure step. A pressure step increment required from 5 to 36 h to reach equilibrium or near-equilibrium conditions, depending on the pressure of the wetting fluid and soil type. After the drainage rate reduced to near zero, the non-wetting fluid pressure was incrementally increased for the next drainage step. Initially the soil was fully saturated with the wetting fluid, and then subsequently drained to a saturation corresponding to slightly in excess of the nonwetting fluid entry pressure, thereby ensuring continuity of the non-wetting fluid throughout the soil core at the onset of wetting fluid displacement24. Both cumulative outflow from the base of the soil core and pressures of both the wetting and non-wetting fluid at the center of soil core were measured continuously using pressure transducers (T 1, T 2 and T 3 in Fig. 1) coupled to a data logger. Cumulative drainage of the wetting fluid (h w|outflow) as a function of time was measured by the wetting fluid pressure in the burette. The two tensiometers enabled calculation of the capillary pressure (P c) as a function of time from the wetting fluid (P w) and the non-wetting fluid (P nw) pressures in the center of the core. Additional details about the experimental methods are presented by Liu et al.22. 2.2 Governing equations Assuming that the porous medium is nondeformable (constant porosity) and that cross-product permeability terms associated with the viscous drag tensor25 can be neglected, the general form of the two-fluid flow equations (without source–sink terms) is described by the two-fluid, volume-averaged momentum and continuity equations26 qw ¼ ¹

kw ·[ⵜPw þ rw g] mw

(1a)

482

J. Chen et al. ⳵(rw Sw ) þ ⵜ·(rw qw ) ¼ 0 ⳵t

(1b)

k qnw ¼ ¹ nw ·[ⵜPnw þ rnw g] mnw

(1c)

f

f

⳵(rnw Snw ) þ ⵜ·(rnw qnw ) ¼ 0 ⳵t

(1d)

(2)

Assuming one-dimensional vertical flow and that the wetting fluid is incompressible, substitution of eqn (1a) into eqn (1b) yields    ⳵vw ⳵ kw ⳵Pw ¼ þ rw g (3) ⳵t ⳵z mw ⳵z When replacing fluid pressures (P w and P nw) and capillary pressure (P c ¼ P nw ¹ P w) with pressure head, defined by hi ¼ Pi /rH2 O g (i ¼ w, nw or c), and defining the hydraulic conductivity of fluid i (K i, L/T) by ki rH2 O g i ¼ w, nw Ki ¼ mi transient flow of the wetting fluid is described by      ⳵hnw ⳵hw ⳵ ⳵hw Kw ¹ þ1 ¼ Cw ⳵z ⳵t ⳵t ⳵z

T

ÿ

In eqns (1a), (1b), (1c) and (1d), the subscripts w and nw denote the wetting and non-wetting fluids, respectively; P i (i ¼ w,nw) denotes pressure (M L ¹1T ¹2); S i (i ¼ w,nw) is the degree of fluid saturation relative to the porosity f (volumetric fluid content vi ¼ fSi ); q i is the flux density vector (L/T); g is the gravitational acceleration vector; m i and r i denote dynamic viscosity (M L ¹1 T ¹1) and density (M L ¹3), respectively; k i ¼ k rik is the effective permeability tensor (L 2), where k is the intrinsic permeability (L 2) and k ri ¼ k ri(S i) is the relative permeability. Assuming that the porous media is nondeformable implies that Sw þ Snw ¼ 1 or vw þ vnw ¼ f

yields the flow equation for the non-wetting fluid:

(4)

(5)

where the capillary capacity C w ¼ dv w/dh c is negative 27. If the non-wetting fluid is compressible (e.g. air is the nonwetting fluid), substitution of eqn (1c) into eqn (1d) yields    ⳵(rnw vnw ) ⳵ rnw knw ⳵Pnw þ rnw g ¼ (6) ⳵z ⳵t ⳵z mnw For air as the non-wetting fluid, it is furthermore assumed that its density is a linear function of the pressure head 28, or   ro, nw (7) hnw rnw ¼ ro, nw þ ho where r o,nw and h o are the reference density and pressure head (at atmospheric pressure), respectively. The ratio r o,nw/h o is defined as the compressibility l (l ¼ 1.24 ⫻ 10 ¹6 g/cm 4). Subsequently, using relationships eqns (4) and (7), and writing water potential in terms of pressure head27


⳵hnw ⳵h þ rnw Cw w f ¹ vw l ¹ rnw Cw ⳵t ⳵t    ⳵ ⳵hnw r þ nw ¼ r K rH 2 O ⳵z nw nw ⳵z

C

ð8Þ

Eqns (5) and (8) can be solved simultaneously for the unknown pressure heads h w and h nw, and thus for h c as well. The upper boundary condition at the top of the soil core (z top) to be used in the inverse modeling of two-fluid flow is described by a zero wetting fluid flux and by a sequence of imposed non-wetting fluid pressures, defined by P nw(T j), or qw (ztop , t) ¼ 0

(9a)

hnw (ztop , t) ¼ hnw (Tj ) ¼

Pnw (Tj ) j ¼ 1, 2, …, M rH2 O g

(9b)

where M is the total number of pressure steps used during time period T j in the multi-step experiment. At the lower boundary of the core (z bottom), the flux of the non-wetting fluid is zero since the entry value of the ceramic plate for the non-wetting fluid is greater than any of the imposed pressures. The wetting fluid pressure at z bottom is determined by the height of the wetting fluid in the burette, which is a function of time as the wetting fluid drains into the burette (Fig. 1), or qnw (zbottom, t ) ¼ 0 hw (zbottom , t) ¼

(10a)

rw h (t) rH2 O wjoutflow

(10b)

Finally, the initial condition for the multi-step method is given by the static pressure of each fluid at time t 1: r z (11a) hw (z, t1 ) ¼ hw (zbottom , t1 ) ¹ w rH 2 O hnw (z, t1 ) ¼ hnw (ztop , t1 ) þ

rnw (z ¹ z) rH2 O top

(11b)

where z is assumed positive upwards and z ¼ 0 at the bottom of the ceramic plate. h w(z bottom,t 1) and h nw(z top,t 1) denote the boundary pressure head values at the start of the transient outflow experiment. 2.3 Constitutive relationships—parametric models Parameter estimation requires the functional description of the capillary pressure–saturation, h c(S w), and permeability functions, k i(S w), in eqns (5) and (8). The choice of a suitable parametric model involves collection of available candidate models, followed by a selection of the appropriate models using techniques to be described later. The constitutive relationships considered below are listed in Table 2.

P

a

Pressure–saturation and permeability functions Perhaps one of the most commonly used models of capillary pressure–saturation function is the van Genuchten29 equation (VG) Sew ¼ [1 þ (avg hc )n ] ¹ m

Sew ¼ (he =hc )l , hc ⬎ he

(13a)

Sew ¼ 1, hc ⱕ he

(13b)

where h e and l are characteristic soil parameters with l characterizing the pore-size distribution and h e equal to the non-wetting fluid entry pressure. The VG model is asymptotically equal to the BC model in the dry range when h c becomes large, as then S ew ⬇ [(1/a)/h c] mn, such that 1/a ¼ h e and mn ¼ l. The lognormal distribution (LN) model31,32 is based on the assumption that soils are represented by a lognormal pore-radius distribution yielding a two-parameter model for the capillary pressure–saturation function 



lnhm ¹ lnhc (14) j   p „ x y2 exp ¹ where Fn (x) ¼ ð1= 2pÞ ¹ ⬁ 2 dy is the normal distribution function that can also be calculated from p Fn (x) ¼ 1=2erfc(x= 2). In eqn (14), the parameters h m and j are physically based and are related to the geometric mean and variance of the pore-size distribution. The Brutsaert (BR) capillary pressure–saturation function 33 is of the form Sew ¼ Fn

b Sew ¼ b þ hgc

include Burdine4 and Mualem34

T

2 ZS 3 e dSe 6 0 27 h 7 6 Burdine (B) : krw ¼ S2ew 6 Z1 c 7, 4 dSe 5

(12)

where S ew denotes the effective saturation of the wetting fluid, Sew ¼ ðvw ¹ vwr Þ/(v ws ¹ v wr), where v ws and v wr are the saturated and residual wetting fluid saturation, respectively; a vg and n are fitting parameters, that are inversely proportional to the non-wetting fluid entry pressure value and the width of pore-size distribution, respectively. In subsequent analysis, we assume that m ¼ 1 ¹ 1/n, and that the effective saturation of the non-wetting fluid (S en) is derived from Sen ¼ 1 ¹ S ew. The Brooks and Corey (BC) model30 is based on empirical observations that log(S ew) is linearly related to log(h c), yielding a power form of S ew(h c)

483

2 Z1 3 dSe 6 S 27 6 e h 7 krn ¼ (1 ¹ Sew )2 6 Z1 c 7 4 dSe 5 0 h2 c

2 ZS 32 e dSe 6 0 7 h 7 6 Mualem (M) : krw ¼ Shew 6 Z1 c 7 , 4 dSe 5 0 hc 2 Z1 32 dSe 6 S 7 6 e h 7 krn ¼ (1 ¹ Sew )h 6 Z1 c 7 4 dSe 5 0

C

0 h2c

P

ð16Þ

ð17Þ

hc

Combining the van Genuchten capillary pressure–saturation eqn (12) with the Mualem (VGM) model yields permeability functions as defined by 35     1 m 2 kw h m ¼ Sew 1 ¹ 1 ¹ Sew krw ¼ (18a) k   1 2m knw h m ¼ (1 ¹ Sew ) 1 ¹ Sew : krnw ¼ k

(18b)

Similarly, substituting eqn (13) into eqn (17) and integrating yields the BCM formulation 34 h þ 2 þ 2=l krw ¼ Sew

(19a) 2 1þ

6 krn ¼ (1 ¹ Sew )h 41 ¹ Sew

3 1 2 l7 5

(19b)

Similarly, the BCB model is derived by substituting eqn (13) into eqn (16) and integrating30

(15) 3þ

where b and g are fitting parameters characterizing the porous medium. While the capillary pressure–saturation function can be considered a static soil property, the permeability function is a hydrodynamic property describing the ability of the soil to conduct a fluid. Capillary pressure–permeability prediction models were developed from conceptual models of flow in capillary tubes combined with pore-size distribution knowledge derived from the capillary pressure–saturation relationship. Typical representations of this type of model

krw ¼ Sew

2 l

(20a) 2 1þ

6 krn ¼ (1 ¹ Sew )2 41 ¹ Sew

3 2 l7 5

(20b)

Comparison of eqns (18) and (19), indicates that the BCB model differs from the BCM model only by the tortuosity parameter h. The VGB model is derived by substitution of eqn (12)

a

484

J. Chen et al.

into eqn (16) and integrating to obtain36   1 m )m krw ¼ S2ew 1 ¹ (1 ¹ Sew

(21a)

  1 m m krn ¼ (1 ¹ Sew )2 1 ¹ Sew

(21b)

Similarly, the LNM model is derived by substitution of eqn (14) into eqn (17) and integrating to obtain 32    2 krw ¼ Shew Fn Fn¹ 1 (Sew ) þ j

(22a)

   2 krn ¼ (1 ¹ Sew )h 1 ¹ Fn Fn¹ 1 (Sew ) þ j

(22b)

where F¹1 n (S ew) is the inverse of F n(S ew). No closed-form Burdine or Mualem permeability models are available for the BR capillary pressure–saturation model eqn (15). However, through simplifying eqn (15), Brutsaert33 obtained a closed form Burdine model (BRB) krw ¼ S2ew [1 ¹ (1 ¹ Sew ) krn ¼ (1 ¹ Sew )

2 1¹ g ]

2 3¹ g

Finally, combining the Gardner krw ¼ e

(23a) (23b) 37

equation

¹ ag hc

(24a) 21

with eqn (17) and using h ¼ 0, Russo pressure–saturation relation (GD)   1 ¹ a g hc 1 Sew ¼ 1 þ ag hc e 2 2

derived a capillary

(24b)

after which also the non-wetting permeabilitity relationship can also be derived (GDM) 0 12 1 ¹ ag hc (24c) krn ¼ @1 ¹ e 2 A In subsequent analyses, for convenience in modeling and due to lack of information to the contrary, we set the h value equal to 0.5 for both the wetting and non-wetting permeability expressions34, despite possible physical reasoning suggesting that h values for wetting and nonwetting fluids may differ38. 2.4 Numerical solution of governing equations and optimization method The governing eqns (3) and (6), or eqns (5) and (8) if fully expressed in pressure head form, in combination with the boundary and initial conditions eqns (9), (10) and (11), with the constitutive relationships described in Table 2, determine the mathematical model of the experimental system. No analytical solution is possible because of the nonlinearity of the constitutive functions. Therefore, we adapted the two-fluid flow numerical model of Dr J.L. Nieber,

University of Minnesota (personal communication) to simulate the transient behavior of two-fluid phase flow as measured in the multi-step outflow experiments. We used this numerical scheme, including the modified Picard linearization algorithm of the mixed-form governing equations and a lumped finite-element approximation39,28. The inverse parameter estimation method is a nonlinear optimization problem, where a vector b containing the unknown parameters of the constitutive relationships is estimated by minimizing an objective function O(b), containing deviations between observed and predicted system response variables for the considered time period. For a given true predictive model (eqns (5) and (8)), formulation of the objective function (OF) using maximum likelihood (ML) considerations yields parameter estimates with non-zero values of the objective function being attributed to measurement errors40,12. Using the notation v as the vector containing the elements of measured flow variables (capillary pressure head, outflow) with the subscript m and s denoting measured and simulated values, respectively, the observation error is given by e ¼ v m ¹ v s(b), with the assumption that model errors are zero. When observation errors are assumed to describe a multivariable normal distribution with zeromean and covariance matrix V ¼ E(v m ¹ v s)(v m ¹ v s) T, the ML estimation results in the formulation of the negative log likelihood function, L 45 n 1 1 L ¼ ln2p þ ln(detV) þ eT V ¹ 1 e (25) 2 2 2 which after removing the constant terms leads to the objective function, O(b), to be minimized: O(b) ¼ eT V ¹ 1 e

(26)

Thus, for a general covariance matrix V, eqn (26) represents a general least squares problem, where the inverse of the error covariance matrix is the weighting matrix that includes measurement accuracy and correlation. In this study, we assume that the measurement errors of e are independent of each other and that their variances are proportional to the magnitudes of the mean values of particular measurements, so that V is a diagonal matrix12. In addition, weighting factors were chosen to be inversely proportional to their mean measured values of cumulative drainage (Q), capillary pressure head (h c) and initial water content v(h c,1,t 1), yielding a weighted least squares (WLS) problem of the form O(b) ¼ WQ

N X

{qj [Qm (tj ) ¹ Qs (tj , b)]}2

j¼1

þ Whc

M X

{qk [hc, m (tk ) ¹ hc, s (tk , b)]}2

k¼1

þ Wv

L X

{q1 [vm (hc, 1 , t1 ) ¹ vs (hc, 1 , t1 , b)]}2

ð27Þ

l¼1

or using W ¼ V ¹1, O(b) ¼ [vm ¹ vs (b)]T W[vm ¹ vs (b)] ¼ eT We

(28)

T

C

P

a

Pressure–saturation and permeability functions

485

Table 2. Two-fluid capillary pressure and permeability models Parameters VGM & VGB (van Genuchten– Mualem/Burdine)

v ws,v wr,k a vg n h

a

Capillary Pressure Function Sew ¼

T

Permeability Function (1) VGM: (m ¼ 1 ¹ 1/n)

1 ½1 þ (avg hc )n ÿm

C

1

h m m 2 krw ¼ Sew ½1 ¹ ð1 ¹ Sew Þ ÿ 1

m 2m krn ¼ (1 ¹ Sew )h ½1 ¹ Sew ÿ

P

a

(2) VGB: (m ¼ 1 ¹ 2/n) 1

2 m Þm ÿ krw ¼ Sew ½1 ¹ ð1 ¹ Sew 1

m m krn ¼ (1 ¹ Sew )2 ½1 ¹ Sew ÿ

 BCM & BCB (Brook and Corey– Mualem/Burdine)

v ws,v wr,k he l h

Sew ¼

he hc

l (3) BCM: h þ 2 þ 2=l krw ¼ Sew 1 þ 1=l 2 krn ¼ (1 ¹ Sew )h ½1 ¹ Sew ÿ

(4) BCB: 3 þ 2=l krw ¼ Sew 1 þ 2=l krn ¼ (1 ¹ Sew )2 ½1 ¹ Sew ÿ

LNM (Lognormal Distribution–Mualem)

v ws,v wr,k hm j h

h krw ¼ Sew {Fn [Fn¹ 1 (Sew ) þ j]}2

Sew ¼ Fn ½lnðhm =hc Þ=jÿ

krn ¼ (1 ¹ Sew )h {1 ¹ Fn [Fn¹ 1 (Sew ) þ j]}2 !   x2 1 x exp ¹ dx ¼ erfc p ¹⬁ 2 2 2

Zx

1 Fn (x) ¼ p 2p b b þ hgc

2 krw ¼ Sew ½1 ¹ ð1 ¹ Sew Þ1 ¹ 2=g ÿ

BRB (Brutsaert–Burdine)

v ws,v wr,k b g

Sew ¼

GDM (Gardner–Mualem)

v ws,v wr,k ag

Sew ¼ ð1 þ 12ag hc Þe ¹ 2 ag hc

a

krn ¼ ð1 ¹ Sew Þ3 ¹ 2=g 1

krw ¼ e ¹ ag hc 1

krn ¼ ð1 ¹ e ¹ 2 ag hc Þ2

In parameter optimization the following parameters were fixed: v ws (measured value) and h ¼ 0.5.

In eqn (27), the subscripts m and s denote measured and simulated values, respectively, and N, M and L refer to the number of observations of Q, h c and v, respectively. The objective function, O(b), is normalized by the weighting factors: WQ ¼ 1

(29a)

12 N 1X Qm (tj ) C B C B N j¼1 C B Whc ¼ B C M X A @ 1 hc, m (tk ) M k¼1

(29b)

0

12 N 1X Qm (tj ) C B C B N j¼1 C B Wv ¼ B C L X A @1 vm (hc, l , tl ) L l¼1 0

(29c)

where L ¼ 1 denotes v m(h c,1,t 1), which is the measured initial volumetric wetting fluid content value corresponding with the initial capillary pressure head value (h c,1) at hydraulic equilibrium before the first pressure step was applied and estimated from cumulative outflow after saturation. Weighting factors q j, q k and q l were arbitrarily set to 1, 1 and 15, respectively. The large value of q l provides a disproportionally large weight to v m(h c,1,t 1), thereby forcing

486

J. Chen et al.

T

C

P

a

Fig. 2. Comparison of measured with optimized h c and Q values of Lincoln soil for (a) air–water, (b) oil–water and (c) air–oil systems.

the optimized capillary pressure–saturation curve to pass through this measured point of the capillary pressure curve. The objective of the optimization procedure is to determine the parameter vector b that minimizes eqn (27). The objective function O(b) is a nonlinear function of b, so that minimization calculations must be carried out iteratively until pre-defined convergence criteria were satisfied. A commonly applied criterion is based on the relative decrease of RMS (root of mean squared residuals) value given by r O(b) (30) RMS ¼ MþN þL

The Levenberg–Marquardt optimization algorithm41 was applied to minimize eqn (27) using the optimization program of Eching and Hopmans9, which was originally developed by Kool et al.42, but was modified to interface with the two-fluid flow model. The parametric models of Table 2 were evaluated based on: (1) comparison of their respective RMS values; (2) posteriori testing of the uniqueness of parameter sets; and (3) model simplicity (in terms of number of fitting parameters). Criteria that consider the number of fitting parameters include the Akaike information criterium (AIC) and the Bayesian information criterium (BIC), which are defined

Pressure–saturation and permeability functions

487

Table 3. RMS values from parameter optimization of air–water (aw), air–oil (ao) and oil–water (ow) soil systems for the seven listed constitutive relationships Parametric model VGM VGB BCM BCB LNM BRB GDM a

Lincoln soil

T

Columbia soil

aw

ao

ow

Total

aw

ao

ow

Total

C

1.66 1.62 5.01 5.19 2.02 1.56 2.46

3.33 4.05 9.68 9.95 4.28 4.55 5.88

2.67

7.66 – 20.75 21.81 9.10 8.59 12.78

2.92 2.57 5.85 5.94 3.54 2.95 5.05

3.85 4.57 9.51 9.61 4.29 2.67 4.68

3.16

9.93 – 19.44 20.89 11.97 9.27 15.19

P

a

6.06 6.66 2.80 2.48 4.44

a

4.08 5.34 4.14 3.65 5.46

No convergence.

by 21,40 AIC ¼ N{ln(2p) þ ln[SSR=(N ¹ p)] þ 1} þ p

(31)

and BIC ¼ (AIC ¹ 2p) þ plnN

(32)

where N and p are the total numbers of observation points and adjustable parameters, respectively, and SSR denotes the sum of squared residuals, numerically equal to O(b). The most parsimonious model minimizes both AIC and BIC, i.e. with their values for the considered fitting models being equal, the model with the least number of fitting parameters is the most acceptable. Hence, the criteria in eqns (31) and (32) penalize for additional fitting parameters. 2.5 Scaling methodology Scaling of interfacial surface tension values provides an additional means by which to evaluate the accuracy of the optimized capillary pressure–saturation functions. In addition to a rigid porous media matrix, the Leverett scaling approach6 assumes that fluids are held in the porous matrix by capillary forces only. In this approach the dimensionless function, J(S ew) is introduced to describe the similarity relationship between two arbitrary soil systems  1  1 Pc, 1 k1 2 Pc, 2 k2 2 ¼ (33) J(Sew ) ¼ j 1 f1 j 2 f2 where j and k denote interfacial tension and intrinsic permeability, respectively. In eqn (33), (k/f) ¹1/2 can be regarded as a microscopic length that depends on the pore geometry only. Therefore, for the same soil matrix with different fluid pairs 1 and 2, eqn (33) reduces to:   j2 (34) h (S ) hc, 2 (Sew ) ¼ j1 c, 1 ew Thus, using eqn (34), the capillary pressure–saturation relationship for a particular fluid pair (e.g. air–water) can be used to predict those for other fluid pairs using the appropriate interfacial tension values. Correct application

of eqn (34) further requires that: (i) the soil is completely water-wet; (ii) the oil phase is the intermediate wetting fluid between the water and air; and (iii) there is no fluid entrapment. It is also assumed that the contact angle is independent of fluid type, the veracity of which is under debate43. The optimized capillary pressure h c(S w) functions of air– oil, and oil–water systems were scaled using the air–water system as a reference. Hence, scaled h c(S w) curves for the air–oil (ao) and oil–water (ow) systems were determined through multiplication of the optimized h c value by j aw/j ao and j aw/j ow, respectively, based on the j values listed in Table 1.

3 RESULTS AND DISCUSSION The inverse parameter estimation procedure was applied to multi-step outflow data for the three two-fluid pairs (air– water, air–oil, and oil–water) and both soils, and parameters were optimized for all seven capillary pressure–saturation and permeability relationships listed in Table 2. Figs 2 and 3 illustrate the comparison of measured with optimized cumulative outflow and capillary pressure heads for the air–water (aw), air–oil (ao), and oil–water (ow) systems of the Lincoln and Columbia soil, respectively, using the VGM constitutive relationships. These are examples of the type of agreement between optimized and measured values achieved with four of the constitutive functions, VGM, LNM, BRB, and GDM. Optimized simulations using the remaining three relationships (VGB, BCM, and BCB) fitted the multi-step data poorly, and are not shown. The sum of RMS values for each two-fluid flow system, constitutive model and soil are listed in Table 3. The relatively small RMS values in Table 3 for the VGM, LNM, BRB and GDM models, as compared to those for the VGB, BCM and BCB models, further supported our observations of the quality of fit between measured and simulated results. The parametric models considered here required six (GDM) or seven parameters. The choice of how many and which parameters to optimize was based on several considerations. First, it must be recognized that as the number of parameters to be optimized increases, the parameter estimation procedure generally leads to better prediction of

a

488

J. Chen et al.

T

C

P

a

Fig. 3. Comparison of measured with optimized h c and Q values of Columbia soil for (a) air–water, (b) oil–water and (c) air–oil systems.

measured flow variables using these optimized parameters. However, this occurs at the expense of the uniqueness of the optimized solution and a corresponding increase in the uncertainty of the parameter estimates. Thus, it is preferred to minimize the total number of fitting parameters. If a parameter can be measured independently with relatively high accuracy, it should be fixed rather than obtained through optimization. Therefore, the saturated wettingfluid content (v ws) obtained experimentally for each soil core from cumulative outflow and initial fluid saturation22 was fixed to measured values of 0.45 cm 3 cm ¹3 (Columbia) and 0.33 cm 3 cm ¹3 (Lincoln). To further reduce the number of optimizable parameters, we set h ¼ 0.534 for the permeability functions. Though the intrinsic permeability

(k) could be measured independently, we found that fixing k overconstrained the permeability relationship such that its shape was determined by only a single fitting parameter. Moreover, differences in pore geometry between the larger pores (defining soil structure) and soil matrix pores (defined by soil texture) resulting in bi-modal pore-size distributions may preclude use of a intrinsic permeability value determined from measurements of K for a saturated soil2,44. Consequently, the constitutive relationships (Table 2) were estimated using four fitting parameters (k, v r, and two additional capillary pressure parameters), with the exception of the GDM model that is described by three fitting parameters (k, v r and a g). Others9,15 have demonstrated that this chosen set of fitting parameters are most sensitive to the

Pressure–saturation and permeability functions

489

Table 4. Final averaged optimized parameters with their uncertainties expressed by the NSD in parentheses (%) for hydraulic functions VGM, LNM, BRB, and GDM of Lincoln (Lin) and Columbia (Col) soils Soil Lin

Model VGM

LNM

BRB

GDM Col

VGM

LNM

BRB c

GDM

Parameter a 3

Air–water ¹3

v r (cm cm ) k (cm 2 ⫻ 10 ¹9) a vg (cm ¹1) n (–) v r (cm 3 cm ¹3) k (cm 2 ⫻ 10 ¹9) h o (cm ¹1) j (–) v r (cm 3 cm ¹3) k (cm 2 ⫻ 10 ¹9) b (m ¹1) g (–) v r (cm 3 cm ¹3) k (cm 2 ⫻ 10 ¹9) a g (cm ¹1) v r (cm 3 cm ¹3) k (cm 2 ⫻ 10 ¹9) a vg (cm ¹1) n (–) v r (cm 3 cm ¹3) k (cm 2 ⫻ 10 ¹9) h o (cm ¹1) j (–) v r (cm 3 cm ¹3) k (cm 2 ⫻ 10 ¹9) b (m ¹1) g (–) v r (cm 3 cm ¹3) k (cm 2 ⫻ 10 ¹9) a g (cm ¹1)

0.0210 (14.0) 24.8 (7.9) 0.0189 ( ⬍ 0.001) 2.811 (1.4) 0.0085 (–) 38.4 (10.0) 66.0 (1.2) 0.77 (1.3) 0.025 (10.0) 38.4 (6.2) 0.004 (8.5) 2.42 (1.1) 0.008 (27.2) 50.9 (12.1) 0.0526 (1.3) 0.0913 (8.7) 5.3 (9.7) 0.0100 ( ⬍ 0.001) 2.145 (2.3) 0.072 (11.7) 6.03 (13.5) 164.8 (3.4) 1.11 (2.4) 0.134 (3.0) 27.8 (17.2) 0.0031 (18.1) 2.120 (1.8) 0.135 (3.4) 10.3 (17.1) 0.0209 (1.8)

Oil–water

Air–oil b

0.00715 (–) 9.2 (10.0) 0.0358 (2.8) 3.137 (2.8) ⬍ 0.001 (–) b 16.9 (14.0) 31.4 (1.6) 0.68 (1.6) ⬍ 0.001 (–) b 13.6 (7.0) 0.017 (9.0) 2.67 (1.3) ⬍ 0.001 (–) b 13.8 (15.7) 0.1011 (1.7) 0.0723 (8.3) 8.0 (10.6) 0.0239 ( ⬍ 0.001) 2.037 (1.5) 0.068 (10.0) 5.17 (14.5) 53.3 (2.2) 0.73 (2.1) 0.093 (4.0) 9.7 (9.3) 0.0032 (18.2) 2.68 (1.3) 0.066 (8.3) 6.4 (15.1) 0.0661 (2.4)

⬍ 0.001 (–) 10.9 (17.0) 0.0529 (1.9) 3.002 (2.5) ⬍ 0.001 (–) b 4.2 (14.8) 24.5 (0.85) 0.382 (3.2) ⬍ 0.001 (–) b 5.0 (14.0) 0.002 (47.0) 3.44 (4.7) ⬍ 0.001 (–) b 7.2 (20.6) 0.1450 (3.8) 0.0613 (9.5) 6.2 (10.4) 0.0253 ( ⬍ 0.001) 2.694 (1.8) 0.083 (8.8) 19.6 (19.0) 72.8 (3.5) 1.03 (2.2) 0.13 (3.4) 81.6 (24.0) 0.02 (17.5) 2.1 (2.2) 0.145 (2.9) 55.1 (27.0) 0.0485 (3.3)

a

v s is fixed for both soils; v s ¼ 0.33 cm 3 cm ¹3 for Lincoln and 0.45 cm 3 cm ¹3 for Columbia soil. Optimized parameter is close to zero, thereby providing a meaningless NSD value. c Not all optimizations converged for air–oil system. b

experimental data obtained from the multi-step outflow experiment when cumulative outflow measurements are supplemented with capillary pressure measurements. Since the VGB model failed to converge for the oil–water system data of both soils, and the BC-type models (BCM and BCB) resulted in RMS values twice as large as that for the remaining models, we eliminated these three models from further evaluation. A possible reason for the poor predictive behavior of the BC-type models in the multi-step outflow application is that experimental capillary pressures values were always larger than the non-wetting-fluid entry

pressure, thereby reducing the significance of the displacement pressure parameter of these models in the optimization. It was surprising to find that the GDM model, with only three fitted parameters is able to describe the soil’s capillary pressure and permeability relations almost as well as the other parametric models, which have one additional fitting parameter. The uniqueness of the optimized parameter sets was tested using three different sets of initial parameter values. Initial parameter estimates were chosen to cover a wide range, but were somewhat constrained to preserve their physical meaning. For example, minimum values of both

Table 5. AIC and BIC values for selected models Soil

Lincoln

Columbia

Model

VGM LNM BRB GDM VGM LNM BRB GDM

Air–water

Air–oil

T

Oil–water

AIC

BIC

AIC

BIC

AIC

BIC

1773 1929 1717 2136 1218 1311 1215 1483

1790 1943 1734 2149 1232 1325 1229 1493

1758 1841 1876 2039 1265 1311 1136 1351

1773 1856 1891 2050 1278 1325 1150 1361

1647 1688 1550 1960 1246 1377 1320 1508

1663 1703 1565 1972 1260 1391 1334 1518

C

P

a

490

J. Chen et al.

T

C

P

a

Fig. 4. Optimized constitutive functions (VGM) of Lincoln soil corresponding to the parameters listed in Table 4: (a) individual capillary pressure–saturation functions; (b) scaled capillary pressure–saturation functions; and (c) relative permeability functions.

permeability and residual water content were always positive. Final parameter values are listed in Table 4. The final parameter set was taken as that having the minimum RMS value, or was computed from mean parameter values if the optimized parameter sets had similar minimum RMS values. The uncertainty in the optimized parameter values in Table 4 is expressed by the normalized standard deviate (NSD), which was calculated from (j i/b i) ⫻ 100%, where j i is the standard deviation of parameter b i as estimated from the parameter covariance matrix12. Using this type of expression, the uncertainties of the parameters can be evaluated by direct comparison of NSD values. Larger NSD values indicate increasing parameter uncertainty and decreasing parameter sensitivity. In general, it can be concluded from the results in Table 4 that k has the largest uncertainty and smallest sensitivity. Consequently, although

k is a function of the porous medium only, optimized k values vary widely between fluid pair and hydraulic function type. Finally we point out that not all optimizations converged for the BRB model. In addition to visual comparison of predicted and measured cumulative outflow and capillary pressure head values, and evaluation of RMS and NSD values, a third approach to discriminate between parametric models is by comparison of the respective AIC and BIC values, accepting those models with minimum AIC or BIC values. These values are listed in Table 5 for the VGM, LNM, BRB, and GDM models. Although the GDM model has one parameter less than all other models, AIC and BIC values are still largest because of its high RMS values (Table 3). The advantage of the reduced number of fitting parameters of the GDM model is partly offset by the large number of data

Pressure–saturation and permeability functions

491

T

C

P

a

Fig. 5. Optimized constitutive functions (VGM) of Columbia soil corresponding to the parameters listed in Table 4: (a) individual capillary pressure–saturation functions; (b) scaled capillary pressure–saturation functions; and (c) relative permeability functions.

values. Hence we conclude a slight advantage for the VGM, LNM and BRB models. Finally, as a means of testing the modified Leverett scaling relationships using the optimized parameters of the VGM model, we included the scaled capillary pressure functions in Fig. 4 (Lincoln) and Fig. 5 (Columbia), together with the optimized capillary pressure–saturation and permeability curves. The capillary pressure–saturation curve of the air–water system was used in applying the scaling Table 6. Interfacial tension and a vg ratios for oil–water and air–oil systems relative to air–water system

Columbia Lincoln

j ow/j aw

a aw/a ow

j ao/j aw

a aw/a ao

0.38 0.53

0.42 0.51

0.35 0.38

0.40 0.37

relationships of eqn (34). The scaled curves match extremely well at high saturations. Deviations at low saturations are similar to that observed by Demond and Roberts43, who suggested that such deviations were caused by limitations of the surface tension scaling approach associated with uncertainties of the contact angle value. Similar results comparing scaled relationships and optimized relationships were obtained for the remaining constitutive models. For the modified Leverett scaling theory to apply, the a vg values (VGM model) of a fluid pair must be inversely proportional to the corresponding interfacial tension values. Table 6 summarizes the a vg and j ratios of the respective fluid pairs, and indicates that both ratios are almost identical, thereby further validating the parameter estimation approach for two-fluid outflow experiments.

492

J. Chen et al.

4 CONCLUSIONS The soil capillary pressure–saturation and permeability– saturation functions are important for subsurface flow modeling applications. We demonstrate the feasibility of using the inverse parameter estimation method for twofluid flow systems in multi-step outflow experiments. Though it is impossible to remove all inherent uncertainty stemming from the nonlinearity and complexity of fluid transport in soil systems, a posteriori analysis of the inverse problem indicated that the inverse parameter estimation of the constitutive functions from transient multi-step outflow experimental data for two-fluid systems is a well-posed problem. The good agreement between the experimentally measured and optimally simulated system responses indicates that the inversely estimated constitutive relationships do indeed characterize the two porous systems investigated. Model discrimination for inverse parameter estimation is evaluated based on RMS and AIC and BIC values, parameter uniqueness, model convergence, and comparison of the scaled relationships of the optimized functions. Overall, both the commonly applied VGM and the physically-based LNM models show high estimation accuracy and low RMS values, while yielding similar AIC and BIC values. Although the BRB model yielded low RMS values, it showed convergence problems for the oil–water systems, similar to the VGB model. Despite the GDM model having only three fitting parameters, AIC and BIC values were larger than any of the three other models, whereas the resulting RMS values were only slightly larger than those for the three other acceptable models. It appears that the VGM, LNM, BRM, and GDM models are equally applicable, with slight individual nuances, in the inverse parameter estimation of multi-step outflow data. The BC-type models showed a poor performance, likely caused by the absence of measured non-wetting-fluid entry values in the objective function. Further evaluation of the optimum model may depend on the considered soils and fluids, and experimental conditions. A relatively broad range of selected model options also indicates the flexibility and robustness of the presented inverse parameter estimation approach. Advantages inherent to the parameter estimation method include that the measurements are simple and accurate, the simultaneous estimation of capillary pressure and permeability functions for the same soil sample, and the consistency of the application in computer modeling. In summary, we conclude that the multi-step outflow method is an accurate and efficient method for determination of two-fluid flow capillary pressure–saturation and permeability functions.

ACKNOWLEDGEMENTS This study was financed by the EPA (RFA 93-G4) under cooperative agreement No. CR-822204. We are especially grateful to Dr John L. Nieber at the University of Minnesota,

for supplying us with his transient one-dimensional twophase flow code that we adapted for our study.

T

REFERENCES

C

1. Dane, J.H., Oostrom, M. and Missildine, B.C. An improved method for the determination of capillary pressure–saturation curves involving TCE, water and air. J. Contam. Hydro., 1992, 11, 69–81. 2. Demond, A.H. and Roberts, P.V. Estimation of two-phase relative permeability relationships for organic liquid contaminants. Water Resour. Research, 1993, 29, 1081–1090. 3. Klute, A. and Dirksen, C., Conductivities and diffusivities of unsaturated soils. In Methods of Soil Analysis, Agron. Monogr. 9, Part 1, 2nd edn, ed. A. Klute. ASA and SSSA, Madison WI 53711, 1986, pp. 687–734. 4. Burdine, N.T. Relative permeability calculations from pore size distribution data. Petroleum Trans. Am. Inst. Mining Eng., 1953, 198, 71–77. 5. Mualem, Y. and Dagan, G. Hydraulic conductivity of soils: Unified approach to the statistical models. Soil Sci. Soc. Amer. J., 1978, 42, 392–395. 6. Leverett, M.C. Capillary behavior in porous solids. Trans. AIME, 1941, 142, 152–169. 7. Miller, E.E. and Miller, R.D. Physical theory for capillary flow phenomena. Journal of Applied Physics, 1956, 27, 324–332. 8. Dane, J.H. and Hruska, D. In-situ determination of soil hydraulic properties during drainage. Soil Sci. Soc. Am. J., 1983, 47, 619–624. 9. Eching, S.O. and Hopmans, J.W. Optimization of hydraulic functions from transient outflow and soil water pressure data. Soil Sci. Soc. Am. J., 1993, 57, 1167–1175. 10. Eching, S.O., Hopmans, J.W. and Wendroth, O. Unsaturated hydraulic conductivity from transient multistep outflow and soil water pressure data. Soil Sci. Soc. Am. J., 1994, 58, 687–695. 11. Hornung, U., Identification of nonlinear soil physical parameters from an output–input experiment, in numerical treatment of inverse problems in differential and integral equations, ed. P. Deuflhard and E. Hairer. Birkhauser, Boston MA, 1983, pp. 227–237. 12. Kool, J.B. and Parker, J.C. Analysis of the inverse problem for transient unsaturated flow. Water Resources Research, 1988, 24, 817–830. 13. Parker, J.C., Kool, J.B. and Van Genuchten, M.Th. Determining soil hydraulic properties from onestep outflow experiments by parameter estimation: II. Experimental studies. Soil Sci. Soc. Amer. J., 1985, 49, 1354–1359. 14. Russo, D., Bresler, E., Shani, U. and Parker, J.C. Analysis of infiltration events in relation to determining soil hydraulic properties by inverse problem methodology. Water Resources Research, 1991, 27, 1361–1373. 15. Toorman, A.F., Wierenga, P.J. and Hills, R.G. Parameter estimation of hydraulic properties from one-step outflow data. Water Resources Research, 1992, 28, 3021–3028. 16. 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. 17. Zachmann, D.W., Duchateau, P.C. and Klute, A. Simultaneous approximation of water capacity and soil hydraulic conductivity by parameter identification. Soil Sci., 1982, 134, 157–163. 18. Carrera, J. and Neuman, S.P. Estimation of aquifer parameters under transient and steady state conditions: 2. Uniqueness, stability, and solution algorithms. Water Resources Research, 1986, 22, 211–227.

P

a

Pressure–saturation and permeability functions 19. Gardner, W.R. Calculation of capillary conductivity from pressure plate outflow data. Soil Sci. Soc. Amer. Proc., 1956, 20, 317–320. 20. van Dam, J.C., Stricker, J.N.M. and Droogers, P. Inverse method to determine soil hydraulic functions from multistep outflow experiments. Soil Sci. Soc., Am. J., 1994, 58, 647–652. 21. Russo, D. Determining soil hydraulic properties by parameter estimation: on the selection of a model for the hydraulic properties. Water Resources Research, 1988, 24, 453–459. 22. Liu, Y.P., Hopmans, J.W., Grismer, M.E. and Chen, J.Y. Direct estimation of air–oil and oil–water capillary pressure and permeability relations from multi-step outflow experiments. J. Contam. Hydrology, 1998, 32, 21–42. 23. Schroth, M.H., Istok, J.D., Ahearn, S.J. and Selker, J.S. Geometry and position of light nonaqueous phase liquid lenses in water-wetted porous media. J. Contam. Hydrol., 1995, 19, 269–287. 24. Hopmans, J.W., Vogel, T. and Koblik, P.D. X-ray tomography of soil water distribution in one-step outflow experiments. Soil Sci. Soc. Am. J., 1992, 56, 355–362. 25. Whitaker, S.A. The closure problem for two-phase flow in homogeneous porous media. Chemical Engineering Science, 1994, 49, 765–780. 26. Bear, J., Dynamics of fluids in porous media. Dover Publications, New York, 1972. 27. Touma, J. and Vauclin, M. Experimental and numerical analysis of two-phase infiltration in a partially saturated soil. Tranport in Porous Media, 1986, 1, 27–55. 28. Celia, M.A. and Binning, P. A mass conservative numerical solution for two-phase flow in porous media with application to unsaturated flow. Water Resources Research, 1992, 28, 2819–2828. 29. 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. 30. Brooks, R.H. and Corey, A.T., Hydraulic properties of porous media. Colorado Univ., Hydrology paper No. 3, 1964. 31. Kosugi, K. Three-parameter lognormal distribution model for soil water retention. Water Resour. Research, 1994, 30, 891–901. 32. Kosugi, K. Lognormal distribution model for unsaturated soil hydraulic properties. Water Resources Research, 1996, 32, 2697–2703.

493

33. Brutsaert, W. Some methods of calculating unsaturated permeability. Trans. ASAE, 1967, 10, 400–404. 34. Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 1976, 12, 513–522. 35. 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 Resources Research, 1989, 25, 2187–2193. 36. van Genuchten, M.Th., Calculating the unsaturated hydraulic conductivity with a new closed-form analytical model, Res. Rep. 78-WR-08, Dept. of Civil Eng., Princeton Univ., Princeton NJ, 1978. 37. Gardner, W.R. Some steady state solutions of unsaturated moisture flow equations with application to evaporation from water table. Soil Sci., 1958, 85, 228–232. 38. Finsterle, S. and Pruess, K. Solving the estimation-identification problem in two phase flow modeling. Water Resources Research, 1995, 31, 913–924. 39. Celia, M.A., Bouloutas, E.T. and Zarba, R.L. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resources Research, 1990, 26, 1483–1496. 40. Carrera, J. and Neuman, S.P. Estimation of aquifer parameters under transient and steady state conditions: 1. Maximum likelihood incorporating prior information. Water Resources Research, 1986, 22, 199–210. 41. More, J.J., The Levenberg–Marquardt algorithm: Implementation and theory. In Lecture Notes in Mathematics, 630, ed. G.A. Watson. Springer-Verlag, New York, 1977. 42. Kool, J.B., Parker, J.C., and van Genuchten, M. Th., ONESTEP: A nonlinear parameter estimation program for evaluating soil hydraulic properties from one-step outflow experiments. Virginia Agricultural Experiment Station, Bulletin 85-3, 1985. 43. Demond, A.H. and Roberts, P.V. Effect of interfacial forces on two-phase capillary pressure–saturation relationships. Water Resour. Research, 1991, 27, 423–437. 44. Simunek, J. and van Genuchten, M.Th. Parameter estimation of soil hydraulic properties from the tension disc infiltrometer experiment by numerical inversion. Water Resources Research, 1996, 32, 2683–2696. 45. Bard, Y., Nonlinear Parameter Estimation. Academic Press, Orlando FL, 1974.

T

C

P

a