WATER RESOURCES RESEARCH, VOL. 33, NO. 9, PAGES 2049 –2063, SEPTEMBER 1997

New piecewise-continuous hydraulic functions for modelingpreferential flow in an intermittent-flood-irrigated field
B. P. Mohanty,1 R. S. Bowman,2 J. M. H. Hendrickx,2 and M. T. van Genuchten1
Abstract. Modeling water flow in macroporous field soils near saturation has been amajor challenge in vadose zone hydrology. Using in situ and laboratory measurements, wedeveloped new piecewise-continuous soil water retention and hydraulic conductivity
functions to describe preferential flow in tile drains under a flood-irrigated agricultural
field in Las Nutrias, New Mexico. After incorporation into a two-dimensional numericalflow code, CHAIN_2D, the performance of the new piecewise-continuous hydraulicfunctions was compared with that of the unimodal van Genuchten-Mualem model andwith measured tile-flow data at the field site during a number of irrigation events. Modelparameters were collected/estimated by site characterization (e.g., soil texture, surface/subsurface saturated/unsaturated soil hydraulic property measurements), as well as by local
and regional-scale hydrologic monitoring (including the use of groundwater monitoring
wells, piezometers, and different surface-irrigation and subsurface-drainage measurement
systems). Comparison of numerical simulation results with the observed tile flow indicatedthat the new piecewise-continuous hydraulic functions generally predicted preferential flowin the tile drain reasonably well following all irrigation events at the field site. Also, thenew bimodal soil water retention and hydraulic conductivity functions performed better
than the unimodal van Genuchten–Mualem functions in terms of describing the observed
flow regime at the field site.

1.

Introduction

Preferential flow under saturated [Logsdon, 1995; Andreiniand Steenhuis, 1990; Everts and Kanwar, 1990; Jardine et al.,
1989; White, 1985; Beven and Germann, 1982] as well as unsaturated [Booltink, 1993; Jardine et al., 1990; Wilson et al., 1990;
Phillips et al., 1989; Bowman and Rice, 1986; White, 1985; Beven
and Germann, 1982; De Vries and Chow, 1978] conditions
through macropores, cracks, and other nonmatrix domainsthat coexist with the soil matrix domain has been found to be
more a rule rather than an exception in structured field soils.
Multidomain models are being developed [Zurmuhl andDurner, 1996; Gwo et al., 1996; Hutson and Wagenet, 1995;
Gerke and van Genuchten, 1993a, b] to account for this nonideal flow behavior and related solute transport processes invariably saturated soil. Multidomain flow/transport models
typically need soil hydraulic conductivity (K i (h i ), i 5 1,
2, z z z , n) and soil water retention ( u i (h i ), i 5 1, 2, z z z , n)relationships as a function of the pressure head, h i , for each of
the n domains, including terms accounting for the interaction
or exchange of water or solutes between the different domains.
Although multidomain approaches are often based on aconceptual analysis of the underlying flow and transport processes [Jardine et al., 1990], most available techniques for measuring soil hydraulic properties (conductivity and retention)
can neither distinguish between the different flow domains andtheir relative contribution to flow [Luxmoore et al., 1990], nor
be used to determine the between-domain exchange terms
1

U.S. Salinity Laboratory, USDA-ARS, Riverside, California.
New Mexico Institute of Mining and Technology, Socorro.

2

Copyright 1997 by the American Geophysical Union.
Paper number 97WR01701.
0043-1397/97/97WR-01701$09.00

during variably saturated flow, except perhaps when using
some type of inverse procedure built upon many simplifying
assumptions [Durner, 1994]. For purposes of multidomain
flow-transport modeling, once the hydraulic functions (K-hand u -h) are known for one domain (the bulk soil), differentu -h and K-h relationships may be derived for different domains by simply multiplying or dividing the one-domain functions with certain factors. These factors could be selectedbased upon the statistics of the fractional cross-sectional areasof the domains being considered. For example, Hutson andWagenet [1995] created domain-specific hydraulic functions fortheir five-domain flow and transport model by applying lognormally distributed scaling factors (having a mean of 1.0 anda standard deviation of 0.25) to the bulk soil K(h) and u (h)curves, thus creating an ensemble of five different hydraulicdomains of equal area. Their approach in turn resulted inmean K(h) and u (h) curves for the five domains similar to
those for the single-domain soil. Although this approach ofdefining characteristic functions for different hydraulic domains may be mathematically appropriate, the task of testingthese functions is formidable in view of current experimental
limitations.
Some researchers [e.g., Anderson and Hopmans, 1994; Singhet al., 1991a, b; Warner et al., 1989; Edwards et al., 1988; Bouma,
1981; Bouma and Dekker, 1978] attempted to quantify macroporosity, including pore dimensions and pore abundance indifferent soil horizons, by direct measurement using computer
tomography, AUTOCAD, image analysis, resin impregnation,or some other method. These studies were aimed at modelingmacropore flow on the basis of the number and geometry ofpores, as was suggested by Edwards et al. [1979]. On the other
hand, Ela et al. [1992] suggested that knowledge of the number
and size of visible macropores by itself may be insufficient tomodel water infiltration through macropores. A more practical

2049

2050

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

Figure 1. A close view of the pore system in the bulk soil
under disc infiltrometer during in situ hydraulic conductivity
measurement. K(h) near saturation is the composite value forthe bulk soil including micropores, mesopores, macropores,
and cracks. Thickness of irregular lines indicates pore size.

indirect approach was proposed by Watson and Luxmoore[1986] and later advocated also by Wilson et al. [1992] andMohanty et al. [1996], among others. They used the differenceof the hydraulic conductivity measured at 0-cm and 3-cm soil
water tensions using disc infiltrometers to define the macropore conductivity and using Poiseuille’s equation to estimate
the effective macroporosity of certain (.1-mm) diameter
pores as recommended by Luxmoore [1981]. Although thismethod in essence is a lumped approach that cannot yield theinteraction/exchange terms in multidomain models, themethod does permit better measurement of the hydraulicproperties of soils near saturation.
The accuracy of soil water content measurements is generally limited by the error of the measurement method. Unfortunately, it is experimentally very difficult to measure the retention characteristics of a soil accurately and reproducibly atvery low suctions. It is further questionable whether retention
curves can properly represent the pore-size distribution of aporous medium at high saturation where the nonwetting phase(air) is not continuous [Corey, 1992]. This means that hydraulicconductivities predicted from retention data are inherently
inaccurate when approaching saturation, even if one assumes
that the predictive conductivity model itself is error-free andapplies to all pore sizes. As a consequence, hydraulic conductivity estimation from soil water retention functions that represent a certain pore space near saturation (e.g., h less than 1cm) are inherently unreliable, regardless of whether unimodal
or multimodal functions are used. This is generally also thecase when the parameter n in van Genuchten’s retention equation (m 5 1 2 1/n) is close to unity [Durner, 1994].
Disk infiltrometry methods involving ponded and tension
infiltrometers [e.g., Ankeny et al., 1988; Prieksat et al., 1992] are
now increasingly used for measuring hydraulic conductivity insitu at low soil water tensions (i.e., close to saturation). One
advantage of these methods is that negative soil water pressures at the soil-infiltrometer interface can be maintained very

close to zero and be decreased in small increments to get
precise K-h functions near saturation [Jarvis and Messing,
1995; Messing and Jarvis, 1993; Clothier and Smettem, 1990].
Using this macroscopic approach, one important assumption isthat the soil water pressure head across the disc surface willremain uniform and thus in equilibrium with the soil medium
for all “conducting” pore sizes and flow domains, including thelarger macropores. Hence the soil water pressure head maintained in the infiltrometer represents the weighted average forall “conducting” flow domains that coexist below the disc (Figure 1). Since this method assumes a uniform head (i.e., water
maintained at a certain positive or negative soil water pressure) or uniform flux (i.e., water application rate) across thesoil surface, flooding/ponded irrigation (as in our field site)
could provide an attractive condition to test the performanceof alternative new hydraulic conductivity functions at and nearsaturation. Note, however, that no bypass flow, as that which
may occur under nonequilibrium conditions through macropores, is involved in this approach (i.e., all small pores at any
given depth have to be filled before the water is transported
only or mainly through the large pores).
The goal of this study hence is to develop physically realistic,
yet practical (based on actual measurements) soil water retention and hydraulic conductivity functions for describing subsurface water flow and related preferential solute transport
processes in a macroporous agricultural field under irrigation.
Main focus is on the development and testing of new piecewise-continuous hydraulic functions based on in situ and laboratory measurements of the soil water retention and hydraulicconductivity using small increments in h. The new functions
will be used to predict observed preferential water flow into tile
drains following several irrigation events at a field site in Las
Nutrias, New Mexico. Furthermore, performance of the new
bimodal hydraulic functions will be compared with the traditional unimodal van Genuchten-Mualem model. A follow-up
paper will address solute transport during preferential flow at
the same field site.

