Advances in Water Resources 30 (2007) 1973–1980
www.elsevier.com/locate/advwatres

Veriﬁcation of numerical solutions of the Richards equation
using a traveling wave solution
Vitaly A. Zlotnik a,*, Tiejun Wang a, John L. Nieber b, Jirka Šimunek cbaDepartment of Geosciences, University of Nebraska-Lincoln, NE 68588, United States
Department of Bioproducts and Biosystems Engineering, University of Minnesota, St. Paul, MN 55108, United Statesc
Department of Environmental Sciences, University of California-Riverside, CA 92507, United States

Received 31 March 2006; received in revised form 20 March 2007; accepted 23 March 2007
Available online 19 April 2007

Abstract
Eﬀorts to ﬁnd solutions that can be used for veriﬁcation of numerical techniques for solving the Richards equation have generated awealth of approximate and exact analytical solutions. Coeﬃcients of this equation involve two highly non-linear functions related to thesoil water potential, the unsaturated hydraulic conductivity, and the soil water content. The known exact solutions for realistic ﬂow
geometries are commonly limited to simpliﬁed descriptions of unsaturated hydraulic properties, while the approximate solutions involve
various simpliﬁcations that require additional veriﬁcation. We present a technique, referred to as the ‘‘launch pad’’ technique, which isbased on the traveling wave solution to generate an exact solution of the boundary value problem for the Richards equation. The technique that is applicable to any descriptor of unsaturated hydraulic properties is illustrated on an application involving the inﬁltration ofwater into soils with properties described by Brooks–Corey and van Genuchten models. Examples of veriﬁcation are presented for
HYDRUS-1D, a popular numerical computer code for solving the Richards equation.
 2007 Elsevier Ltd. All rights reserved.
Keywords: Richards equation; Soil descriptors; Traveling wave; Numerical solution; Code veriﬁcation

1. Introduction
The Richards equation for soil moisture movement
occupies a very important place in modern theoreticaland applied hydrology. Most commonly, tools for solving various hydrologic problems for vadose zone ﬂows
utilize numerical methods based on either the ﬁnite-diﬀerence or ﬁnite-element techniques. Recently, Scanlon et al.
[18] provided an inter-comparison of the most commonlyused codes in the US (HELP, HYDRUS-1D, SoilCover,
SHAW, SWIM, UNSAT-H, and VS2DTI) by simulating
the water balance of surﬁcial sediments in a semiarid
region. The purpose of this comparison was to assess
eﬀects of diﬀerent approaches in treatment of the hydrologic assumptions for ﬂuxes at the upper boundary (inﬁl*Corresponding author. Tel.: +1 402 472 2495; fax: +1 402 472 4917.
E-mail address: vzlotnik@unl.edu (V.A. Zlotnik).

0309-1708/$ - see front matter  2007 Elsevier Ltd. All rights reserved.
doi:10.1016/j.advwatres.2007.03.008

tration and evaporation), the lower boundary (drainage,etc.), and aquifer parameters. Vanderborght et al. [23]reviewed the veriﬁcation aspects of both ﬂow and transport model components and included several additional
codes that are used in the European Union for both academic and applied purposes. A common feature of all
numerical codes, whether those on the list above or others, is the lack of good benchmark tests for veriﬁcationof accuracy of solution methods for transient ﬂow
conditions.
The motivation of our study was to ﬁnd an exact but
simple transient solution for soil moisture dynamics thatcan be used for any descriptor of unsaturated hydraulicproperties and simple but realistic initial and boundary
conditions, including relations described by tabular data.
The requirement that the solution be valid for any descriptor of unsaturated hydraulic properties is especially important because various available codes use a variety of1974

V.A. Zlotnik et al. / Advances in Water Resources 30 (2007) 1973–1980

