Published online May 16, 2007

SOIL PHYSICS

Uniqueness of Soil Hydraulic Parameters
Determined by a Combined Wooding
Inverse Approach
N. Lazarovitch*
The Wyler Dep. of Dryland Agriculture
Jacob Blaustein Inst. for Desert Research
Ben-Gurion Univ. of the Negev
Sede Boqer
84990 Israel

A. Ben-Gal
Inst. of Soil, Water and Environmental
Sciences
Agricultural Research Organization
Gilat Research Center
D.N. Negev
85280 Israel

J. Šimůnek
Dep. of Environmental Sciences
Univ. of California–Riverside
Riverside, CA 92521

U. Shani

Knowledge of soil hydraulic properties is essential for a proper understanding and evaluation of physical and chemical processes within the vadose zone involved in variably saturated
water flow and transport of water-dissolved salts and pollutants. Soil hydraulic properties areoften expressed using functional relationships between the soil hydraulic conductivity (K),water content (θ), and matric potential Ψ). Our objectives were (i) to combine the Wooding’s
analytical solution for steady-state infiltration from a circular pond and the inverse determination of parameters from transient infiltration events into a coupled method for the in situ
estimation of soil hydraulic properties, and (ii) to develop a simple semiautomatic device forthe in situ estimation of the soil hydraulic functions. The experimental method consists ofmeasurements of transient (short-term) and steady-state (long-term) infiltration flow rates
from a set of rings having different radii, each positioned sequentially at the same locationon the soil surface. A shallow water depth is maintained within a 1- to 2-mm range over thesoil surface with an electrode set. The flow rate is determined by continuous weighing of awater reservoir. The flow is monitored and controlled by a laptop computer, which also automatically calculates the soil hydraulic properties from collected data. The coupled method
starts with the application of the Wooding’s analytical solution to obtain estimates of thesoil hydraulic properties using steady-state fluxes. These estimates are then finalized using
numerical inversion of the transient data. The coupled method was evaluated using numerically generated data. Unique and fast reproduction of soil hydraulic properties for generated
data was obtained. The method’s applicability was also tested using field experiments for two
soils. The technique was found to be sound and the device simple to operate.

Faculty of Agricultural, Food andEnvironmental Quality Sciences
The Hebrew Univ. of Jerusalem
Rehovot, IsraelKnowledge of soil hydraulic properties is essential for proper
understanding and evaluating of the physical and chemical processes involved in transport of water, dissolved salts, andpollutants within the vadose zone (Si et al., 1999; Al-Jabri et al.,2002). Soil hydraulic functions describing hydraulic properties
are highly nonlinear, and thus their laboratory or field measurements are tedious and time consuming, and involve considerable uncertainty (Shani and Or, 1995; Hopmans et al., 2002;
Bodhinayake et al., 2004). There are many in situ methods for
direct estimation of the soil hydraulic functions (Hopmans etal., 2002). Examples include the crust method, where steady fluxis imposed at the soil surface (Hillel and Gardner, 1970), theSoil Sci. Soc. Am. J. 71:860–865
doi:10.2136/sssaj2005.0420
Received 26 Dec. 2005.
*Corresponding author (lazarovi@bgu.ac.il).
© Soil Science Society of America
677 S. Segoe Rd. Madison WI 53711 USA
All rights reserved. No part of this periodical may be reproduced ortransmitted in any form or by any means, electronic or mechanical,
including photocopying, recording, or any information storage andretrieval system, without permission in writing from the publisher.
Permission for printing and for reprinting the material contained
herein has been obtained by the publisher.

860

instantaneous profile method (Watson, 1966), and unit gradientinternal drainage (Libardi et al.,1980). Although relatively simple in concept, such direct methods have a number of limitations
that significantly limit their use in practice. The main limitationis that these methods tend to be very time consuming due to theneed to adhere to relatively strict initial and boundary conditions(Šimůnek and van Genuchten, 1996). Another limitation, especially in soils with low matric heads, is that steady-state conditions are reached only after a significant amount of time due toflow dependence on hydraulic conductivity (Zhang, 1998).
Recently, parameter estimation methods have been applied todetermine soil hydraulic properties (Russo et al., 1991; Hopmans
et al., 2002). Inherent to parameter estimation is the assumptionthat soil hydraulic properties may be described by a relatively simple deterministic model containing a small number of unknown
parameters (Russo et al., 1991; Shani and Gordin-Katz, 1998).
Contrary to the previously mentioned direct methods, parameter
estimation with a predefined hydraulic model does not impose
any strict requirements on the initial and boundary conditions of
the measured system. Experimental methods based on the parameter estimation approach usually require less time and labor thandirect methods and thus enable a larger number of measurements.
This is especially important for the in situ characterization of largeor heterogeneous sites (Russo et al., 1991).
SSSAJ: Volume 71: Number 3 • May–June 2007