2.

Field Site and Hydrology Monitoring

Our experimental field site is located on a 24-ha commercial
agricultural farm, situated on an alluvial flood plain of the Rio
Grande, near Las Nutrias, New Mexico. The basin lies within
the Rio Grande rift, a series of north-south trending half grabens oriented parallel to the course of the Rio Grande. The
climate is arid to semiarid and has a wide range of temperatures and rainfall. The parent material for the entire valley isalluvium. Plate 1 shows the major soil classes across the 6-ha
center bench of the field, including our study section. In general, the field site consisted of silty clay loam sediments withmoderate to poor drainage properties, underlain by fine sands,
and with no impeding strata to a depth of about 7 m [Natural
Resource Conservation Service (NRCS), 1992]. Surface horizons, however, contained visible root channels, worm holes,
and cracks, thereby composing a complex network of preferential flow paths. The field was relatively flat and divided into
three benches, namely, east (240 3 270 m), center (240 3
270 m; Figure 2), and west (420 3 270 m), separated by two
berms. To sustain agricultural productivity, the field site was
equipped with surface irrigation and subsurface (tile) drainage
systems to reduce salinity in the crop-root zone. The subsurface tile drainage system consisted of four lateral drain lines
installed approximately at a depth of 1.2 m. The tile drains

Plate 1. Soil series map of the center bench at the Las Nutrias field site, in New Mexico. Study section contains four major soil
series that in turn have different soil textural classes.

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS
2051

2052

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

Figure 2. Three-dimensional view of the field site including important hydrology monitoring facilities andthe locations of surface and subsurface site characterization studies.

connected to a main drain line leading to an off-site openchannel drain on the west. Irrigation return flows draining
from the field, as well as periodic build up of shallow groundwater flowing in from off-site areas, were collected by thissystem. A southern section of the center bench drainage system was isolated by installing two manholes along a single
drain line at the east (up) and the west (down) ends of thebench. The manholes allowed regular sampling of water quantity and quality entering and leaving that (study) section of thedrain line. The manholes were equipped with Signett flow
measurement systems including data loggers for monitoringflow rates at 5-min intervals. We used this isolated section of
the center bench for collecting soil hydraulic parameters andfor monitoring flow in the vadose zone. Horizontal and verticalpressure gradients at different (2-, 5-, and 7-m) depths in theunderlying aquifer were determined regularly using piezometers and observation wells along the periphery of the field site.
For this purpose, eight nested piezometers at 2-, 5-, and 7-m
depths and 33 monitoring wells at 2-m depth (Figure 2) were
installed around the field. Three rain gauges were located on
the center bench to monitor the average precipitation input ofindividual rainfall events. A concrete-lined canal on the northside of the field supplied irrigation water to all three benches
through a number of control gates using flood-irrigation practice. Irrigation water moved from the north end of the field tothe south end in about 4 –7 hours and remained ponded in thefield for several hours with simultaneous infiltration. After therequired amount of water was delivered to the field, irrigation
supply was cut off, thus allowing continued variably saturated
flow and redistribution within the soil profile. A 25-cm diameter circular weir of the type developed by Samani [1993] wasused in the irrigation canal to measure water input to thecenter bench during each irrigation event. The farmer scheduled the water input events to different benches depending on
the consumptive needs of different crops at different growthstages. Table 1 shows the irrigation dates and amounts including the bench sequence during 1994 and 1995. In general, allthree benches received water on the same day, although irrigation in three benches was spread occasionally over 2-day
periods. Table 2 shows the crop rotation for 1994/1995.

3.

Site Characterization

Surface infiltration is by far one of the most important variables to determine the behavior of any field-scale flow andtransport process. Disc infiltrometers including ponded- andTable 1. Irrigation Application Dates and Volumes
Applied to the Center Bench

Date

Total
Elapsed
Time,
Hours

Total
Volume,L

Total
Depth,mm

Bench
Sequence

June 8, 1994
June 27, 1994
Aug. 8, 1994
Sept. 28, 1994
April 10, 1995
April 26, 1995
May 15, 1995
June 5, 1995
June 19, 1995
July 10, 1995
July 28, 1995
Aug. 8, 1995
Sept. 7, 1995
Sept. 26, 1995
Oct. 23, 1995

4.5
7.0
6.0
7.5
6.0
6.5
6.0
6.0
7.5
7.0
7.4
4.8
4.2
4.5
7.2

4.14E106
4.56E106
5.10E106
5.88E106
6.47E106
5.37E106
6.77E106
6.17E106
6.18E106
6.84E106
6.65E106
4.98E106
4.08E106
5.07E106
7.78E106

61.22
67.43
75.41
86.93
95.76
79.48
100.0991.23
91.46
101.26
98.41
73.62
60.37
75.03
115.03

E-C-W
E-C-W
E-C-W
W-C-E
E-W-C
E-W-C
E-W-C
E-C-W
E-C-W
E-C-W
W-C-E
E-W-C
W-C-E
W-C-E
W-C-E

Bench sequence indicates the order of water application during aparticular irrigation event. E, east bench; C, center bench; W, west
bench.

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

tension-infiltrometers were used to measure the K-h relationship in the near-saturation pressure head range (i.e., low soil
water tension range). To incorporate field-scale variability inhydraulic properties resulting from different soil types andpore-size distribution across the field, in situ hydraulic conductivity measurements at 104 selected sites along four transects(as shown in Plate 1) were made. Since the flow regime in thevadose zone and in the shallow groundwater aquifer at thecenter of the study section was two-dimensional with flow lines
converging to the tile line, we implemented a more rigorous
parameter estimation effort across the north-south transect
orthogonal to the tile line, which was oriented in an east-west
direction. On the basis of the tile-drainage design, the particular tile supposedly served a 75-m-wide and 240-m-long section of the center bench. To better define the subsurface flow
field of the section, we selected 61 sites along the north-southtransect located symmetrically across the tile line. Sequential
grid spacings of 2 3 0.5 m, 1 3 1 m, 25 3 1.5 m, 1 3 12 m, and1 3 15 m were used on both the north and south sides of thetile line. Each of the remaining three transects contained 15
measurement locations. All hydraulic property measurements
along the east-west transect were taken directly above thedrain line with a regular spacing of 15 m. The measurements
along the diagonal transects were made at 18-m intervals. All
measurements were later used for estimating field-average(^ &) hydraulic properties using simple (geometric or arithmetic) averaging schemes. Surface infiltration measurements
at 0-, 30-, 60-, and 150-mm soil water tensions were made atthese sites in the summer of 1994, using experimental procedures as described in a similar previous study by Mohanty et al.
[1994, 1996]. To minimize any site disturbance, the watersupply tower of the tension infiltrometer was refilled during theexperiment using the procedure of Ankeny [1992]. In addition
to measurements at the aforementioned soil water tensions, 10
sites across the center bench representing different soil types
were selected for conducting infiltration measurements between 0- to 30-mm soil water tensions at 2- to ;5-mm tension
increments. Ankeny [1992] cautioned that infiltration for tensions below 10 mm should be interpreted carefully because of
the possible inaccuracies in the applied tension owing to bubbling. As such, we interpreted the infiltration rate in this narrow range of tensions only in a qualitative manner so as toreveal relative trends near saturation. The data were used todevelop new piecewise-continuous hydraulic conductivity functions whose fitting parameters could be further adjusted duringmodel calibration. For the above intensive infiltrometer measurements, we used thin layers (1 ; 2 mm) of contact sand tominimize any interfacial effect of the sand on the infiltration
rate close to saturation as reported by Everts and Kanwar[1993]. The actual tensions at the soil surface as reported here
were estimated by assuming unit-gradient water flow throughthe sand layer and by reducing the tension at the base of theinfiltrometer by the average sand depth.
Following the surface infiltration measurements, 5-cmdiameter and 5-cm-long undisturbed soil cores were collected
directly below the disc from the same 10 sites representing
different soil series and textures across the field. The soil cores
were subsequently used for measuring the water retention andhydraulic conductivity for a wider range (0 to 1700 mm) of soil
water tensions. The laboratory method of Richards [1965] andthe multistep outflow approach of Gardner [1956] were used
for measuring soil water retention and hydraulic conductivity
functions, respectively. Similar to our approach in the field, we