descriptors. It is desirable that such a solution be derivablewith minimal computational eﬀort.
During the last four decades a large number of analytical approaches have been developed (e.g., Barry et al. [2],Broadbridge and Rogers [3], Fleming et al. [8], Parlange[11,12], Philip [13–15], Pullan [16], Ross and Parlange[17], Smith et al. [21], Srivastava and Yeh [22], Warrick[24], and Witelski [25]). Most of these analytical solutions
have been developed and tested for speciﬁc unsaturated
hydraulic properties. Only Philip [13–15] and Ross andParlange [17] developed solutions for general descriptorsof unsaturated hydraulic properties and one-dimensional
ﬂow in a semi-inﬁnite domain. There are also diﬀerences
between these solutions.
Philip [13] provided a solution for inﬁltration withconsatant soil moisture content at the surface into asemi-inﬁnite soil proﬁle. The solution is given as an expansion series for square-root-time. A disadvantage of thesolution is that series converge only for ﬁnite time range[11,13] and the evaluation of series coeﬃcients is complicated by solving an inﬁnite system of integral diﬀerential
equations. The complexity of numerical evaluation makes
it prohibitive for systematic use, and in fact, characteristics
of only one soil (Yolo light clay) found use in the Philip’s
solution [13–15], with rare exceptions [26]. However, Philip
proved that after large time this solution has a self-similar
front, called it ‘‘solution at inﬁnity’’, and derived its analytical form [15]. The front velocity for this solution is a function of soil descriptors and boundary conditions. The
transition time between small and large time solutions
has not been established [15]. This solution (or so-called
traveling wave [25]) can be considered as a solution of initial boundary value problem with given soil moisture content/head at z = ±1.
Ross and Parlange [17] proposed using a self-similarsolution for a speciﬁed ﬂux at the upper boundary in asemi-inﬁnite domain to test numerical schemes. However,the link with the Philip [15] solution has not been established. In [17], the front velocity was assigned arbitrarily
for a wide range of parameters. It was also shown that surface ponding could occur for some ﬂux values, after which
this solution requires adjustments.
Traveling waves have been well known in ﬂuid mechanics previously [1]. In hydrology, applications of travelingwave solutions are of interest in experimental studies(e.g., [19]) and inﬁltration theory (e.g., [6,7]). However,the practical ramiﬁcations of traveling waves and validation of numerical codes in particular have been underestimated [21].
In this article we combine approaches of Philip [15]and Ross and Parlange [17] for developing a simple andﬂexible method that yields an exact solution that satisﬁes
the initial-boundary value problem for the Richards
equation. This approach has the following traits: (1) it
is valid for any descriptor of the hydraulic properties of
unsaturated soils as in [15,17]; (2) the solution utilizes
velocity c, which is deﬁned [15] by soil properties andboundary and initial conditions, and is valid for arbitrary
time; (3) the solution is presented for either speciﬁed pressure or speciﬁed ﬂux upper boundary condition. We alsopresent technical aspects associated with the evaluation of
integrals involved in the assessment of the traveling wave
proﬁle.
These solutions were applied to soils with the most general unsaturated soil hydraulic properties used in vadoze
zone hydrology such as the Brooks–Corey and van
Genuchten. As an example, we use this method for veriﬁcation of one of the most popular computer codes for simulation of soil moisture ﬂow in vadose zone, HYDRUS-1D[20]. This code has been used in numerous applications;nevertheless traveling wave solutions are shown to be avaluable instrument in clarifying some aspects of the code
accuracy. This means that such solutions can be a goodbenchmark for any test.
2. Traveling wave solution in inﬁnite and semi-inﬁnite spatial
domains
2.1. The problem statement for semi-inﬁnite domain
We consider vertical inﬁltration into an initially moist
bare soil without surface saturation, using the Richards
equation written in water content-based form


oh
o
oh
oKðhÞ
¼
DðhÞ
;
ð1Þ