Two main approaches are often used for in situ estimation of hydraulic parameters of near-surface soils via infiltration
experiments (Shani and Gordin-Katz, 1998). The first approach
involves inverse determination of parameters from transient
infiltration events (Kool et al., 1985; Parker et al., 1985; Russo
et al., 1991; van Dam et al., 1992; Warrick, 1993; Šimůnek andvan Genuchten 1996, 1997; Schwartz and Evett, 2002). The
second approach applies Wooding’s approximate analytical solution (Wooding, 1968) for steady-state infiltration from a circularpond with a constant pressure head at the soil surface (White and
Sully, 1987; Shani et al., 1987; Smettem and Clothier, 1989).
In Wooding’s-based approaches, the matrix flux potential, φ[L2 T−1], and the saturated hydraulic conductivity, Ks [L T−1],
are obtained either from multiple steady-state infiltration eventsusing disks of several sizes (Shani et al., 1987; Smettem andClothier, 1989), or from a series of infiltration events from asingle disk with multiple tensions (Ankeny et al., 1991). The
Wooding’s-based approach is usually limited to two-parameter hydraulic conductivity models, such as Gardner’s (1958).
Additional information is required when more complex hydraulic models are desired, i.e., models that use more parameters tobetter describe the hydraulic relationships.
Two hydraulic parameters can also be accurately estimatedfrom infiltration experiments using the inverse method. For models using larger number of hydraulic parameters, however, thisapproach is limited due to problems associated with convergenceof parameter estimation method, nonuniqueness of optimized
parameters, and computational efficiency of their overall optimization scheme, especially when several hydraulic parameters mustbe estimated simultaneously (Russo et al., 1991; van Genuchten
and Leij, 1992; Šimůnek and van Genuchten, 1996). Specifically,Šimůnek and van Genuchten (1996) used a quasi-three-dimensional (axisymetrical two-dimensional) numerical code for identifying soil hydraulic parameters from unsaturated flow data froma tension disk infiltrometer. They showed that instantaneous or
cumulative infiltration rates alone do not provide a unique solution while solving for more than two parameters. A similar conclusion was reached by Russo et al. (1991) for ponded infiltration.
In this study, we developed a coupled parameter estimation
method based on the numerical solution of the Richards equation
with Wooding’s analytical solution for flow from a circular surface
source, to better use available information for the determination ofmultiple-parameter hydraulic models. Infiltration data from a surface disk infiltrometer with several radii are used to simultaneously
obtain both the transient information needed for the numericalinverse approach and the steady-state infiltration rate required
for the Wooding’s analytical solution. Furthermore, we have constructed a simple semiautomatic device for the in situ estimationof the soil hydraulic conductivity and retentivity.

MATERIALS AND METHODS
Determination of Soil Hydraulic Parameters
Inverse Numerical Solution
Inverse determination of soil hydraulic parameters from transient infiltration from a surface disk source involves a numerical solution of the flow
equation augmented with parameterized soil hydraulic functions and suitable initial and boundary conditions. Initial estimates of optimized parameters are iteratively improved during the optimization process that minimizes

the sum of squared deviations between measured and simulated variables
until a desired degree of precision is obtained (Šimůnek et al., 1998).
The governing equation for radial symmetric isothermal Darcian
flow in a variably saturated isotropic rigid porous medium is given by
the following modified form of the Richards (Richards, 1931) equation (e.g., Warrick, 1992):∂θ 1 ∂ ⎜⎛ ∂ψ ⎞⎟ ∂ ⎛⎜ ∂ψ ⎞⎟ ∂K⎟⎟ +=⎟ + ⎜K⎜rK∂t r ∂r ⎝⎜ ∂r ⎠⎟ ∂z ⎜⎜⎝ ∂z ⎠⎟ ∂z[1]

where r is a radial coordinate [L], z is the vertical coordinate [L] positive upward, K is the hydraulic conductivity [L T−1], and t is time [T].
Equation [1] is solved numerically for the following initial and boundary
conditions that are appropriate for infiltration from a circular source:ψ (r , z , t ) = ψ it =0ψ (r , z , t ) = ψ 0 = 0−∂ψ (r , z , t )∂z[2]
0 < r < r0 ,= 1 r > r0 , z = 0