2053

Table 2. Crop Rotation Schedule for the Center Bench
Crop

Planting Date

Winter wheat
Sorghum-Sudan crossAlfalfa

Oct. 10, 1993
May 13, 1994
Sept. 25, 1994

HarvestDate(s)not harvested
July 19, 1994
May 29, 1995
June 2, 1995
July 5, 1995
Aug. 6, 1995
Sept. 6, 1995
Oct. 15, 1995

were mostly interested in relative trends rather than absolutevalues of the retention data near saturation. While the actual
numbers may be imprecise, they could be adjusted by subsequent calibration of selected parameters in the new piecewisecontinuous retention function so as to better represent the flow
scenario at the field site. For this purpose the in situ andlaboratory experiments were carried out as carefully as possible using very small increments in the tension near saturation.
For subsurface horizons, water retention and hydraulic conductivity properties at four representative sites in the field(Figure 2) were measured using the instantaneous profile (IP)method [Watson, 1966]. Particle size estimates to a depth of
1.5 m at 200 sites on a regular square grid across the center
bench were used to locate the representative profiles on thebasis of the occurrence of similar textural layering in different
horizons. Since the groundwater table at the field site wasapproximately 40 ; 50 cm lower during spring than in thesummer, we conducted a 3-week long instantaneous profile
experiment during the spring of 1995. Each IP experimental
plot was equipped with a neutron probe (at the center) and 18
tensiometers installed at different depths up to a depth of
1.2 m. Soil water contents and pressure heads at different
depths were measured at 30-min intervals at the beginning of
the experiment; frequency was decreased gradually to several
hours/days on the basis of the rates of change in the observed
data. The data were analyzed to define the soil water retention
and hydraulic conductivity functions at different depths (horizons) across the soil profile for subsequent use as input parameters for the numerical flow simulations.

4. Initial and Time-Dependent Boundary
Conditions at the Field Site
In late February 1994 soil water content measurements weremade on a parallel transect, 1 m off the original north-southtransect used for the surface infiltration measurements (Figure2). The same transect was followed for defining the initial(water table) conditions across the two-dimensional finite element grid for numerical simulation of the flow field at the fieldsite. Soil samples were collected at 30-cm increments up to adepth of 120 cm using Veihmeyert soil sampling tubes. Thesesamples were later used for gravimetric measurements of thesoil water content in the laboratory. Because of uncertainity
about the possible effects of freezing at the field site during thewinter of 1994/1995, we repeated the soil water content measurements along the same transect in February 1995. The
spring soil water content data were later used to define initialconditions for the flow simulations. These initial conditions
were thought to be appropriate for estimating and simulating
surface-contributed flow into the subsurface tile drains since2054

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

Figure 3. Error bar (geometric mean and standard error) plots of the ^field-averaged& soil hydraulic conductivity (K-h) function for the surface horizon by superimposing the in situ and laboratory measurements:(a) semilog plot and (b) log-log plot. Inflection points were found between 1.5 and 3 cm (1), and between 3
cm and 30 cm (2). Inflection point 1 near saturation indicates that the flow-governing process changes fromnon– capillary dominated to capillary dominated.
the regional groundwater flow contribution was relatively small
immediately after the spring irrigation events. In addition tointermittent irrigation and rainfall events, preestablished empirical relationships for the geographical region and crops[Mapel et al., 1985] were used for calculating time-dependent
evapotranspiration rates from the top soil layers near the soilatmosphere boundary.
In addition to the soil water content measurements along thetransect, the observation wells and piezometers along theboundary of the field were monitored once a month to estimate
the water table head and its gradient associated with regional
groundwater flow below the field site. Water table heights wereinterpolated to define time-dependent hydrologic boundary
conditions for the two-dimensional (or pseudo threedimensional) numerical flow domain. Slug tests were conducted in the piezometers to define the hydraulic conductivityof the aquifer at deeper depths. A thick clay lens at a depth of
7 m below our field site [NRCS, 1992] served as the lower
no-flow boundary for the numerical simulations. We believethat the assumption of an impermeable boundary is reasonable
since this study focused on short-term preferential flow phenomena in the tile drains resulting from surface irrigation
events.

5. Soil Water Retention and Hydraulic
Conductivity Functions
The hydraulic conductivity data obtained with the ponded
infiltrometer, tension infiltrometer, and laboratory outflow ex-

periments were superimposed across the tension range of 0 to
1700 mm for all measurement sites. Considering the lack ofany single instrument for measuring the hydraulic conductivity
across the entire soil water tension range from wet to dry, thesuperimposition approach of Clothier and Smettem [1990] appears to be very practical for characterizing the hydraulic conductivity function of any field soil exhibiting preferential flownear saturation. However, we emphasize that the resulting
hydraulic conductivity function is still based on the use of
average soil water tensions across the entire pore-size range forthe soils under investigation. In other words, the hydraulicconductivity function represents the bulk soil rather than specific flow domains as required in the multiporosity modelingapproach. While the latter approach may be more suitable
when different pore groups can be quantified, for example, ina controlled column experiment, the lumped method seems
more practical for addressing large-scale in situ flow problems.
Note that the disc (infiltrometer) surface area (182.4 cm2) isassumed to encompass the representative elementary volume(REV) of soil for all possible pore (flow) domains at the fieldsite. Figure 3 shows semilog and log-log plots of the fieldaverage K-h curve based on 82 measurements across the study
section of the center bench. Interestingly, the ^field-averaged&
hydraulic conductivity curve for the surface soil showed threeturning (or two inflection) points within the measured tension
range, that is, at approximately 15, 30, and 300 mm. The sharpchanges in K across the inflection point between 15 and 30 mm
suggests the presence of relatively fast (preferential) flow phe-

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

nomena near saturation. When the K-h functions were plotted
for the individual measurement sites (not shown here), we
encountered similar inflection points although their positions
on the curves varied somewhat from site to site. Also, we foundthat the standard errors of the hydraulic conductivity in thetension range of 0 (K 0 ) to 20 mm (K 20 ) were much larger incomparison with standard errors in K h at higher tensions (h)(Figure 3).
Several studies [Jarvis and Messing, 1995; van Genuchten,
1994; Wilson et al., 1992; Campbell and Shiozawa, 1992; van
Genuchten et al., 1991; Clothier and Smettem, 1990; Seyfried andRao, 1987] have inferred that because of its large variability,
the saturated hydraulic conductivity (K s or K 0 as per ourconvention) is not a good physical parameter for fitting or
scaling hydraulic functions of field soils. For example, Booltink
et al. [1991] found that K 0 varied from 10 to 5500 cm/day andthat these two values generated two distinctly different K-hfunctions near saturation. In spite of these findings, K 0 is still
popularly used as a matching point in most hydrologic studies
since this parameter can be more easily measured. This use ofK 0 can lead to extreme underpredictions or overpredictions ofrelated hydraulic parameters (e.g., a in Gardners’ model or aand n in the van Genuchten–Mualem model), thereby also
affecting the accuracy of numerical simulations using thosedata. By using tension infiltrometers for measuring K h at different soil water tensions (h) in situ, we explored other more
appropriate K h values as possible matching points for the hydraulic conductivity functions. Such an approach had been
advocated by, for example, Campbell and Shiozawa [1992] andLuckner et al. [1989] as an alternative to the traditional approach of fitting K at zero tension. Clothier and White [1981]
and Smettem and Collis-George [1985] found somewhat similar
approaches useful for determining ponding time and for infiltration modeling, respectively, although the exact position of
the matching point may remain somewhat arbitrary.
Similar to the Figure 3, Figure 4 shows measured soil watercontents ( u -h) across the soil water tension range. Although
the soil water content data near saturation are not very precise,one may still notice a very small but conceptually important
increase in u in the tension range of 30 to 0 mm. This result
confirms similar findings by Wilson et al., [1992] and Germann
and Beven [1981]. Moreover, our result qualitatively supports
Hoogmoed and Bouma [1980], who claimed that lateral absorption from macropores into the soil matrix during bypass flow isin the order of only a few percent. Figure 4 (field-averaged soil
water retention) and Figure 3 (field-averaged hydraulic conductivity) over the same (30- to 0-mm) tension range, indicatethat while the hydraulic conductivity showed a sevenfold (or550%) increase, water retention increased only very little. Thisis consistent with similar findings by R. W. Skaggs (personalcommunication, 1995), Clothier and Smettem [1990] and Seyfried and Rao [1987]. We believe that the soil water retentionfunction of the macroporous soil combines several features ofboth the rigid soil matrix and swell-shrink macropores as was
suggested by Klute [1986]. Mohanty et al. [1996], Booltink andBouma [1993], and Watson and Luxmoore [1986] found forsimilar field studies that approximately 80 –90% of total infiltration occurs close to saturation between tensions of 0 and 30
mm. All of the above findings seem to indicate that in thisnarrow soil water tension range, water is mostly gravity driven
and partially mobile and that this water can increase the flux
many fold. At the same time, total water retained in the bulk
soil appears to remain relatively constant in this range. The