ot oz
oz
oz
where hydraulic conductivity K(h), diﬀusivity D(h) =K(h)(dh/dh), and retention curve h(h) are non-linear
functions, h is the volumetric water content, h is soil water
potential, t is time, and z is a vertical coordinate oriented in
the direction of gravity. Examples of the most commonlyused soil moisture characteristics are given in Table 1.
The initial soil moisture distribution h0(z) in semi-inﬁnite domain (z > 0) is non-uniform
hðz; 0Þ ¼ h0 ðzÞ;ð2Þz > 0;and the transient upper boundary condition is as follows:
hð0; tÞ ¼ h1 ðtÞ;

ð3Þ

t>0

Our task is to ﬁnd an exact solution that satisﬁes the initialboundary value problem (1)–(3).
2.2. Analytical solution for inﬁnite domain
Philip [14,15] obtained a solution for the Richards equation as a traveling wave that moves from z = 1 toz = +1 with some velocity c. The new variablehðz; tÞ ¼ Hð1Þ; wheref ¼ z  ct

ð4Þ

and the wave velocity c is deﬁned by h = h1 at z = 1 andh = h0 at z = +1 as follows:c¼

Kðh1 Þ  Kðh0 Þ
h 1  h0ð5Þ

V.A. Zlotnik et al. / Advances in Water Resources 30 (2007) 1973–1980

1975

Table 1
Soil moisture characteristics
Characteristics

FormulaParameters

van GenuchtenhðhÞ  hr
1¼ Se ¼;½1þ j ahjn mhs  hrh<0m 2KðhÞ ¼ K s S le ½1  ð1  S 1=me Þ 1
m¼1 ; n>1nBrooks–CoreyhðhÞ  hr¼ Se ¼hs  hrj hhe jn ;1;

hr and hs residual and saturated soil water contents, respectivelyKs
saturated hydraulic conductivity
Seeﬀective water saturationaparameter related to air-entry valuen, mpore size distribution indiceslpore connectivity parameterheinverse of air-entry valuenlh < heh P hepore-size distribution index
pore-connectivity parameterKðhÞ ¼ K s S 2=nþlþ2e

Introduction of variable f into (1) leads to the ordinary differential equation


dH
d
oH
dKðHÞ
¼
DðHÞ
ð6Þ
c

d1
d1
o1
d1
The exact solution in an implicit form for f(H), can be obtained by integrating Eq. (6)1 ¼ F ðHÞ
¼ ðh1  h0 Þ
Z h1
DðuÞ
du;

ðh1  h0 ÞKðuÞ  ðh1  uÞKðh0 Þ  ðu  h0 ÞKðh1 Þhh0 6 h 6 h1 ;ð7Þwhere h1 ¼ h1  D1 and h0 ¼ h0 þ D0 , and D0 and D1 aresmall positive values. These values regularize the integral
F(H), because in general the integrand may have singularities at both integration limits h = h0 and h = h1. Eq. (7)can be integrated in a closed form for simple descriptorsof K(h) and h(h) [17,25], but the commonly used descriptors(Table 1) must be treated numerically.
2.3. Solution evaluation
Theoretically, the wave front length (i.e., the distance
between points with h0 and h1) is inﬁnite, stretching fromz = 1, where h(z, t) = h1 to z = 1 where h(z, t) = h0 at
any given moment of time. This is reﬂected in Eq. (7) as follows from the integrand singularities. With h0 = const andh1 = const > h0, this wetting front moves at velocity c, without changing shape. To make this solution of practical value,
it is necessary to consider only the range of soil moisture contents between h0 and h1 to avoid integrand singularities.