ψ (r , z , t ) = ψ i

r 2 +z 2 → ∞

z =0

[3]
[4]
[5]

where ψi is the initial soil pressure head [L], and ψ0 is the inlet pressure headat the soil surface in the infiltration ring [L], which is equal to zero for our
specific conditions, i.e., water flow from a soil surface disk source.

Hydraulic Model
The van Genuchten–Mualem hydraulic model (Mualem, 1976;
van Genuchten, 1980) was chosen to demonstrate our coupledapproach. This model was selected because, although it is widely used,
there is no simple in situ method to estimate its parameters. The functional form of the van Genuchten–Mualem hydraulic model isSe =θ − θrn −1
n −m= ⎡⎢1 + (−αψ) ⎤⎥ ; m =⎦nθs − θ r ⎣m ⎤2⎡
K (ψ) = K s Se 0.5 ⎢1 − 1 − Se1/ m ⎥
⎣⎢
⎦⎥

(

)

[6]
[7]

where Se is the effective fluid saturation [dimensionless], θr and θsdenote the residual and saturated water contents [L3 L−3], respectively; and α [L−1], n [dimensionless], and m (= 1 − 1/n) [dimensionless] are empirical shape parameters.

Estimation of Saturated Hydraulic Conductivity andMatric Flux Potential using Wooding’s Approach
Wooding’s (1968) analytical solution of Eq. [1–5] is represented
in the following equation that relates the steady-state soil water flux,
q∞ [L T−1], from a shallow disk to the soil hydraulic properties:
q∞ =

Q∞
πr0

2

= Ks +

4φ
πr0

[8]

where Q∞ [L3 T−1] is the source discharge for steady-state conditions,r0 [L] is the infiltration disk radius, and φ is the matrix flux potential[L2 T−1] (Gardner, 1958):0

φ = ∫ K (ψ)dψ

[9]

ψi

SSSAJ: Volume 71: Number 3 • May–June 2007861

Equation [8] is linear and thus estimates of Ks and φ can be determined from a linear regression of q∞ vs. r0−1. The value of Ks can be
estimated from the intercept of the linear regression, and φ can beobtained as

φ=

bπ
4

[10]

where b is the slope of the fitted curve [L2 T−1]. Values for q∞(r0−1)
can be obtained from multiple infiltration measurements using varied disk diameters (Shani et al., 1987; Smettem and Clothier, 1989;
Thony et al., 1991; Al-Jabri et al., 2002).

Coupled Approach for Estimation of Soil
Hydraulic Parameters
Combining the two approaches (e.g., numerical parameter estimation and Wooding’s analytical solution) allows a more robust andmore unique estimation of optimized parameters. In this case, the Ks
and φ values, obtained using the Wooding approach, serve to reduce
the number of inversely sought parameters and to constrain optimized
parameters. The objective function, Φ, for the numerical approach isK2Φ (β) = ∑ Wi ⎡⎢⎣ I i * ( t i ) − I i (β, t i )⎤⎥⎦

[11]

i =1

where Wi is the weight of a particular measured point, Ii* [L] is themeasured cumulative infiltration at time ti, and Ii [L] stands for thesimulated cumulative infiltration obtained using parameter set β:I (t ) =

t

t

0

0

J
1
1
d
q
t
t
=
∫ ()
∫ ∑ i j (t )dt
πr02 tπr02 t j =1

[12]

where t0 is the starting time of infiltration [T], q(t) is the instantaneous ring infiltration rate [L3 T−1], J [dimensionless] is the numberof nodes representing the infiltration ring in the numerical scheme,
and ij is the infiltration rate [L3 T−1] in node j.

Numerical Simulation of Infiltration
The HYDRUS-2D software package (Šimůnek et al., 1999) wasused for the numerical solution of the Richards equation. A transportdomain was selected such that the outer boundaries did not affect theflow field inside of the domain. The transport domain (50 by 50 cm)
was discretized into unstructured finite element mesh with 3505 tri-