2055

Figure 4. Error bar (arithmetic mean and standard error)
plot of the ^field-averaged& soil water retention ( u -h) function
for the surface horizon measured in the laboratory. Turning
point at 30 cm indicates the bubbling pressure.
very small increase from ( u h )* to (u0) (in Figure 4) may be
only an experimental artifact, which, if true, suggests that watercontent at this tension range could be equated to ( u h )* for allpractical purposes, even for macroporous field soils. However,
because only one data point between ( u h )* and (u0) was available, a reasonable statistical analysis could not be performed toverify our claim. This situation could be approximated numerically by imposing a slightly higher value at 0 tension than at h*and using linear interpolation between h* and 0.
Closed-form van Genuchten–Mualem type models [van Genuchten, 1980] for the soil water retention and hydraulic conductivity have been quite successful in describing variably saturated flow as long as the soil has a relatively unimodal poresize distribution [e.g., van Genuchten and Nielsen, 1985]. Whenthe soil has a multimodal pore-size distribution, these models
become inaccurate close to saturation [Durner, 1994, and references therein]. One reason for this failure is a disruption in
the otherwise smooth pore-size distribution of a macroporous
soil at or near the larger pores. Also, gravity flow through thelarger pores during near-saturated conditions will dominate
capillary flow which was the basis of Mualem’s model. Morerecently, several researchers [Mallants et al., 1995; Durner,
1994; Othmer et al., 1991; Peters and Klavetter, 1988] haveachieved some success in modeling the hydraulic properties ofbimodal or multimodal soils using sums of a number of van
Genuchten–Mualem type functions. They found that these
multimodal hydraulic functions did perform better than unimodal hydraulic functions at intermediate soil water suctions( pF 1 to 5). However, the sum functions generally did a betterjob only in optimizing either the retention or the hydraulicconductivity function near saturation, but not both [Durner,
1994]. Another group of researchers [Jarvis and Messing, 1995;
Messing and Jarvis, 1993; Wilson et al., 1992] used a number of2056

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

piecewise-continuous Russo-Gardner type [Russo, 1988] functions, joined at different junction points, for describing theretention and hydraulic conductivity of different pore groups.
Although these exponential models (for the hydraulic conductivity) could somewhat better fit the actual data close to saturation [Messing and Jarvis, 1993], they became less accurate at
midrange soil water suctions. These findings for sum and junction type functions suggest that the soil water retention andhydraulic conductivity functions of multimodal field pore systems should be treated separately [Logsdon, 1995]. Moreover,both the sum and junction approaches have their related mathematical and numerical advantages and limitations [Rossi and
Nimmo, 1994]. Thus, for field soils having multimodal pore size
distributions, it may be more appropriate to invoke a hybrid
sum-junction approach involving two or more piecewisecontinuous functions based on the dominant physical (gravity,
capillarity, or adsorption) forces associated with different conducting pore regions across the soil water tension range.
The above findings by different researchers, as well as ourown measurements, suggest that the van Genuchten retentionfunction in conjunction with Mualem’s hydraulic conductivity
function should be used only with care in a small but criticaltension range near saturation for macroporous field soils. In
other words, the van Genuchten–Mualem models should be
used only for tensions greater than h* (e.g., 3 cm for our casestudy) when capillary flow dominates, in conjunction withother more appropriate models for the soil water retention andespecially the hydraulic conductivity curve for tensions lessthan h* when gravity flow dominates. The threshold tension
h* may well change from soil to soil. Nonlinear optimization[e.g., Wilson et al., 1992] could be used to filter out the inflection points in the soil water retention and hydraulic conductivity curves and hence to distinguish between the different
flow regions. In this study, however, visual examination ofthese functions was assumed to be sufficient, and h *u (for thesoil water retention function) was approximated with h *K (forthe hydraulic conductivity function). To avoid confusion whendefining our new hydraulic functions at 0 and h* tensions, wewill henceforward use K and u with the related tension values
as subscripts.
Hybrid sum-junction hydraulic functions (1)–(5) follow. Fora multimodal pore-size distribution the new soil water retention and hydraulic conductivity functions can be written as

u ~h! 5O

w cu c~h! 5cO F

w c u r,c 1cu s,c 2 u r,c@1 1 ~ a ch! nc# mc

G(1)h # h *u

K~h! 5u 5 u s~c51!

O

w cK c~h!

Owch . h *u

(2)

c

5

c

~1 2 ~ a c h! nc21@1 1 ~ a ch! nc# 2mc! 2@1 1 ~ a ch! nc# mc/ 2~m c 5 1 2 1/n c!

K nc~h! 5 K* 1 K*@exp~h2h*!d 2 1#
K nc~h! 5 K* 1 K*@exp~2h*!d 2 1#
where

(3)
h # h *Kh *K , h # 0

(4)

h.0

(5)

K c the hydraulic conductivity for capillarydominated flow domain c [L T 21 ];K nc the hydraulic conductivity for
non– capillary dominated flow domain nc[L T 21 ];u s,c the saturated water content for capillarydominated flow domain c [L 3 L 23 ];u r,c the residual water content for capillarydominated flow domain c [L 3 L 23 ];h the equilibrium soil water pressure headfor bulk soil (across all flow domains) [L];h *K ' h *u 5 h* the critical or break-point soil waterpressure head where flow changes fromcapillary dominated to non– capillarydominated flow or vice versa [L];K* the hydraulic conductivity corresponding
to h* [L T 21 ];d a fitting parameter representing effective
macroporosity or other structural featurescontributing to non– capillary dominatedflow [L 21 ];a c , n c the van Genuchten fitting parameters[van Genuchten, 1980] for the capillarydominated flow domain c [L 21 , –],c the number of capillary dominated flow
domains where for c 5 1 the sum typemultimodal van Genuchten–Mualem
hydraulic functions (equations (1) and(3)) reduce to the original unimodal van
Genuchten–Mualem–type functions;nc the non– capillary dominated flow domain;w c the weighting factor for capillarydominated flow domain c[–], subjected to• w c 5 1 and 0 , w c , 1.
When the soil is only bimodal (as shown by K-h measurements within the adopted tension range for our field site;
Figure 3) and hence has only one capillary dominated and one
non– capillary dominated flow domain, the hydraulic conductivity functions reduce to forms that closely resemble the Fermi
function used by Wilson et al. [1992] and the bimodal function
proposed by Smettem and Kirkby [1990]. However, the important difference between our approach and the Fermi function[Wilson et al., 1992, equation 9] is that our K 0 is a floating
parameter and depends on the (measured) K* and (fitted) d,
while Wilson et al. [1992] fixed their K s and K* to match theobserved K- u data in the macropore region. In our study, as
discussed earlier and unlike Wilson et al. [1992], we assumed
that u is partially mobile and gravity driven and that the change
in u near saturation may possibly be an (laboratory) experimental artifact, so that we could keep u* constant beyond h*for our field soils. We also extrapolated K (in the tension rangeof less than h*) from K* using the nonlinear equation (4). In
contrast, Wilson et al. [1992] adopted a linear interpolation
technique between the two fixed points in their macropore
region.
If no distinct changes in hydraulic conductivity were observed near saturation, we optimized our hydraulic parameters
for the surface and subsurface horizons using the nonlinear
optimization schemes of RETC (for steady state flow experiments [van Genuchten et al., 1991]) and SFIT (for transient
flow experiments [Kool and Parker, 1987]). For situations withbimodal hydraulic characteristics, the breakpoint soil waterMOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