We will use the term ‘‘truncated front’’ for this portionof the solution H(f) or h(z, t) lying in the range ofh0 6 h 6 h1 ; where the most spatial changes in H(f) orh(z, t) occur. The shape of the truncated front is exact
between h0 and h1 and one can select positive D0 ¼
h0  h0 and D1 ¼ h1  h1 as small as desired. Reducing D0
and D1 will increase spatial distance between the endpointsof the solution and expand the extent of a portion of theexact solution, but not change the front shape otherwise.
The function h(z, t) in domains [h0, h0 Þ and (h1 , h1] can be
replaced by h0 and h1 , respectively with any desired accuracy by reducing D0 and D1.
We will illustrate this eﬀect of varying parameters D0
and D1 in the integration of Eq. (7) using its implementation in the Maple V language [5]. Parameters for the
sandy soil moisture characteristics (Table 2) were selectedfrom the HYDRUS-1D database that provides averagesoil hydraulic parameters for diﬀerent USDA texturalclasses [4].The front portions of the wave proﬁle for boththe Brooks–Corey and van Genuchten models are shownin Fig. 1 for the parameters given in Table 3 (h0, h1, andtwo diﬀerent values of D0 = D1). A decrease in D0 = D1
leads to a change in corresponding f values and separation of the two proﬁles. Note that the dramatic diﬀerencein slope and curvature behaviors of proﬁles near h ¼ h0and h ¼ h1 is attributed to non-linearity of the soil
descriptors.
When applying this regularization procedure, theremaining portions of the traveling wave in domains[h0, h0 Þ and (h1 , h1] are approximated by horizontal lines
starting at hðz; tÞ ¼ h1 and hðz; tÞ ¼ h0 . The actual slopesof the wetting front in these sections can be determinedfrom (dh/df) = (df/dh)1, which approaches zero when soil

Table 2
Parameters for soil moisture characteristics
Characteristics

hs

hr

a

he

n

l

Ks (cm/day)van Genuchten
Brooks–Corey

0.43
0.417

0.045
0.02

0.145
–

–
0.138

2.68
0.592

0.5
1

712.8
504

1976

V.A. Zlotnik et al. / Advances in Water Resources 30 (2007) 1973–1980

Fig. 1. Eﬀect of regularization on computation of truncated proﬁles of soil moisture contents for various soil descriptors: (a) Brooks–Corey model,(b) area of detail for Brooks–Corey model, (c) van Genuchten model, (d) area of detail for van Genuchten model.

Table 3
Parameters for calculating soil moisture proﬁles
h0h1

D0 = D1

0.05

0.25

105

108

moisture content approaches h0 or h1, i.e., D0 or D1
approaching zero. For example, slopes of the wetting front
using the Brooks–Corey model are 6.845 · 102 m1 and6.853 · 105 m1 at point h ¼ h0 for D0 = 105 and 108,respectively, and 1.699 · 104 m1 and 1.699 ·107 m1 at points h ¼ h1 for D1 = 105 and 108, respectively (see Fig. 1c). This shows that the remaining part ofthe solution can be replaced by h  h0 for h < h0 andh  h1 for h > h1 with any desirable accuracy by controlling D0 and D1, because the actual soil moisture contents
h in these parts of the proﬁle is almost uniform.
The diﬃculty of using the traveling wave solution in
the original form for practical purposes is that this proﬁleis not associated with any speciﬁc moment of time [21].
However, a remarkable property of the truncated front
is that it oﬀers an opportunity to construct analyticalsolutions in the semi-inﬁnite domain for speciﬁc momentsof time t, if the initial front (e.g., t = 0) position is
known. Further evolution of the proﬁle is achieved bymere translation of this portion of the front with known
velocity c [7]. We will call this initial position of a truncated front a ‘‘launch pad’’. Naturally, the boundary
condition at z = 0 must be consistent with the front
translation. Any initial position of proﬁle, or a launchpad, at t = 0, provides a new solution to initial-boundary
value problems with any desired accuracy as will be shown
below.
2.4. Selection of parameters in the analytical solution
The desired solution accuracy determines selection of D0
and D1. Typically, this accuracy is estimated as a fraction of
the h range. For practical purposes, D0  D1 = ejh0 6 hj,
where e = 103–104 far exceeds the accuracy of ﬁeld andlaboratory measurements of h. Decreasing D0 and D1 shifts
and extends the length of the ‘‘proﬁle fraction’’ as shown in
Fig. 1d. Asymptotic expansion of the proﬁle f = F(H) nearh0 and h1 indicates that further increase in accuracy (reduction in e, D0 and D1) only weakly extends the ‘‘proﬁle frac-

V.A. Zlotnik et al. / Advances in Water Resources 30 (2007) 1973–1980

tion’’ in the parts of the proﬁle near h0 and h1 (as log e).
These parts are trivial and not important for analyses orapplications (see Fig. 1d).
3. Traveling wave solution in semi-inﬁnite-domain
3.1. Solution for uniform initial soil moisture content(transient conditions at the boundary)z > 0;

t¼0

Note that h0 and h1 can approximate h0 and h1, respectively
with any degree of accuracy.
The upper boundary condition can be formulated also
in terms of ﬂux. The ﬂux at the soil surface can be derived
directly from the calculated surface soil moisture content
by applying the deﬁnition of ﬂux on the soil surface andthe traveling wave deﬁnition from Eq. (4)qð0; tÞ ¼ ½ðH  h1 ÞKðh0 Þ  ðH  h0 ÞKðh1 Þ þ 2ðh1  h0 Þ

We consider the launch pad outside the ﬂow domain(z < 0 at t = 0). In this case, the initial soil moisture proﬁleis uniform
hðz; 0Þ ¼ h0 ;

1977

ð8Þ

The truncated proﬁle is set ‘‘outside’’ the ﬂow domain
initially, namely in the domain 1 < z 6 0 similarly asin Ross and Parlange [17]. As time progresses (t0 <t1 < t2. . .), soil moisture content and ﬂux at the upperboundary vary with time, and the proﬁle ‘‘slides’’ into thedomain (Fig. 2a). In our case, the constant translation
velocity c is determined by soil descriptors and initial andboundary conditions in contrast with [17]. The upper transient boundary condition for soil moisture content, h(0, t),must comply with the passage of the front through theground surface z = 0. Therefore,8 >< h0 ; t ¼ 0
hð0; tÞ ¼ F 1 ðctÞ; 0 < t < tS
ð9Þ
>
: 
h1 ; t > t S
After time tS, the soil moisture at the boundary becomessteady at hð0; tÞ ¼ h1  h1 , t > tS (Fig. 2b).
Estimation of time tS can be performed using function
F(h) from Eq. (7) and its inverse function F1(f). Indeed,time tS is found from relationships
fjz¼0 ¼ ct ¼ F ðhð0; tÞÞ; ctS ¼ F ðh1 Þ; and tS ¼ F ðh1 Þ=cð10Þ KðHÞjz¼0 =ðh1  h0 Þ;

Hjz¼0 ¼ hð0; tÞð11ÞNote, that q(0,t) ! K(h1) for large times because H ! h1.
3.2. Non-uniform initial soil moisture proﬁle (constantmoisture content at the boundary)Selection of the launch pad inside the ﬂow domain generates a solution for a speciﬁc non-uniform initial head/soil
moisture distribution. In Fig. 2a, the truncated proﬁle att > t3  tS is practically located within the ﬂow domainz > 0. This proﬁle can be taken as another initial condition
for z > 0 and t = 0, as done recently by Vanderborght et al.
[23]). As time progresses, the proﬁle translates downward
with a constant velocity c.
The initial soil moisture distribution is constructedaccording to the traveling wave solution (Fig. 1) as
described in Section 2.2, that is,hðz; 0Þ ¼ HðfÞ;z > 0;

t¼0

ð12Þ

The upper boundary condition is then one of constant speciﬁc value of head/soil moisture contents and constant inﬁltration ﬂux, and is given by
hð0; tÞ ¼ h1 ;t > 0 or q ¼ Kðh1 Þ;

t>0

ð13Þ

Note that h1 may approximate the use of h1 with any
desired accuracy because D1 is negligibly small. With decrease in D0 and D1, our constructed solution becomesexact.

Fig. 2. Exact solution with uniform initial proﬁle h(z, 0) = h0 and the passage of the front through z = 0 s: (a) Truncated proﬁle is set fully outside of theﬂow domain, Brooks–Corey model, z < 0; (b) The upper transient boundary condition h(0, t), and tS  0.29 day.

1978

V.A. Zlotnik et al. / Advances in Water Resources 30 (2007) 1973–1980

4. Traveling wave application for veriﬁcation of a numericalcode

Table 4
Parameters for calculating wetting fronts in initially dry soils by travelingwave solution (TWS) and HYDRUS-1D

Comparative studies [18,23] have not distinguished adistinctly superior code among the selected group of themost commonly used codes. Each computer code listed
includes a complex interaction of numerical method, griddesign, and input data preparation including computations
or interpolation of soil descriptors for numerical simulation (e.g., [9,10]).
For example, one of the most well documented andwidely used codes, HYDRUS-1D, either evaluates soil
hydraulic properties directly from analytical models or by
the linear interpolation from internal tables that are generated at the beginning of a numerical simulation [20]. Whenthe internal tables are used, then HYDRUS-1D generates
within a speciﬁed interval of pressure heads [ha,hb] (orlower and upper limits of tension, LLT and ULT, respectively) a table of water content, hydraulic conductivity,and speciﬁc water capacity values from the speciﬁed setof soil parameters. One hundred logarithmically-varying
intervals are used for interpolation. Values of these characteristics are then computed during the iterative solution
process using linear interpolation between entries in thetable. If an argument h falls outside the prescribed interval(ha, hb) or ha = hb, the hydraulic characteristics areevaluated directly from the hydraulic functions. The above
interpolation technique was found to be much fastercomputationally than direct evaluation of the hydraulicfunctions over the entire range of pressure heads, while preserving the precision of calculated results in a majority of
applications. Because of the highly nonlinear relationships
between h, h, and K, the degree of inﬂuence of LLT andULT on simulation of the wetting front generally is notapparent.
Overall code optimization cannot be achieved by mere
increase in accuracy of each code component (numericalmethod, discretization, and interpolation). One advantage

Characteristics

h0 in TWS

h0 in HYDRUS-1D

h1

D0 = D1

van Genuchten
Brooks–Corey

0.045
0.02

0.045001
0.02001

0.25
0.25

105
105

of using analytical methods is that they provide an opportunity for sensitivity analysis and quantify inﬂuence of all
steps of the numerical solution on simulation results inentirety. In this aspect, ﬂexibility of the traveling wave
solution is very attractive as a benchmark test.
The proposed traveling wave solution for uniform initialsoil moisture proﬁles was compared with simulations performed using HYDRUS-1D. Usually, simulation ofadvancement of a steep wetting front into initially dry proﬁle is a stringent test of a numerical solution of the Richards equation. Both Brooks–Corey and van Genuchten
models were tested with hydraulic parameters given in
Table 2 and parameters for calculating traveling wave inTable 4. Note that the initial soil moisture content was
set to be slightly higher than the residual soil moisture contents for convergence. Simulations for 2-day inﬁltration
events are presented in Fig. 3. For van Genuchten model,
h0 (corresponding to h0) and h1 (corresponding to h1) are1.46 · 104 cm and 8.47 cm, respectively, and for
Brooks–Corey model, h0 and h1 are 4.25 · 108cm and18.22 cm, respectively.
The ﬁrst simulations used the default values of LLT =106 cm and ULT = 104 cm in HYDRUS-1D. The results(Fig. 3) show that when comparing positions of wetting
fronts at h  h0, the simulated wetting front fromHYDRUS-1D are about 6.3 and 20.1 cm/day faster thanvelocities of the traveling wave for the Brooks–Corey andvan Genuchten models, respectively (Note that the calculated discrepancy of the position of wetting fronts is also
dependent on spatial grid discretization in HYDRUS-

Fig. 3. Comparisons of numerical and analytical solutions for the uniform initial proﬁle problems: (a) Brooks–Corey model and (b) van Genuchten
model.

V.A. Zlotnik et al. / Advances in Water Resources 30 (2007) 1973–1980

1D). When soil hydraulic properties are evaluated directlyfrom analytical models (by specifying ULT larger than theboundary value jh1j, corresponding to h1 (e.g., ULT =1.0 cm), the discrepancy of wetting fronts is signiﬁcantly
reduced to 2.3 and 1.62 cm/day for Brooks–Corey andvan Genuchten models, respectively. The results conﬁrm
the notion that evaluating the soil hydraulic properties byinternal interpolation instead of direct use of the analytical
functions requires a tradeoﬀ between computational speed
and accuracy.
5. Conclusions
We have revisited the exact traveling wave solutions for
the Richards equation developed by Philip [14,15] andRoss and Parlange [17]. The self-similarity of travelingwave solutions allows to generate diﬀerent exact solutions
for inﬁltration problems for arbitrary ﬁnite times using
the ‘‘launch pad’’ concept. After computing a major portion of the exact solution called a ‘‘truncated proﬁle’’, thelaunch pad for this truncated proﬁle can be set initially inany part of the ﬂow domain. If the launch pad for the truncated proﬁle is initially located fully within the ﬂow domain(z > 0), the solution for inﬁltration with constant moisture
content and ﬂux at the boundary is generated. If the launch
pad for the truncated proﬁle is located outside of thedomain (z < 0), this approach generates a solution for inﬁltration into the initially uniform soil moisture proﬁle with acertain transient boundary condition, when the surface soil
moisture content is consistent with the traveling wave
proﬁle. We presented the procedure for evaluating thetruncated proﬁle and analytical expression of this timedependent boundary condition. This approach is valid
for arbitrary soil characteristics and for speciﬁc soil-moisture content-based or ﬂux-based boundary conditions at
the surface.
The practical use of the launch pad technique in semiinﬁnite vertical soil proﬁle was illustrated by assessmentof the eﬀects of computational parameters of a commonlyused code, HYDRUS-1D for Brooks–Corey and van
Genuchten soil descriptors. Overall veriﬁcation of numerical codes with interplay of numerical method, discretization, and interpolation can be facilitated by using ﬂexible
traveling wave solutions as benchmark tests. Such an
approach can be utilized for numerical codes with any
realizable set of unsaturated hydraulic property
descriptors.
Acknowledgements
We thank G.I. Barenblatt (UC Berkeley) and T. Witelski (Duke University) for fruitful discussion. This study
was partially supported by NSF’s Biocomplexity in theEnvironment Program BE03-22067 (V. Zlotnik); and theTerrestrial Sciences Program of the Army Research Oﬃce
(J. Šimunek).

1979

References[1] Barenblatt GI. Scaling, self-similarity, and intermediate asymptotics. Cambridge: Cambridge University Press; 1996.
[2] Barry DA, Parlange JY, Sander GC, Sivaplan M. A class ofexact solutions for Richards equation. J Hydrol 1993;142:
29–46.
[3] Broadbridge P, Rogers C. Exact solutions for vertical drainage andredistribution in soils. J Eng Math 1990;24:25–43.
[4] Carsel RF, Parrish RS. Developing joint probability distributions ofsoil water retention characteristics. Water Resour Res
1988;24:755–69.
[5] Char BW, Geddes KO, Gonnet GH, Leong BL, Monagon MB,
Watt SM. Maple V language reference manual. Springer-Verlag;
1991.
[6] DiCarlo DA. Modeling observed saturation overshoot with continuum additions to standard unsaturated theory. Adv Water Res
2005;28:1021–7.
[7] Egorov AG, Dautov RZ, Nieber JL, Sheshukov AY. Stability
analysis of gravity-driven inﬁltrating ﬂow. Water Resour Res
2003;39:1266. doi:10.1029/2002WR001886.
[8] Fleming JF, Parlange JY, Hogarth WL. Scaling of ﬂux and watercontent relations: comparison of optimal and exact results. Soil Sci
1984;137:464–8.
[9] Haverkamp R, Vauclin M, Touma J, Wierenga PJ, Vachaud G. A
comparison of numerical simulation models for one-dimensional
inﬁltration. Soil Sci Soc Am J 1977;41:285–93.
[10] Milly PCD. Advances in modeling of water in the unsaturated zone.
Transp Porous Media 1988;3:491–514.
[11] Parlange JY. Theory of water-movement is soils: 2. One dimensional
inﬁltration. Soil Sci 1971;111:170–4.
[12] Parlange JY. Theory of water-movement is soils: 8. One dimensional
inﬁltration with constant ﬂux at the surface. Soil Sci 1972;114:1–
4.
[13] Philip JR. The theory of inﬁltration: 1. The inﬁltration equation andits solution. Soil Sci 1957;83:345–57.
[14] Philip JR. The theory of inﬁltration: 2. The proﬁle at inﬁnity. Soil Sci
1957;83:435–48.
[15] Philip JR. Theory of inﬁltration. Adv Hydrosci 1967;5:215–305.
[16] Pullan AJ. The quasilinear approximation for unsaturated porous
media ﬂow. Water Resour Res 1990;26:1219–34.
[17] Ross PJ, Parlange JY. Comparing exact and numerical solutions ofRichards equation for one-dimensional inﬁltration and drainage. Soil
Sci 1994;157:341–4.
[18] Scanlon BR, Christman M, Reedy RC, Porro I, Simunek J,
Flerchinger GN. Intercode comparisons for simulating water balanceof surﬁcial sediments in semiarid regions. Water Resour Res
2002;38:1323. doi:10.1029/2001WR001233.
[19] Selker J, Parlange JY, Steenhuis T. Fingered ﬂow in two dimensions,2, predicting ﬁnger moisture proﬁle. Water Resour Res
1992;28(9):2523–8.
[20] Simunek J, Sejna M, van Genuchten MT. The HYDRUS-1D
software package for simulating the one-dimensional movement ofwater, heat, and multiple solutes in variably-saturated media, Version
3.0, HYDRUS Software Series 1, Department of Environmental
Sciences, University of California-Riverside, Riverside, CA, 270,
2005.
[21] Smith RE, Smettem KRJ, Broadbridge P, Woolhiser DA. Inﬁltration
theory for hydrologic applications. Washington, DC: American
Geophysical Union; 2002.
[22] Srivastava R, Yeh TCJ. Analytical solutions for one-dimensional,transient inﬁltration toward the water table in homogeneous andlayered soils. Water Resour Res 1991;27:753–62.
[23] Vanderborght J, Kasteel R, Herbst M, Javaux M, Thiery D,
Vanclooster M, et al. A set of analytical benchmarks to test
numerical models of ﬂow and transport in soils. Vadose Zone J
2005;4:206–21.

1980

V.A. Zlotnik et al. / Advances in Water Resources 30 (2007) 1973–1980[24] Warrick AW. Soil water dynamics. Oxford University Press;
2003.
[25] Witelski TP. Motion of wetting fronts moving into partially pre-wet
soil. Adv Water Res 2005;28:1133–41.

[26] Zlotnik VA, Murashko MG. Comparison of two methods forcalculation of one-dimensional vertical inﬁltration. In: Soil Physics
Applications for Soil Reclamation, Trudy Agrophysicheskogo Instituta, vol. 31, USSR: Leningrad; 1973. p. 234–7 (in Russian).