angular finite elements with significantly smaller triangles immediately surrounding the infiltration ring and increasingly larger triangles
as distance from the ring increased. Boundary conditions followed Eq.
[3–4], whereas the initial conditions were specified in terms of soil
water content and were slightly higher than the residual values.
An optimization procedure implemented in the HYDRUS-2D code
was used to match the modeled cumulative infiltration with measured
data. Detailed description of the basic procedure for parameter optimization has been given by Šimůnek and van Genuchten (1996). The hydraulic model described by Eq. [6] and [7] includes five parameters: θs, θr, Ks,
α, and n. Of these parameters, the saturated water content (θs) needs to
be independently measured, while the residual water content (θr) needs
to be independently estimated, possibly by using pedotransfer functions(Schaap et al., 2001). The added information from the Wooding analysis
allows estimation of Ks and φ, which in turn relates parameters α with n
and is derived from Eq. [6], [7], and [9] as follows:
0nφ = K s ∫ ⎢⎡1 + (−αψ ) ⎥⎤⎣⎦−∞n −1 ⎞2⎧⎪ ⎡⎫⎪ n ⎟⎟⎟n ⎤ −1 ⎪⎪1 − ⎢1 + (−αψ) ⎥ ⎬ ⎟⎟ dψ⎜1 − ⎨⎦ ⎭⎪⎪ ⎣⎪⎪ ⎟⎟⎟⎜⎜⎜ ⎩⎝⎠0.5(1−n ) ⎛
⎜⎜
n
⎜

[13]

Following the estimation of Ks and φ using Eq. [6] from multiple
radii infiltration runs, the inverse optimization includes the following
steps: (i) fixing Ks to that obtained by the Wooding analysis; (ii) selection of an initial value of n; (iii) evaluation of α from Eq. [13] by root
finding, i.e., using a nonlinear root-finding procedure that involves abracketing and bisection method (Press et al., 1986); (iv) calculationof I(t) using the HYDRUS-2D code; (v) evaluation of the objectivefunction Eq. [11]; and (vi) updating the n value using the Levenberg–
Marquardt nonlinear minimization procedure (Marquardt, 1963).
This iterative procedure (Steps iii–vi) continues until required convergence criteria are satisfied (Marquardt, 1963).
To evaluate robustness, solution uniqueness, and computational
efficiency of the coupled parameter estimation method, direct numerical simulation of a 0.5-h infiltration event from a soil surface disk source(r0 = 0.0185 m) was used to generate infiltration data for two soils withknown hydraulic properties. The effect of the selection of initial parameter values was evaluated by assuming two different sets of initial estimatesand by studying the rate and accuracy of the convergence procedure to the
known values. The initial estimates included parameters that yielded both
higher and lower hydraulic conductivities than the true values. Hydraulicproperties were selected from the HYDRUS-2D database with θr, θs, α, n,and Ks of 0.078, 0.43, 3.6 m−1, 1.56, and 2.88 × 10−6 m s−1 for a loamy
soil and 0.057, 0.41, 12.4 m−1, 2.28, and 4.05 × 10−5 m s−1 for a loamy
sand soil. Initial estimated values of n were 1.1 and 2.1 for the loamy soil
and 2.1 and 3.1 for the loamy sand soil. Initial estimates of α values were
calculated from Eq. [13] using the selected n values. Calculated α values
were 0.38 and 6.24 m−1 for the loamy soil and 9.8 and 1.56 m−1 for the
loamy sand soil.

Device Design

Fig. 1. Schematic of a device for measuring infiltration rates into soil fromrings: A, container; B, electric valve; C, ring; D, electrode unit; E,
suspended scale; F, laptop computer; G, controller.

862

Implementation of the coupled method relies on in situ measurementof infiltration via rings of varied sizes. An illustration of a semiautomatic
device for conducting these measurements is depicted in Fig. 1. Water
flows from a container (A) via a low-pressure, ½″, two-way electric solenoid(ACL, Caponago, Italy) (B) to the ring (C) positioned at the soil surface.
Depth of the water within the ring is maintained by an electrode unit (D)that controls the electric valve. One electrode is attached vertically in the
middle of the ring and the height of its tip determines the water level above
the soil. The second electrode is grounded through the ring wall. The solenoid (on the water supply tube) (B) is open when there is no electric curSSSAJ: Volume 71: Number 3 • May–June 2007

rent between the two electrodes, i.e., water is not shorting them. The solenoid is switched to the closed position when water connects the electrodes,
enabling an electric current. Initially, there is no water within the ring, thesolenoid is open, and water flows inside via the supply tube. When waterreaches the electrode tip, the electric circuit between the electrodes is closed,
the solenoid closes, and water supply to the cylinder stops. As the waterlevel decreases (due to infiltration), the electrode inside the ring is exposed
and water supply resumes. Using this method and needle electrodes, which
decrease the surface tension between the water and the electrode, water headfluctuations can be reduced to 1 mm. The energy requirement for this unit
is mostly for the solenoid and any commercial 12-V battery-operated solenoid will suffice. Flow rate is measured by monitoring the weight change of
the water supply using a suspended s-type load cell (XLS2-HSS, Load Cell
Central, Monroeton, PA) (E). When the flow reaches steady state, measurements are terminated. A laptop computer (F) that is attached to the scale via
a controller (MicroLogix 1000, Allen-Bradley, Milwaukee, WI) (G) records
weights and calculates the flow rate in real time. The same computer canalso analyze the data and determine the soil hydraulic properties in situ.