2057

Table 3. ^Field-Averaged& Parameters of Bimodal Piecewise-Continuous and Unimodal van Genuchten–Mualem Hydraulic
Functions for Different Soil Horizons at the Field Site in Las Nutrias, New Mexico
Bimodal Piecewise-Continuous Hydraulic Functions
Depth,cm^ u s,i51 &,cm3 cm23^ u r,i51 &,cm3 cm23^ a &,
cm21

^n&
(–)

^K*&,
cm/s

^h*&,
cm

^ d &,
cm21

0–40
40–100*100–700†

0.475
0.459
0.430

0.110
0.050
0.090

0.015
0.021
0.083

1.60
1.40
2.12

5.55E204
4.16E204
2.73E203

3
3
0

0.92
0.70
n/a

Unimodal van Genuchten-Mualem Hydraulic Functions (High/Low)
Depth
(cm)

^ u s &,
cm3 cm23^ u r &,
cm3 cm23^ a &,
cm21

^n&
(–)

^K s &,
cm/s

0–40
40–100100–700†

0.484/0.475
0.464/0.459
0.430/0.430

0.080/0.120
0.045/0.050
0.090/0.090

0.004/0.015
0.010/0.021
0.083/0.083

1.44/1.55
1.25/1.38
2.12/2.12

4.08E203/9.50E204
3.50E203/6.80E204
2.73E203/2.73E203

The ^field-averaged& parameters were obtained by fitting the van Genuchten–Mualem hydraulic functions (equations (1) and (3)) or the
inverted exponential function (equation (4)) to (geometric or arithmetic) mean hydraulic data (e.g., as shown in Figures 4 and 5).
*As hydraulic parameters near saturation were not precise (using IP experiment) at deeper depths, ^h*&, ^ d & parameters were adopted andcalibrated from the surface horizon.
†As soil type varied from fine sand to sand in this depth range, flow was assumed to be well presented by unimodal van Genuchten–Mualem
model.

pressure head (h*) was determined by visual inspection, followed by separate nonlinear optimization of the hydraulic parameters for the capillary and noncapillary domains. Table 3
shows average values of the different hydraulic parametersoptimized from field-averaged hydraulic functions (Figures 3
and 4) at different depths for our field site. Since the adopted
IP experiment did not permit precise measurement of thesubsurface hydraulic properties very close to saturation (asopposed to using the disc infiltrometer), we adjusted the hydraulic conductivity functions near saturation for the C1 (silty
clay loam–silty clay) horizon by model calibration.

6.

Numerical Modeling

The new hybrid sum-junction-type hydraulic functions wereincorporated into the CHAIN_2D computer model [Simunek
and van Genuchten, 1994] to enable two-dimensional transient
flow simulations at the field site. We assumed that water flow,
overall, could be described with Richards’ equation (for Darcytype flow under isothermal and variably saturated conditions)in conjunction with our bimodal (or multimodal) piecewisecontinuous hydraulic functions. The soil was assumed to be
rigid and nonhysteretic. For two-dimensional numerical simulations, the north-south cross section across the center bench
was divided into 9350 nodes and 9126 quadrilateral elements(as shown in Figure 5) by subdividing the soil sampling andsurface infiltration measurement grid (Figure 2). On the basisof the surface infiltration measurements and the IP experimental results, the soil profile was divided into three soil horizons(characterized by three ^field-averaged& hydraulic functions).
These three horizons were further subdivided into a total of
169 layers for numerical simulation purposes. Finer discretizations were used near the soil surface, across the horizon interfaces, and around the subsurface tile drain to handle intermittent (abrupt) changes in local fluxes and hence the pressure
gradients. In addition, we kept our maximum vertical discretization relatively small (Dz 5 1 cm) in the upper two horizons(of depth of 100 cm). Finer discretization improved the performance of numerical solutions, in particular when piecewise-

continuous hydraulic functions were used in upper two horizons (i.e., Dz , uh*u). For the simulations we used a timevariant evapotranspiration rate assuming a normalized water
uptake distribution across the top 40 cm of the soil profile inaccordance with the approach of Feddes et al. [1978]. Measured
soil water contents before the crop growing season (i.e., February 1994 or 1995) at different nodal points and interpolated
water table contours from the groundwater monitoring wells
along the field boundary were used to define the initial conditions in the two-dimensional flow domain. Time-dependent
boundary conditions, such as intermittent irrigation/precipitation rates at atmospheric boundary nodes; groundwater headvalues along hydrologic boundary nodes; and no-flow conditions along impermeable bottom boundary nodes were invoked. Since gradients of the groundwater table (Table 4) were
found to be very small, we simplified the two hydrologic boundaries G A D and G B C , Figure 5) by assuming impermeableboundaries. The tile drain at the center of the study section was
treated as a boundary node surrounded by four regular square
elements with adjusted hydraulic conductivities using the electric analog approach of Vimoke et al. [1963] and Fipps et al.
[1986].
The initial value problem was solved assuming variably saturated flow [Richards, 1931], with the initial distribution of thepressure head within the flow domain, V, being a linear function of depth, and implementing system-independent, timedependent specified fluxes including intermittent irrigation/
rainfall and variable evapotranspiration rates (as shown in
Figure 5). In addition to system-independent boundary conditions, CHAIN_2D considers system-dependent boundary conditions that cannot be defined a priori. For example, the potential fluid flux across the soil-air interface is controlled
exclusively by external conditions (e.g., irrigation/rainfall).
However, the actual flux depends also on the prevailing (transient) soil moisture conditions. Thus the soil surface boundary
may change from prescribed flux to prescribed head-type conditions (and vice versa). The governing equation and associated initial and boundary conditions are hence as follows.

2058

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

Figure 5. Typical/adapted finite element mesh including system-independent (time-dependent) and systemdependent flux/head-type atmospheric boundary, no-flux bottom and hydrologic boundaries, tile drain boundaries, and root water uptake sink term for the two-dimensional flow simulation in the vertical cross-section of
the field (study section) at Las Nutrias, New Mexico. The cross section of unit length is assumed to be
representative over the 240-m-long study section.

Governing Equation
­u
­
5
­t
­ xi

FS

K K ijA­h
1 K izA
­ xj

DG

2S

(6)

where S is a sink term [T 21 ] for root water uptake, x i (i 5 1,
2) are spatial coordinates [L] (i.e., x 1 5 x, and x 2 5 z), t is
time [T], and K A
ij are components of a dimensionless anisotropic conductivity tensor, K A . Because no horizontal K measurements at deeper depths were made at the field site, we
assumed K A to be isotropic in this study. Other terms in (6)were defined before.
Initial Conditions
h~ x, z, t! 5 h 0~ x, z!