Field Tests

Table 1. Textural properties and the residual (θr) and saturated (θs) volumetric water contents of the studied soils.
Site
Hazerim
Rehovot

Clay

θr

θs

———— % ————
75
12.5
12.5
95
4
1

0.05
0.03

0.48
0.4

Sand

Silt

WIA procedure except that it did not use the matrix flux potential constraint (obtained from the Wooding analytical solution)
on optimized parameters α and n. A convergence to n = 1.56,
the value used in the direct solution, was attained with the WIA
code regardless of the initial estimate of n for the loamy soil. A
poor convergence was achieved with the SIA, with neither casereaching the true value of n. For the loamy sand soil, the WIA
again produced a precise convergence (n = 2.279) for both initial estimates, while the SIA produced a successful convergence
only when smaller value of n was used as the initial estimate.
The SIA diverged when a higher initial estimate of n was used.
The calculation time with the WIA was approximately fourtimes shorter than with the SIA, and the number of iterations
needed to converge to the true value of optimized parameters
was smaller. Convergence schemes for higher and lower initial
estimates of n using both WIA and SIA are given for the loamy
soil in Fig. 2 and for the loamy sand soil in Fig. 3. Dashed lines
represent the assembly of all possible log Φ(n) values that obey
the φ value within the presented n range. Note that any n yields asingle value of α (Eq. [13]) due to the constraint imposed by thematrix flux potential φ. Data points represent the actual convergence paths of particular runs. With the WIA, all convergencepaths fall directly on the dashed lines. Depending on the optimization algorithm, the convergence path, while always descending
toward the minimum, can shift from one side to the other of
the true value of n while converging. Since the φ constraint on
optimized parameters is not used in the SIA code, log Φ(n) cantake any value within the domain, and hence the convergence is
uncertain and tedious. The φ constraint decreases the number of
optimized parameters in the WIA method compared with the
SIA method by directly linking parameters α and n.

Infiltration trials were made with two soils: Hazerim sandy loam
and Rehovot sand. Particle size distribution and saturated and residual
water contents are summarized in Table 1. The soil surface was clearedof large clods and stones and smoothed using a brush before placementof the ring and the electrode unit.
Three rings of different radii were used for each soil: 0.0185, 0.025,
and 0.05 m for the Rehovot sand and 0.0185, 0.025, and 0.075 m forHazerim sandy loam. Shallow water depth was maintained within 1 to
2 mm above the soil surface by the electrode unit. Each measurement
set started with the smallest ring placed on dry soil and, after steadystate flow was reached, progressed successively to the next larger ring.
Steady state was assumed when discharge changes were <0.5% during
an interval of 5 min. Only data from the first ring were used for thetransient analysis since they provided the greatest time rate of change ininfiltration rate and therefore maximum information for the parameter
estimation method. A single run was chosen as representative and used
for hydraulic property determination. The average initial water content
was measured at several locations and depths (down to 30 cm) using
time domain reflectometry (TDR-100, Campbell Scientific, Logan,
UT). The residual water content was determinedin the laboratory by oven drying undisturbed soil Table 2. Comparison between inverse procedures with numerically generated
data† for two soils.
cores of known volume and by estimation of bulk
density. The average initial water content was foundInitialFinal
Calculation
Case
Soil
Iterations
to be slightly higher than the residual water contimennααtent. The saturated water content was estimated
s
m−1
m−1
using the TDR measurement from the saturated
Wooding inverse approach with Ψi = –10 m, θr = 0.078, θs = 0.43, Ks = 2.88 × 10−6 m s−1
1
Loam
0.38
1.1
3.59
1.56
1037
5
zone inside the largest ring taken after the last infil2
Loam
6.24
2.1
3.59
1.559
1112
4
tration experiment.

RESULTS AND DISCUSSION
Numerical Tests
Initial estimates, predicted hydraulicparameters vs. values that were used in the
direct computation, computation time, andthe number of iterations for the coupled
Wooding inverse approach (WIA) and thestandard inverse approach (SIA) are summarized in Table 2 for two soils (loam and loamysand). The SIA approach was similar to the
SSSAJ: Volume 71: Number 3 • May–June 2007

Standard inverse approach with Ψi = –10 m, θr = 0.078, θs = 0.43, Ks = 2.88 × 10−6 m s−1
3
Loam
0.38
1.1
3.23
1.508
4203
14
4
Loam
6.24
2.1
5.36
1.865
5638
12
Wooding inverse approach with Ψi = –5 m, θr = 0.057, θs = 0.41, Ks = 4.05 × 10−5 m s−1
5
Loamy sand
9.8
1.9
12.4
2.279
1021
3
6
Loamy sand
1.56
3.1
12.4
2.279
1704
4
Standard inverse approach with Ψi = –5 m, θr = 0.057, θs = 0.41, Ks = 4.05 × 10−5 m s−1
7
Loamy sand
9.8
1.9
12.4
2.279
5565
14
8
Loamy sand‡ 1.56
3.1
–
–
–
–
†α and n, empirical shape parameters; Ψi, initial soil water pressure head; θr and θs, residual
and saturated volumetric water content, respectively; Ks, saturated hydraulic conductivity.
‡ Unsuccessful run.

863

Fig. 2. Log of the objective function (Φ) as a function of the empiricalshape parameter n for the hypothetical case of the loamy soil.
Data points connected by lines represent optimization Cases 1
through 4 from Table 2. The dashed line represents the assemblyof all possible log Φ that obey the matric flux potential (φ) value
within the presented n range.

Field Tests
Experimentally measured, steady-state soil water fluxes, q∞,
as a function of the reciprocal ring radii, r0−1, for two soils, using
the prescribed methodology are depicted in Fig. 4. Linearity, as
projected in Eq. [6], is clearly evident. Increasing the size of theinfiltration ring resulted in decreasing the steady-state flux. It canbe seen from Fig. 4 that both Ks (the intercept) and φ (proportional
to the slope) decreased for the finer textured Hazerim sandy loam
compared with the coarser textured Rehovot sand. From the linear
regression of the plot, Ks of the Rehovot sand is 1.33 × 10−4 m s−1
and φ is 6.78 × 10−6 m2 s−1, and for the Hazerim sandy loam Ks is5.22 × 10−6 m s−1 and φ is 4.04 × 10−6 m2 s−1. The coefficient ofdetermination (r 2) for the linear regression was 1 for the Rehovot
sand and 0.9994 for the Hazerim sandy loam. These high valuesof r 2 demonstrate the robustness of Wooding’s solution. Although
two rings alone are sufficient to satisfy the analytical solution, we
recommend using three rings to ensure consistency.
Simulated and measured cumulative infiltration, I, of representative case for Hazerim sandy loam as a function of time withvarious r0 is presented in Fig. 5. The transient stage of the infiltra-

Fig. 3. Log of the objective function (Φ) as a function of the empirical shape parameter n for the hypothetical case of the loamy
sand soil. Data points connected by lines represent optimization
Cases 5 through 7 from Table 2. The dashed line represents theassembly of all possible log Φ that obey the matric flux potential(φ) value within the presented n range.
864

Fig. 4. Measured (symbols) steady-state flux (q∞) as a function of
the reciprocal of the ring radius (r0−1) for Rehovot sand andHazerim sandy loam. Solid lines represent the best fit.

tion for the smallest ring radius was used for the numerical inverse
solution and, together with parameters obtained from Fig. 4(steady-state infiltration), parameters of the van Genuchten (1980)
hydraulic model were estimated. Excellent agreement between
the measured and fitted cumulative field infiltration curves was
obtained. Comparison of hydraulic parameters obtained from the
proposed combined method with those from the field dripper
method (Shani et al., 1987) or from retention curve data determined by the least-square optimization technique of the RETC
program (van Genuchten et al., 1991) is summarized in Table 3.

CONCLUSIONSWe propose a method and a device for the in situ estimationof soil hydraulic properties. The method consists of measurementsof transient (short-term) and steady-state (long-term) flow rates
from a set of rings having several different radii, each positioned
at the soil surface. Determination of Ks and φ is accomplishedusing Wooding’s analytical solution for the steady-state flow rates
and then α and n determined by using the numerical inversion of
the transient flow rates. This combination of numerical parameter
estimation and Wooding’s analytical solution allowed robust andunique estimation of optimized parameters for the sand and sandy
loam soils evaluated in field trials. The Ks and φ values, obtainedusing Wooding’s approach, serve to reduce the number of inversely
required parameters and to constrain optimized parameters. By