t50[V

(7)

where s[–] are prescribed (irrigation/rainfall water flux) functions of x, z, and t; and n i are the components of the outward
unit vector normal to the boundary G AB ; E is the maximum
potential rate of infiltration or evaporation under the current
atmospheric conditions; h is the pressure head at the soil
surface; and h A and h s are, respectively, minimum and maximum pressure heads allowed under the prevailing soil conditions. When h s was set to zero, E and h A were calculated using
relationships defined in Feddes et al. [1974].
Hydrologic Boundary G AD and G BC and Lower Impermeable
Boundary G CDFS

2 K K ijA

where h0 is a prescribed function of x and z in the flow region, V.
Boundary ConditionsFS­h
2 K K ijA
1 K izA­ xj

DG

n i 5 s 1~ x, z, t!
~ x, z! [ G AB

(8)t [ t flux

System-Dependent Atmospheric Boundary G AB

S

DG

ni 5 0(11)~ x, z! [ G AD ~ G BC ~ G CD;t

K drain 5 KC d

(12)

Tile Drain Boundary

System-Independent Atmospheric Boundary G AB

uK K ijA­h
1 K izA­ xj

D­h
1 K izA n iu # E­ xj~ x, z! [ G AB

(9)

and
hA # h # hs

~ x, z! [ G AB

(10)

where K drain is the adjusted conductivity [L T 21 ] and C d is acorrection factor [–] to be calculated from the ratio of theeffective drain diameter, d e [L], of the drain and the side
length, D [L], of the square formed by finite elements surrounding the drain node [Vimoke and Taylor, 1962].
Since this study focused on modeling preferential flow thatwould occur during or relatively soon after a large water input
event (e.g., flood irrigation), we used hourly time steps andlimited the real simulation time to approximately 100 hours
from the start of an irrigation event. Also we adopted 0.1 hour
(6 min) initial time steps for our numerical simulations toMOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

closely follow the 5-min intervals of flow measurements in thefield. For the numerical simulations we used a du convergencecriterion [Huang et al., 1996] in conjunction with the “massconservative” modified Picard iteration method [Celia et al.,1990] to enable relatively fast and robust convergence of thenumerical solution. Other details of the numerical procedures
adapted in CHAIN_2D are given by Simunek and van Genuchten [1994].

7.

Simulation Results

Simulated tile flow using the new ^field-averaged& bimodal
piecewise-continuous soil water retention and hydraulic conductivity functions (Table 3) are presented in Figure 6 andcompared against the field-observed values. We also compared
the simulation results with results obtained using ^fieldaveraged& unimodal van Genuchten–Mualem–type soil water
retention and hydraulic conductivity functions. To present aqualitative comparison, we used two unimodal functions in thesurface horizon based on ^K 0 mm& through ^K 30 mm& being
included in the parameter optimization (unimodal_high), and^K 0 mm& through ^K 30 mm& being excluded from the parameter
optimization (unimodal_low). The short-term simulation runs
were made separately by setting up the initial conditions andboundary conditions at the start of each simulation based onfield measured data. Our simulations considered four flow
scenarios: (1) the west bench irrigated prior to the center
bench, (2) the center bench irrigated prior to the west bench,(3) irrigation during periods with a relatively low groundwatertable (i.e., initially below the tile line), and (4) irrigation with arelatively high groundwater table (i.e., initially above the tileline). From the observed tile-drain flow it is clear, that following an irrigation event, water moves rapidly downward throughthe soil profile. Depending on initial conditions, and water
input amount and intensity, different irrigation events required
10 –20 hours for the water to reach the water table and subsequently the tile line. After a few hours of irrigation the tile flowrate suddenly increased and reached a peak that was almost
instantly followed by a similar sharp drop, suggesting preferential flow through the field macropore network. This situation
was followed by a recession phase dominated by prolonged
matrix flow and water redistribution in the profile. In general,
predicted flow using the bimodal hydraulic functions matchedthe observed flow rates reasonably well and far better thanwhen the unimodal hydraulic functions were used.
A recent hypothetical study by Zurmuhl and Durner [1996]
gave useful insight into the hydraulics of a similar type ofbimodal porous medium. They found that transient soil waterpressure head distributions versus depth changed far more
abruptly for a bimodal soil than for a unimodal soil; our experimental and modeling results qualitatively agree with theirfindings. Our new hydraulic functions appear especially usefulfor flood-irrigation events when the water potential in the bulk
soil can be represented more or less by a single value (as in unit
gradient condition), unlike the multidomain or multiporosity
approach. For this same reason the method may lose its advantage for lower-rate water application events. Furthermore,
hysteretic, nonisothermal unsaturated flow, and threedimensional groundwater flow are some of the missing characteristics of our modeling approach; these processes likely
also contributed to the significant differences between observed and simulated flux values.
A comparison of the simulation results (Figure 6) for the2059

Table 4. Mean Hydraulic Head and Hydraulic Gradients in
Shallow Groundwater Aquifer Under the Center Bench
Field in Las Nutrias, New Mexico
Date

Mean Head,m

Mean Gradient,
m/m

Feb.
March
April
May
June
July
Aug.
Sept.
Oct.
Nov.
Dec.

1994
98.79
99.39
99.24
99.17
99.20
99.28
zzz
99.17
98.95
98.90
98.85

4.30E204
5.10E204
1.00E203
6.50E204
7.20E204
8.00E204zzz
6.30E204
8.00E204
7.20E204
7.20E204

Jan.
Feb.
March
April
May
June
July
Aug.
Sept.
Oct.

1995
98.78
98.75
98.95
zzz
99.01
99.10
99.18
99.26
99.20
zzz

6.90E204
8.00E204
5.80E204zzz
8.30E204
8.70E204
7.40E204
5.10E204
5.80E204
zzz

Ellipses denote no measurement. Reference elevation is 100.00 m at
the soil surface.

four flow scenarios does reveal some interesting spatiotemporal phenomena governing tile flow in the experimental area.
We emphasize here that no calibration of field-measured hydraulic functions was attempted during our field-scale-flow
simulations. Peak flow predictions for scenario (1) slightly underestimated the observed peak. This is because west bench
irrigation prior to the center bench irrigation reduced thegroundwater table gradient in the east-west direction, thus
forcing the subsurface flow field (with larger capture zone) toconverge earlier to the tile drain. On the other hand, the higher
east-west gradient in the groundwater table for scenario (2)caused a relatively delayed lower peak in the observed tile flow
as compared to the model predictions. Furthermore, for scenario (3) tile flow showed atypical negative flow afterwards.
We believe that this may have been an artifact of temporarily
reversed flow from the downgradient (west) bench to the upgradient (center) bench; this could not be registered by the
one-way tile-flow valve located in the west sump. A detailed
discussion of some of these unusual flow processes near the
manholes, and their governing factors, is given by B. P. Mohanty et al. (Preferential transport of nitrate to a tile drain in
an intermittent-flood-irrigated field: Model development andexperimental evaluation, submitted to Water Resources Research, 1997). Comparison of tile flow predictions for scenarios(3) and (4) shows better model performance under low
groundwater table conditions at the field site. Although theexact reason for this behavior is still unknown, we believe thatfor the initially high groundwater table conditions, regional
groundwater contribution to the tile flow may not have been
properly accounted for using the adopted hydrologic boundary
conditions in the two-dimensional model. A three-dimensional
flow model may be a more suitable for our field site. We
suspect also that temporal variability in the soil hydraulic prop-

2060

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

Figure 6. Comparison of simulated tile-drain flow (from the study section) using the piecewise-continuous
hydraulic functions (bimodal soil porosity distribution) and the van Genuchten–Mualem (unimodal soil
porosity distribution) hydraulic functions with observed data (net tile drain flow below the center bench) at thefield site. (a) Scenario 1, west bench irrigated prior to the center bench; (b) scenario 2, center bench irrigatedprior to west bench; (c) scenario 3, irrigation during (initially) low groundwater table condition; and (d)scenario 4, irrigation during (initially) high groundwater table condition. Time 0 corresponds to the midnightof the beginning day of center-bench irrigation.

erties may have played an important role in this behavior. A
follow-up paper will address some of these issues, including
spatial variability considerations using similar media scaling
schemes.
Most importantly, from the different simulation runs weconclude that using the piecewise-continuous hydraulic functions generally minimized the overpredictions or underpredictions of flow in the different flow domains as is typically obtained with the one-domain approach. This study shows thatmatrix flow (during higher soil water tension) can be kept
almost unchanged during calibration of the macropore flow

domain (near saturation) and vice versa. We also found thatthe new piecewise-continuous hydraulic functions outperformed the unimodal function (unimodal_high with equal K s )in terms of better predicting the timing and magnitude of thepreferential flow peaks.

8.

Concluding Remarks

Capillary flow– based unimodal soil water retention and hydraulic conductivity functions should be complemented withother functions to correctly simulate noncapillary dominatedMOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS

Figure 6.

flow near saturation in field soils. In this study we measured therelative trend of soil water retention and hydraulic conductivity
functions near saturation and established piecewise-continuous functions to predict flow across a wide range of soil watertensions. The new hydraulic functions were used to predict thefield-scale flow under a flood-irrigated field in Las Nutrias,
New Mexico. A comparison of the performance of the newpiecewise-continuous (bimodal) function with that of the (unimodal) van Genuchten–Mualem soil water retention and hydraulic conductivity functions showed that the new functionsindeed performed better for the field soil. Unresolved issues of
our new approach include (1) precision in differentiating thematching points (h*, K*, u *) for different types of soil, (2)precision of measurement of u near saturation for different
types of soil, and (3) description and inclusion of site-specificmultidimensional regional hydrology in the mathematical

2061(continued)

model. To resolve some of these issues and to make the approach more robust, independent mathematical and more controlled (laboratory, field) experiments are warranted. Finally,
our study also suggests that field-scale hydrologic model development always should be corroborated with field measurements.
Acknowledgments. This project was partly funded by the USDACSRS through special water quality grant 9301322. Thanks to graduate
students Jeff Chaves, Tracy Roth, Bob Reedy, and Gunter Hilmes at
New Mexico Institute of Mining and Technology, Socorro, New Mexico, for field and laboratory data collection.

References
Anderson, S. H., and J. W. Hopmans (Eds.), Tomography of Soil Water
Root Processes, Spec. Publ. 36, Soil Sci. Soc. Am., Madison, Wis.,
1994.

2062

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONSAndreini, M. S., and T. S. Steenhuis, Preferential paths of flow underconventional and conservation tillage, Geoderma, 46, 85–102, 1990.
Ankeny, M. D., Methods and theory for unsaturated infiltration measurements, in G. Topp et al. (editors), Advances in Measurement ofSoil Physical Properties: Bringing Theory into Practice, Spec. Publ. 30,
edited by G. Topp et al., pp. 123–141, Soil Sci. Soc. Am., Madison,Wis., 1992.
Ankeny, M. D., T. C. Kaspar, and R. Horton, Design of an automatedtension infiltrometer, Soil Sci. Soc. Am. J., 52, 893– 896, 1988.
Beven, K., and P. Germann, Macropores and water flow in soils, Water
Resour. Res., 18, 1311–1325, 1982.
Booltink, H. W. G., Morphometric methods for simulation of water
flow, Ph.D. thesis, Univ. of Wagennigen, Wagennigen, Netherlands,
1993.
Booltink, H. W. G., and J. Bouma, Sensitivity analysis on processes
affecting bypass flow, Hydrol. Processes, 7, 33– 43, 1993.
Booltink, H. W. G., J. Bouma, and D. Gimenez, A suction crust
infiltrometer for measuring hydraulic conductivity of unsaturated
soil near saturation, Soil Sci. Soc. Am. J., 55, 566 –568, 1991.
Bouma, J., Soil morphology and preferential flow along macropores,
Agric. Water Manage., 3, 235–250, 1981.
Bouma, J., and L. W. Dekker, A case study on infiltration into dry clay
soil, I, Morphological observations, Geoderma, 20, 27– 40, 1978.
Bowman, R. S., and R. C. Rice, Transport of conservative tracers in
the field under intermittent flood irrigation, Water Resour. Res., 22,1531–1536, 1986.
Campbell, G. S., and S. Shiozawa, Prediction of hydraulic properties ofsoils using particle-size distribution and bulk density data, in Indirect
Methods for Estimating the Hydraulic Properties of Unsaturated Soils,
edited by M. T. van Genuchten, F. J. Leij, and L. J. Lund, pp.
317–328, Univ. of Calif., Riverside, 1992.
Celia, M. A., E. T. Bouloutas, and R. L. Zabra, A general massconservative numerical solution for the unsaturated flow equation,
Water Resour. Res., 26, 1483–1496, 1990.
Clothier, B. E., and K. R. J. Smettem, Combining laboratory and fieldmeasurements to define the hydraulic properties of soil, Soil Sci.
Soc. Am. J., 54, 299 –304, 1990.
Clothier, B. E., and I. White, Measurement of sorptivity and soil water
diffusivity in the field soil, Soil Sci. Soc. Am. J., 45, 241–245, 1981.
Corey, A. T., Pore-size distribution, in Indirect Methods for Estimating
the Hydraulic Properties of Unsaturated Soils, edited by M. T. van
Genuchten, F. J. Leij and L. J. Lund, pp. 37– 44, Univ. of Calif.,
Riverside, 1992.
De Vries, J., and T. L. Chow, Hydrologic behavior of a forested
mountain soil in coastal British Columbia. Water Resour. Res., 14,935–942, 1978.
Durner, W., Hydraulic conductivity estimation for soils with heterogenous pore structure, Water Resour. Res., 30, 211–223, 1994.
Edwards, W. M., R. R. van der Ploeg, and W. Ehlers, A numerical
study of the effects of non-capillary sized pores upon infiltration, Soil
Sci. Soc. Am. J., 43, 851– 856, 1979.
Edwards, W. M., L. D. Norton, and C. E. Remond, Characterizing
macropores that affect infiltration into nontilled soil, Soil Sci. Soc.
Am. J., 52, 483– 487, 1988.
Ela, S. D., S. C. Gupta, and W. J. Rawls, Macropore and surface sealinteractions affecting water infiltration into soil, Soil Sci. Soc. Am. J.,
56, 714 –721, 1992.
Everts, C. J., and R. S. Kanwar, Estimating preferential flow to a
subsurface drain with tracers, Trans. ASAE, 33, 451– 457, 1990.
Everts, C. J., and R. S. Kanwar, Interpreting tension-infiltrometer datafor quantifying soil macropores: Some practical considerations,
Trans. ASAE, 36, 423– 428, 1993.
Feddes, R. A., E. Bresler, and S. P. Neuman, Field test of a modified
numerical model for water uptake by root systems, Water Resour.
Res., 10, 1199 –1206, 1974.
Feddes, R. A., P. J. Kowalik, and H. Zaradny, Simulation of Field Water
Use and Crop Yield, John Wiley, New York, 1978.
Fipps, G., R. W. Skaggs, and J. L. Nieber, Drains as a boundary
condition in finite elements, Water Resour. Res., 22, 1613–1621, 1986.
Gardner, W. R., Calculation of capillary conductivity from pressure
plate outflow data, Soil Sci. Soc. Am. Proc., 20, 317–320, 1956.
Gerke, H. H., and M. T. van Genuchten, A dual-porosity model for
simulating the preferential movement of water and solutes in structured porous media, Water Resour. Res., 29, 305–319, 1993a.
Gerke, H. H., and M. T. van Genuchten, Evaluation of a first-orderwater transfer term of variably saturated dual-porosity flow models,
Water Resour. Res., 29, 1225–1238, 1993b.
Germann, P. F., and K. Beven, Water flow in soil macropores, I, An
experimental approach, J. Soil Sci., 32, 1–13, 1981.
Gwo, J. P., P. M. Jardine, G. V. Wilson, and G. T. Yeh, Using a
multiregion model to study the effects of advective and diffusive
mass transfer on local physical nonequilibrium and solute mobility ina structured soil, Water Resour. Res., 32, 561–570, 1996.
Hoogmoed, W. B., and J. Bouma, A simulation model for predicting
infiltration into a cracked clay soil, Soil Sci. Soc. Am. J., 44, 458 – 461,
1980.
Huang, K., B. P. Mohanty, and M. T. van Genuchten, A new convergence criterion for the modified picard iteration method to solve thevariably saturated flow equation, J. Hydrol., 178, 69 –91, 1996.
Hutson, J. L., and R. J. Wagenet, A multiregion model describing
water flow and solute transport in heterogenous soils, Soil Sci. Soc.
Am. J., 59, 743–751, 1995.
Jardine, P. M., G. V. Wilson, R. J. Luxmoore, and J. F. McCarthy,
Transport of inorganic and natural organic tracers through an isolated pedon in a forested watersheds, Soil Sci. Soc. Am. J., 53,
317–323, 1989.
Jardine, P. M., G. V. Wilson, and R. J. Luxmoore, Unsaturated solutetransport through a forest soil during rain storm events, Geoderma,
46, 103–118, 1990.
Jarvis, N. J., and I. Messing, Near-saturated hydraulic conductivity insoils of contrasting texture measured by tension infiltrometers, Soil
Sci. Soc. Am. J., 59, 27–34, 1995.
Klute, A., Water retention: Laboratory methods, in Methods of Soil
Analysis, part I, Agron. Monogr., vol. 9, edited by A. Klute, pp.
635– 662, Am. Soc. of Agron., Madison, Wis., 1986.
Kool, J. B., and J. C. Parker, Estimating soil hydraulic properties fromtransient flow experiments: SFIT user’s guide, Electr. Power Res.
Inst., Palo Alto, Calif., 1987.
Logsdon, S. D., Flow mechanisms through continuous and buried
macropores, Soil Sci., 160, 237–241, 1995.
Luckner, L., M. T. van Genuchten, and D. R. Nielsen, A consistent setof parametric models for the two-phase flow of immiscible fluids inthe subsurface, Water Resour. Res., 25, 2187–2193, 1989.
Luxmoore, R. J., Micro-, meso-, and macro-porosity of soil, Soil Sci.
Soc. Am. J., 45, 671, 1981.
Luxmoore, R. J., P. M. Jardine, G. V. Wilson, J. R. Jones, and L. W.
Zelanzy, Physical and chemical controls of preferred path flow
through a forested hill slope, Geoderma, 46, 139 –154, 1990.
Mallants, D., D. J. Kim, J. Feyen, and P.-H. Tseng, Prediction of
unsaturated hydraulic conductivity using multimodal hydraulic functions, Eos Trans. AGU, 76(46), Fall Meet. Suppl., F229, 1995.
Mapel, C. L., T. W. Sammis, and R. R. Lansford, Estimating cropwater production functions based on transpiration and crop growth
curves through modeling, Rep. 191, N. M. Water Resour. Res. Inst.,
Las Cruces, 1985.
Messing, I., and N. J. Jarvis, Temporal variation in the hydraulicconductivity of a tilled clay soils as measured by tension infiltrometers, J. Soil Sci., 44, 11–24, 1993.
Mohanty, B. P., M. D. Ankeny, R. Horton, and R. S. Kanwar, Spatial
analysis of hydraulic conductivity measured using disc infiltrometers,
Water Resour. Res., 30, 2489 –2498, 1994.
Mohanty, B. P., R. Horton, and M. D. Ankeny, Infiltration and macroporosity under a row crop agricultural field in a glacial till soil, Soil
Sci., 161, 205–213, 1996.
Natural Resource Conservation Service, Soil survey map of Las Nutrias, New Mexico, Lincoln, Nebr., 1992.
Othmer, H., B. Diekkruger, and M. Kutilek, Bimodal porosity andunsaturated hydraulic conductivity, Soil Sci., 152, 139 –150, 1991.
Peters, R. R., and E. A. Klavetter, A continuum model for water
movement in an unsaturated fractured rock mass, Water Resour.
Res., 24, 416 – 430, 1988.
Phillips, R. E., V. L. Quisenberry, J. M. Zeleznik, and G. H. Dunn,
Mechanism of water entry into simulated macropores, Soil Sci. Soc.
Am. J., 53, 1629 –1635, 1989.
Prieksat, M. A., M. D. Ankeny, and T. C. Kaspar, Design an automated, self-regulating, single-ring infiltrometer, Soil Sci. Soc. Am. J.,
56, 1409 –1411, 1992.
Richards, L. A., Capillary conduction of liquids in porous mediums,
Physics, 1, 318 –333, 1931.
Richards, L. A., Physical condition of water in soil, in Methods of Soil

MOHANTY ET AL.: NEW PIECEWISE-CONTINUOUS HYDRAULIC FUNCTIONS
Analysis, Part 1, Agron. Monogr., vol. 9, edited by C. A. Black et al.,pp. 128 –152, Am. Soc. of Agron., Madison, Wis., 1965.
Rossi, C., and J. R. Nimmo, Modeling of soil water retention fromsaturation to oven dryness, Water Resour. Res., 30, 701–708, 1994.
Russo, D., Determining soil hydraulic properties by parameter estimation: On the selection of a model for the hydraulic properties, Water
Resour. Res., 24, 453– 459, 1988.
Samani, Z., Measuring water in trapezoidal canals, J. Irrig. Drain. Eng.,
119, 181–186, 1993.
Seyfried, M. S., and P. S. C. Rao, Solute transport in undisturbed
columns of an aggregated tropical soil: Preferential flow effects, Soil
Sci. Soc. Am. J., 51, 1434 –1444, 1987.
Simunek, J., and M. T. van Genuchten, The CHAIN_2D code forsimulating the two-dimensional movement of water, heat, and multiple solutes in variably-saturated porous media, U.S. Salinity Lab.,
Res. Rep. 136, 1994.
Singh, P., R. S. Kanwar, and M. L. Thompson, Macropore characterization for two tillage systems using resin-impregnation technique,
Soil Sci. Soc. Am. J., 55, 1674 –1679, 1991a.
Singh, P., R. S. Kanwar, and M. L. Thompson, Measurement andcharacterization of macropores by using AUTOCAD and automatic
image analysis, J. Environ. Qual., 20, 289 –294, 1991b.
Smettem, K. R. J., and N. Collis-George, The influence of cylindricalmacropores on steady-state infiltration in a soil under pasture, J.
Hydrol., 79, 107–114, 1985.
Smettem, K. R. J., and C. Kirkby, Measuring the hydraulic propertiesof a stable aggregated soil, J. Hydrol., 117, 1–13, 1990.
van Genuchten, M. T., A closed-form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,892– 898, 1980.
van Genuchten, M. T., New issues and challenges in soil physics research, Transactions of 15th World Congress of Soil Sci., vol. 1, 5–27;
Int. Soc. of Soil Sci., Acapulco, Mexico, 1994.
van Genuchten, M. T., and D. R. Nielsen, On describing and predicting the hydraulic properties of unsaturated soils, Ann. Geophys., 3,615– 628, 1985.
van Genuchten, M. T., F. J. Leij, and S. R. Yates, The RETC code forquantifying the hydraulic functions of unsaturated soils, EPA/600/291/065, 85 pp., Environ. Prot. Agency, Ada, Okla., 1991.

2063

Vimoke, B. S., and G. S. Taylor, Simulating water flow in soil with anelectric resistance network, Res. Rep. 41-65, Soil and Water Conserv.
Res. Div., USDA-ARS, Columbus, Ohio, 1962.
Vimoke, B. S., T. D. Tura, T. J. Thiel, and G. S. Taylor, Improvements
in construction and of resistance networks for studying drainage
problems, Soil Sci. Soc. Am. Proc., 26, 203–207, 1963.
Warner, G. S., J. L. Nieber, I. D. Moore, and R. A. Geise, Characterizing macropores in soil by computed tomography, Soil Sci. Soc.
Am. J., 53, 653– 660, 1989.
Watson, K. K., An instantaneous profile method for determining thehydraulic conductivity of unsaturated porous materials, Water
Resour. Res., 2, 709 –715, 1966.
Watson, K. W., and R. J. Luxmoore, Estimating macroporosity in aforest watershed by use of a tension infiltrometer, Soil Sci. Soc.
Am. J., 50, 578 –582, 1986.
White, R. E., The influence of macropores on the transport of dissolved and suspended matter through soil, in Advances in Soil Science, vol. 3, edited by B. A. Stewart, pp. 95–121, Springer-Verlag,
New York, 1985.
Wilson, G. V., P. M. Jardine, and J. P. Gwo, Modeling the hydraulicproperties of a multiregion soil, Soil Sci. Soc. Am. J., 56, 1731–1737,
1992.
Wilson, G. V., P. M. Jardine, R. J. Luxmoore, and J. R. Jones, Hydrology of a forested hillslope during storm events, Geoderma, 46,
119 –138, 1990.
Zurmuhl, T., and W. Durner, Modeling transient water and solutetransport in a biporous soil, Water Resour. Res., 32, 819 – 829, 1996.
R. S. Bowman and J. M. H. Hendrickx, Department of Geoscience,
New Mexico Institute of Mining and Technology, Socorro, NM 87801.
B. P. Mohanty and M. T. van Genuchten, U.S. Salinity Laboratory,
450 W. Big Springs Rd., Riverside, CA 92507. (e-mail:bmohanty@ussl.ars.usda.gov)(Received October 4, 1996; revised May 10, 1997;
accepted June 9, 1997.)