Fig. 5. Measured (symbols) and fitted (solid line) cumulative infiltration (I) for Hazerim loamy sand as a function of time with various ring radii (r0). Fitted line was calculated with the Wooding–inverse approach.
SSSAJ: Volume 71: Number 3 • May–June 2007

Table 3. Saturated hydraulic conductivity (Ks), empirical shapeparameters α and n, and matrix flux potential (φ) of thesoils under study as measured by different methods.
Method

Ks

α

m s−1

m−1

Hazerim sandy loam
Field dripper
5.22 × 10−6
Retention curve–
Combined method
1.22 × 10−5
Rehovot sandField dripper
1.4 × 10−4
Retention curve–
Combined method
1.33 × 10−4nφm2 s−1

–
1.66
1.59

–
1.644
1.81

4.04 × 10−6–
3.01 × 10−6

–
5.92
10.7

–
2.28
2.66

1.05 × 10−5–
6.78 × 10−6

constraining optimized parameters, more accurate and more complicated hydraulic models involving more parameters can be evaluated. The flux potential that constrains parameters in the inverse
optimization can obviously be obtained using other methods. The
effectiveness of the inverse solution will be increased irrespectiveof how this constrain is obtained. Similar improvement could beobtained if one or more parameters were measured independently
and fixed during the optimization.
The experimental device maintains a shallow water depth
within a 1- to 2-mm range above the soil surface with an electrodeset. The flow rate is determined by continuous weighing of a water
reservoir. The system is mobile and able to be powered from a car
cigarette lighter. Flow is automatically monitored and controlled
by a laptop computer, which also calculates the hydraulic properties. Exact spatial coordinates of the measurement site can be easily
obtained using a global positioning system (GPS). The device is simple to operate and measurements can be taken in remote locations
where the water supply is limited. The method is not destructive,
thus allowing measurement of changes with time within a specificsampling site. Since the boundary water content is saturated, steadystate flux is reached after a short time, compared with methods that
impose tension at the surface (e.g., tension disk infiltrometers). This
allows an increased number of sampling points across a field.
REFERENCES
Al-Jabri, S.A., R. Horton, and D.B. Jaynes. 2002. A point-source method for
rapid simultaneous estimation of soil hydraulic and chemical transport
properties. Soil Sci. Soc. Am. J. 66:12–18.
Ankeny, M.D., M. Ahmed, T.C. Kaspar, and R. Horton. 1991. Simple field
method for determining unsaturated hydraulic conductivity. Soil Sci.
Soc. Am. J. 55:467–470.
Bodhinayake, W., B.C. Si, and K. Noborio. 2004. Determination of
hydraulic properties in sloping landscapes from tension and double-ring
infiltrometers. Vadose Zone J. 3:964–970.
Gardner, W.R. 1958. Some steady-state solutions of unsaturated moisture flow
equations with application to evaporation from a water table. Soil Sci.
85:228–232.
Hillel, D., and W.R. Gardner. 1970. Measurement of unsaturated hydraulicconductivity by infiltration through an impeding layer. Soil Sci.
109:149–153.
Hopmans, J.W., J. Šimůnek, N. Romano, and W. Durner. 2002. Inverse Methods.
p. 963–1008. In J.H. Dane and G.C. Topp (ed.) Methods of soil analysis.
Part 4. Physical methods. 3rd ed. SSSA Book Ser. 5. SSSA, Madison, WI.
Kool, J.B., J.C. Parker, and M.Th. van Genuchten. 1985. Determining soil hydraulicproperties from one-step outflow experiments by parameter estimation: I.
Theory and numerical studies. Soil Sci. Soc. Am. J. 49:1348–1354.
Libardi, P.L., K. Reichardt, D.R. Nilsen, and J.W. Biggar. 1980. Simple field method
for estimating soil hydraulic conductivity. Soil Sci. Soc. Am. J. 44:3–7.
Marquardt, D.W. 1963. An algorithm for least-squares estimation of nonlinear

SSSAJ: Volume 71: Number 3 • May–June 2007

parameters. J. Soc. Ind. Appl. Math. 11:431–441.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity ofunsaturated porous media. Water Resour. Res. 12:513–522.
Parker, J.C., J.B. Kool, and M.Th. van Genuchten. 1985. Determining soil
hydraulic properties from one-step outflow experiments by parameter
estimation: II. Experimental studies. Soil Sci. Soc. Am. J. 49:1354–1359.
Press, W.S., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. 1986.
Numerical recipes in FORTRAN 77: The art of scientific computing.
Cambridge Univ. Press, Cambridge, UK.
Richards, L.A. 1931. Capillary conduction of liquids through porous mediums.
Physics 1:318–333.
Russo, D., E. Bresler, U. Shani, and J.C. Parker. 1991. Analyses of infiltration
events in relation to determining soil hydraulic properties by inverse
problem methodology. Water Resour. Res. 27:1361–1373.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA:
A computer program for estimating soil hydraulic parameters withhierarchical pedotransfer functions. J. Hydrol. 251:163–176.
Schwartz, R.C., and S.R. Evett. 2002. Estimating hydraulic properties of a finetextured soil using a disc infiltrometer. Soil Sci. Soc. Am. J. 66:1409–1423.
Shani, U., and R. Gordin-Katz. 1998. Soil hydraulic models selection based on
in-situ measurements. J. Irrig. Drain. Eng. 124:285–289.
Shani, U., R.J. Hanks, E. Bresler, and C.A.S. Oliveiera. 1987. A simple field
method for estimating the hydraulic conductivity and matric potential–
water content relations of soils. Soil Sci. Soc. Am. J. 51:298–302.
Shani, U., and D. Or. 1995. In situ method for estimating subsurface
unsaturated hydraulic conductivity. Water Resour. Res. 21:1863–1870.
Si, B., R.G. Kachanoski, F. Zhang, G.W. Parkin, and D.E. Elrick. 1999.
Measurement of hydraulic properties during constant flux infiltration:
Field average. Soil Sci. Soc. Am. J. 63:793–799.
Šimůnek, J., R. Angulo-Jaramillo, M.G. Schaap, J.-P. Vandervaere, and M.Th. van
Genuchten. 1998. Using an inverse method to estimate the hydraulic propertiesof crusted soils from tension disc infiltrometer data. Geoderma 86:61–81.
Šimůnek, J., and M.Th. van Genuchten. 1996. Estimating unsaturated soil
hydraulic properties from tension disc infiltrometer data by numericalinversion. Water Resour. Res. 32:2683–2696.
Šimůnek, J., and M.Th. van Genuchten. 1997. Estimating unsaturated soil
hydraulic properties from multiple tension disc infiltrometer data. Soil
Sci. 162:383–398.
Šimůnek, J., O. Wendroth, and M.Th. van Genuchten. 1999. Estimating
unsaturated soil hydraulic properties from laboratory tension disc
infiltrometer experiments. Water Resour. Res. 35:2965–2979.
Smettem, K.R.J., and B.E. Clothier. 1989. Measuring unsaturated sorptivity
and hydraulic conductivity using multiple disk permeameters. J. Soil Sci.
40:563–568.
Thony, J.L., G. Vachaud, B.E. Clothier, and R. Angulu-Jaramillo. 1991. Field
measurements of the hydraulic properties of soil. Soil Technol. 4:111–123.
van Dam, J.C., J.N.M. Stricker, and P. Droogers. 1992. Inverse method for
determining soil hydraulic functions from one-step outflow experiments.
Soil Sci. Soc. Am. Proc. 56:1042–1050.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
van Genuchten, M.Th., and F.J. Leij. 1992. On estimating the hydraulicproperties of unsaturated soils. p. 1–14. In M.Th. van Genuchten etal. (ed.) Indirect methods for estimating the hydraulic properties of
unsaturated soils. Univ. of California, Riverside.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC code forquantifying the hydraulic functions of unsaturated soils. USEPA Rep.
600/2–91–065 (IAGDW1293–3-934). U.S. Salinity Lab., Riverside, CA.
Warrick, A.W. 1992. Models for disc infiltrometers. Water Resour. Res.
28:1319–1327.
Warrick, A.W. 1993. Inverse estimations of soil hydraulic properties withscaling: One-dimensional infiltration. Soil Sci. Soc. Am. J. 57:631–636.
Watson, K.K. 1966. An instantaneous profile method for determining the hydraulicconductivity of unsaturated porous materials. Water Resour. Res. 2:709–715.
White, I., and M.J. Sully. 1987. Macroscopic and microscopic capillary lengthand time scales from field infiltration. Water Resour. Res. 23:1514–1522.
Wooding, R.A. 1968. Steady infiltration from a shallow circular pond. Water
Resour. Res. 4:1259–1273.
Zhang, R. 1998. Estimating soil hydraulic conductivity and macroscopic capillary
length from the disk infiltrometer. Soil Sci. Soc. Am. J. 62:1513–1521.

865