New Issues and Challenges in Soil Physics Research
M. Th. van Genuchten

Volume 1: Inaugural and State of the Art Conferences
pp. 5-27

Transactions, 15th World Congress of Soil Science, Acapulco, Mexico
July 10-16, 1994

Published by the International Society of Soil Scienceand the Mexican Society of Soil ScienceNew Issues and Challenges in Soil Physics Research
M. Th. van Genuchten. U.S. Salinity Laboratory, USDA, ARS, 4500 Glenwood
Drive, Riverside, CA, 92501, USA.
Abstract. The focus of soil physics research, and that of soil science in general, has graduallybroadened from mostly agricultural production issues to more comprehensive studies ofsubsurface water flow and chemical transport geared toward environmental issues. There is anabundance of ongoing and new challenges in soil physics research, notably in the areas ofsubsurface heterogeneity (including preferential flow), multiphase flow, multicomponent
transport, and hydraulic property characterization. Soil heterogeneity at a variety of spatial
scales remains a frustrating problem affecting the measurement, prediction, and managementof subsurface flow and transport processes. More realistic deterministic and stochastic models,
as well as improved sampling methods and instrumental techniques, are slowly becoming
available to make the research more effective. Several examples are given to illustrate soil
physics contributions to the solution of agricultural and environmental problems.
Introduction
Much evidence exists of increased human-induced stresses on the earth’s environment. Manyof these stresses are, directly or indirectly, a consequence of pollutants generated byagricultural, industrial, and municipal activities. Pesticides and fertilizers used in farming
operations continue to contaminate soil and groundwater systems worldwide. Soil erosion
remains a problem in many areas. The accumulation of salts and toxic elements in irrigation
agriculture is a worsening problem in many arid and semi-arid regions of the world. Acid rain
is a cause of pollution in most or all continents. Chemicals migrating from municipal andindustrial disposal sites are similarly creating hazards to the environment. Especiallychallenging are nonaqueous phase liquids including chlorinated hydrocarbons, polychlorinatedbiphenyls and related toxic substances. Radionuclides released into the environment (including
those produced by the Chernobyl disaster) are also increasingly becoming a problem. The
much-debated topics of deforestation, biodiversity, and climatic change, are real and must beof concern to all of us.
Contamination of the subsurface poses special challenges to soil scientists including soil
physicists. Many or most subsurface pollution problems stem directly from activities involving
the unsaturated (vadose) zone between the soil surface and the groundwater table. It is theunsaturated zone where management will offer the best opportunities for preventing or limiting
pollution, or for remediating ongoing pollution problems. Solute residence times ingroundwater aquifers can range from several years to millennia or more. Hence, once
contaminants have entered groundwater, pollution is in many cases essentially irreversible, or
can be remediated only with extreme costs. Consequently, soil and groundwater pollution
problems are best dealt with, or prevented, through proper management of the unsaturated
zone. Because of our expertise in the measurement and modeling of unsaturated zone flow andtransport processes, soil physicists are in a unique position, and have special responsibilities, tocontribute to the prevention or remediation of subsurface pollution problems.

This paper addresses some of the opportunities for soil physicists in the continuing effort to
optimize agricultural production, while also contributing to the solution of environmental
problems. It is critical that our research be part of comprehensive efforts involving other soil
scientists, hydrologists, geologists, civil and environmental engineers, and (micro)climatologists.
In this paper I shall concentrate mostly on flow and transport processes in the unsaturated zone.
Areas of research where progress has been obtained, or where progress is particularly needed,
will be emphasized. These areas include, but are not limited to, flow through deep unsaturated
soil profiles, preferential flow, flow and transport in unsaturated fractured rock, multiphase flow,
better methods and instrumentation for measuring model parameters, the need to improve
communication between soil scientists globally, and the educational implications of a broader
soil physics perspective.
Soil, Physics, and Soil Physics

Soil physics may be defined as the application of the principles of physics to the characterizationof soil properties and the understanding of soil processes, especially those involving thetransport of matter or energy (Sposito and Reginato, 1992). This definition implies that soil
physics is a subdiscipline of both physics and soil science. While useful, the definition does notnecessarily reflect the type of research typically carried out by soil physicists. It may be
instructive to briefly consider typical research topics dealt with by soil physicists over the past
several decades.
Early studies of soil physics from the 1900s to the 1940s generally involved such issues as soil
structure and soil aggregation, soil pore space, field soil water status, capillarity and soil water
retention, evaporation, soil mechanics, soluble salts, soil salinity, diffusion, and heat content.
Many of these early studies focused on the physical properties of soils as a medium for crop
production, including the effects of tillage and compaction on root growth and plant water
uptake. More comprehensive theoretical studies by soil physicists on water, heat and solute
movement in soils started in the 1930s with the formulation of the Richards equation fortransient unsaturated- water flow (Richards, 1931). Similar differential equations wereintroduced for heat and solute movement, and subsequently applied in the 1950s and 1960s to
a flurry of laboratory experiments involving mostly disturbed, homogeneous soil systems.
Laboratory testing of convection-dispersion type solute transport models was especially popular
at that time. Field-scale testing of flow/transport theories and models did not really start until
the mid 1970s in part spurred by the development of better numerical models and theincreased availability of mainframe and personal computers. Process-based simulation models
were subsequently developed for the purpose of integrating available knowledge, and forformulating and evaluating alternative agricultural soil, crop, and water management practices.
While most of these simulations initially focused primarily on water quantity issues (irrigation,
drainage, soil erosion), they increasingly also involved agricultural water quality issues
associated with pesticides, fertilizers, soil salinity, and sewage sludge.
Soil physicists must remain concerned about the physical environment of plants. Agricultural
production unquestionably will remain a critical issue as the world’s population keeps increasing
in an era of limited soil and water supplies. At the same time, however, soil scientists/physicists
have also a responsibility to address soil and groundwater pollution problems caused byagricultural and non-agricultural sources at local, regional, and global scales. It is critical thatthese environmental problems be addressed in a systematic manner involving efforts which crossdisciplinary, institutional, and geographic barriers. The need for a broad multidisciplinary
perspective should be obvious in view of the many physical, chemical and biological processes
known to affect the behavior and transport of water and dissolved constituents in the subsurface
environment. The following examples illustrate the need for a more broadly defined frameworkof soil physics research.
First, soil physics research can no longer be confined to the upper soil horizons between the soil
surface and a depth of 2 m traditionally viewed as the lower boundary of the soil medium (e.g.,
the bottom of the root zone). Soil physicists must be concerned with flow and transport
processes in the entire vadose zone between the soil surface and the groundwater table, even
if the vadose zone is hundreds of meters deep as is often the case in arid and semi-arid areas.
We must be able to quantify how flow and transport processes in the vadose zone affect thequality of underlying or down-gradient groundwater reservoirs. In situations involving shallow
water tables, we must be concerned not only with downward flow and transport, but conversely
also with possible upward flow of water and dissolved constituents. Because of the dangers of(secondary) soil salinization, shallow water table management is a special concern in arid andsemi-arid regions. An integrated soil/aquifer research approach requires interaction withhydrologists and geologists.
Second, soil physicists are increasingly becoming participants in global-scale hydrologic research.
Especially needed (e.g., Wood, 1991) are improved models of land-surface hydrologic processes
in general circulation models for predictions at regional, continental, and global scales, as well
as appropriate field studies over a range of scales insofar as these scales impact hydrologicalprocesses at the larger scales. Soil physicists should be especially interested in the land surface
component of global hydrological models, including the development of realistic physicallybased models for biosphere-atmosphere interactions to estimate the transfer of energy, mass,and momentum between the vegetated surface of the earth and the atmosphere.- With theirtraditional focus on soil water dynamics and microclimatology, soil physicists should be natural
cooperators in research with hydrologists, climatologists and global climate modelers.
Third, soil physicists must be interested in media other than the weathered and fragmented
outer layer of the earth’s terrestrial surface usually referred to as soil (e.g., Hillel, 1980).
Development of realistic models for unsaturated flow and solute transport in fractured rock andother geologic materials is crucial for establishing suitable subsurface waste disposal sites. For
example, much research in several countries has been, or is being, carried out in connection
with the safe disposal of spent nuclear fuel and high-level radioactive waste in unsaturated
fractured rock (e.g., Evans and Nicholson, 1987). The same is true for the performance of clay
liners as barriers to the migration of contaminants from waste disposal sites. Soil physicists can,and should, provide much expertise to these projects.
Fourth, soil physicists must be concerned with fluids other than water, which traditionally has
been the main focus in agricultural research. For example, subsurface pollution by nonaqueous
phase liquids (NAPL’s) is increasingly becoming a problem in many countries. While
multiphase (air-water-oil) fluid systems pose additional problems, soil physicists have long
worked almost exclusively on the measurement and modeling of two-phase air-water systems;
they must provide sorely needed input into multiphase flow research. Soil physicists have made,and should continue to make, important contributions to theoretical studies of NAPL transport,

including the measurement and scaling of multiphase capillary-pressure saturation andpermeability relationships. The same is true for vapor phase transport involving volatile organic
chemicals and other constituents.
Fifth, soil physicists must be equally concerned with agricultural chemicals (fertilizers, pesticides,
soil salinity) as well as with nonagricultural chemicals, including radionuclides, toxic traceelements, and certain volatile organics. Colloid-facilitated contaminant transport should
similarly be a concern to soil physicists.
The above examples demonstrate that soil physicists are or should be involved in a broad arrayof research topics well beyond the traditional boundaries of soil physics. In the remainder ofthis paper I will address a few selected topics on flow and transport in the unsaturated zone,
notably on flow and transport modeling, subsurface heterogeneity, and parameter estimation.
I will also briefly discuss the need for more effective communication between soil physicists
globally, and stress the educational implications of a broader soil physics/soil science research
perspective.
-Water Flow and Solute Transport
Process-based descriptions of water and solute movement in the unsaturated zone are generally
based on the classical Richards equation for variably-saturated water flow and the convectiondispersion equation for solute transport. These equations, summarized in Table 1, have beensolved for a variety of one- and multi-dimensional applications. This paper focuses primarilyon conceptual issues, and hence no attempt is made to provide an exhaustive review of available
literature on vadose zone flow and transport modeling. Recent reviews of analytical andnumerical modeling are given by Giiven et al. (1990) and Sudicky and Huyakorn (1991).
Models based on Eqs. (1) and (2) in Table 1 have provided, and will continue to provide,
indispensable tools in research and management. For example, they may be used for designingor analyzing specific experiments on water and solute movement; for extrapolating information
gained from a limited number of field studies to different soil, crop, and climatic conditions, aswell as to different tillage and water management schemes; for evaluating the comparative
effects of alternative soil and water management practices and chemical application
technologies on crop production and groundwater pollution; for remediation purposes; and forrisk assessment studies involving specific pollution cases. The importance of model
development and use should not be underestimated: models provide a means for organizing andintegrating knowledge, give guidance to an experimental research program, while also being
helpful for many practical management applications.
While models based on Eqs. (1) and (2) are important tools, they are also subject to a large
number of simplifying assumptions which limit their applicability to many problems in the field.
It my be instructive to list some of these limitations and assumptions. For example, theequations assume that (a) the air phase plays a relatively minor role during unsaturated flow,
and hence that a single equation can be used to describe the variably-saturated flow process,(b) Darcy’s equation is valid at both very low and very high flow velocities (including thoseoccurring in structured soils), (c) the osmotic and electro-chemical gradients in the soil water
potential are negligible, (d) the fluid density is independent of the solute concentration, and (e)matrix and fluid compressibilities are relatively small. The equations are further complicatedTable 1. Governing equations for one-dimensional variably-saturated waterflow (Eq. 1) and convective-dispersive solute transport (Eq. 2).

C(h)2 = ;[K(h)g -K(h)] + .s-+eog

I + duw _ aatat-qc)++(2)

where C is the soil water capacity, i.e., the slope of the soil water retention curve, B(h); 0 is thevolumetric water content, h is the soil water pressure head (negative for unsaturated conditions), I
is time, z is distance from the soil surface downward, K is the hydraulic conductivity as a functionof h or 0, s is the solute concentration associated with the solid phase of the soil, c is the soluteconcentration of the fluid phase, p is the soil bulk density, D is the solute dispersion coefficient, S
and 4 are sources and sinks for water and solutes, and 4 is the volumetric fluid flux (Darcy’s law):--q = -K(h) ah + K(h)az

For conditions of steady-state water flow in homogeneous soils, neglecting the source/sink term S
and I$, and assuming linear sorption such that the adsorbed concentration (s) is linearly related to thesolution concentration (c) through a distribution coefficient, k (i.e., s=kc), Eq (2) reduces to thestandard convection-dispersion equation (CDE):
j& =Da’c _ ,,eata22 a~(4)

where v = q/B is the averagepore water velocity, and R = 1 + pkl0 is the solute retardation factor.

by (f) the hysteretic nature of the retention function, @(it), (g) the extreme nonlinearity of thehydraulic conductivity-function, K(k), (h) the lack of accurate and cheap methods for measuring
the unsaturated hydraulic properties, (i) the extreme heterogeneity of the subsurface
environment, and (j) inconsistencies between the scale at which the hydraulic and solutetransport parameters in Eqs. (1) and (2) are usually measured, and the scale at which thepredictive models are being applied. In addition, Eqs. (1) and (2) are formulated for isothermal
soil conditions. In reality, most physical, chemical and microbial processes in the soil are
strongly influenced by soil temperature. This also applies to water flow itself, including theindirect effects of temperature on the unsaturated soil hydraulic properties. Hence, a complete
description of unsaturated zone transfer processes requires also consideration of heat flow andits nonlinear effect on most processes taking place in the soil-plant system.
Other limitations in the use of Eqs. (1) and (2) involve the source/sink terms S and Q in theflow and transport equations. The term S accounts primarily for water uptake by plant roots.
While widely different approaches exist for simulating water uptake (e.g., Molz, 1981), mostmodels currently used for water uptake are essentially empirical and contain parameters that
depend on specific crop, soil, and environmental conditions. Much research is still needed toderive process-based descriptions of root growth and root water uptake as a function of water,
salinity, temperature, and other stresses in the root zone, and to couple these descriptions withsuitable crop growth models. Similarly, the source/sink term C$ in Eq. (2) accounts for nutrient
uptake and/or a variety of chemical and biological reactions and transformations insofar asthese processes are not already included in the sorption/exchange term aps/at. These processes
can be highly dynamic and nonlinear in time and space, especially for nitrogen and pesticide
products. For example, among the nitrogen transformation processes that may need to beconsidered are nitrification, denitrification, mineralization, and nitrogen uptake by plants(Stevenson, 1982). For microbially induced organic and inorganic transformations, thedegradation process should also consider the growth and maintenance of soil microbes.
Several simplifying assumptions are often also invoked when defining the sorption/exchange
term Lips/at in Eq. (2). For example, the use of a constant retardation factor R in Eq. (4)assumes that a linear equilibrium isotherm can be used to describe solute interactions between
the liquid and solid phases of the soil. While the use of a linear isotherm can greatly simplify
the mathematics of a transport problem, sorption and exchange reactions are generally
nonlinear and usually depend also on the presence of competing species in the soil solution.
This, in turn, may require the coupling of Eq. 2) with suitable solution chemistry and/or ion
exchange models (Mangold and Tsang, 1991; (simunek and Suarez, 1994). The equilibrium
assumption itself for sorption or exchange is often also questionable. A number of chemicalkinetic and diffusion-controlled “physical” models have been used to describe nonequilibriumtransport. Attempts to model nonequilibrium transport have usually involved relatively simple
first-order (one-site) kinetic rate equations. More refined nonequilibrium models have invoked
the assumptions of two-site or multi-site sorption, and/or two-region (dual-porosity) transport
involving solute exchange between mobile and immobile liquid transport (e.g., Wagenet, 1983;
Brusseau et al., 1992). Models of this type have resulted in-better descriptions of observedlaboratory and field transport data, in part by providing additional degrees of freedom for fitting
observed concentration distributions. It should be noted that similar problems related tononlinear and nonequilibrium‘sorption also pertain to the transport of organic solutes (MacKay
et al., 1985; Pinder and Abriola, 1986). Depending upon the type of chemical involved, models
to predict organic transport in the unsaturated zone may also need to account for volatilization,
microbial, chemical and photochemical transformations, and often multiphase flow withpartitioning between different fluid phases (McCarty et al., 1984; Widdowson et al., 1988;
Kaluarachchi and Parker, 1989).
The above limitations and assumptions regarding Eqs. (1) through (4) in Table 1 provide many
opportunities for research. Unfortunately, there is ample evidence to suggest that solutions ofthese classical models, no matter how refined to include the most relevant chemical andmicrobiological processes and soil properties, fail to accurately describe transfer processes inmost natural field soils. The one major factor responsible for this failure is the overwhelming
heterogeneity of the subsurface environment. Heterogeneity occurs at a hierarchy of spatial andtime scales (Wheatcraft and Cushman, 1991), ranging from microscopic scales involving timedependent chemical sorption and precipitation/dissolution reactions, to intermediate scales
involving the preferential movement of water and chemicals through macropores or fractures,and to much large scales involving the spatial variability of soils across the landscape. Several
lines of research are being followed to deal with the different types of soil heterogeneity. On
the one hand, subsurface heterogeneity can be addressed in terms of process-based descriptions
which attempt to consider the effects of heterogeneity at one or several scales (kinetic sorption,
preferential flow, field-scale spatial variability). On the other hand, subsurface heterogeneity
is often also addressed using stochastic approaches which incorporate certain assumptions about
the transport process in the heterogeneous system. Much can be learned from both approaches.

Preferential flow in structured media will be discussed below as an example of a process-based
approach, followed by a brief overview of stochastic approaches.
Preferential Flow in Structured Media

Most soils contain macropores or other structural features associated with interaggregate pores,
earthworm and gopher holes, decayed root channels, or drying cracks in fine-textured soils.
Water and dissolved chemicals can move into and through such structured media with widelydifferent velocities, thereby creating highly non-uniform flow fields. This process, often referred
to as preferential flow, macropore flow, or short-circuiting, has been shown to occur not onlyin aggregated field soils (Beven and Germann, 1982; Gish and Shirmohammadi, 1991) andunsaturated fractured rock (Wang, 1991), but also in seemingly homogeneous soils because of
fingering or some other unstable flow process (Hillel, 1993). While preferential flow cansignificantly affect infiltration rates at the soil surface, and unsaturated water flow in general,
its main implications are likely the accelerated movement of surface-applied fertilizers,
pesticides or other pollutants into and through the unsaturated zone.
Transport processes in structured media are often described using dual-porosity, two-region, or
bi-continuum models. Models of this type assume that the medium consists of two interacting
pore regions, one associated with the macropore or fracture network, and one with themicropores inside soil aggregates or rock matrix blocks. Different formulations arise depending
upon how water and solute movement in the micropore region are modeled, and how the
micropore and macropore regions are coupled.
Geometry-Based Models. A rigorous analysis of transport in structured soils can be made whenthe medium is assumed to contain geometrically well-defined cylindrical, rectangular or othertypes of macropores or fractures. Models may be formulated by assuming that the chemical istransported by convection, and possibly by diffusion and dispersion, through the macropores,
while diffusion-type equations are used to describe the transfer of solutes from the larger pores
into the micropores of the soil matrix. Table 2 illustrates such a model for transport througha soil system containing a series of parallel rectangular fractures. Similar models may beformulated for other aggregate geometries. Geometry-based transport models like those in
Table 2 have been successfully applied to laboratory-scale experiments as well as to selected
field studies involving mostly saturated conditions. As an example, Figure 1 shows calculated
and observed Cl effluent curves from a 76-cm long undisturbed column of fractured clayey till.
The extremely skewed (nonsigmoidal) shape of the effluent curve is a direct result of water anddissolved chemical moving mostly through the fractures and bypassing the soil matrix.

While geometry-based models are conceptually attractive, they may be too complicated forroutine applications since structured field soils usually contain a mixture of aggregates of various
sizes and shapes. More importantly, the problem of macropore and fluid flow continuity is noteasily addressed with geometry-based flow models. The main objective of preferential flow
studies should not be to quantify the macropore or fracture network per se, but rather the fluid
flow network contained in some fashion within the fracture network. The issue of fluid flow
continuity is critical, in part also since the macropore domain itself may be subject topreferential flow and channeling during unsaturated conditions. Finally, it is to be expected thatthe preferential flow paths will change with the degree of saturation during unsaturated flow.
Some of these features are more easily considered using first-order models as discussed below.

Table 2. Governing equations for transport through soil containing
parallel rectangular voids.

(5)

where. the subscriptsfand m refer to the interaggregate (fracture f) and intra-aggregate (matrix m)pore regions, respectively, c,(z,x,f) is the local concentration in the aggregate, x is the horizontal
coordinate, and Da is the effective soil or rock matrix diffusion coefficient. Equation (5) describes
vertical convective-dispersive transport through the fractures, while (6) accounts for linear diffusion
in slab of width 2a in the horizontal (x) direction. Equation (7) represents the average concentrationof the immobile soil matrix liquid phase Equations (5) and (6) are coupled using the assumptionof concentration continuity across the fracture/matrix interface:C~(Z,W) = C,(Z,f)(8)

The water contents t9, and 8, in (6) are given in terms of the bulk soil volume, i.e.,9, = w,e,l_Ym = (1 - “$9,(9)where w, is the volume of the fracture pore system relative to that of the total soil pore system.

TIME, t (days)
Fig. 1. Measured breakthrough curve for Cl transport through fractured clayey till (opencircles; data from Grisak et al., 1980). The solid line was obtained with the exact solution(e.g., van Genuchten, 1985) of the model in Table 2 using parameter values given by Grisak
et al. (1980). The dashed line was obtained by ignoring dispersion in the interaggregate
region (D,=O) and allowing the fracture spacing to go to infinity (U--X).

Table 3. Governing equations for transport in structured media assuming first-ordersolute exchange between mobile and immobile liquid regions:where a is a first-order solute mass transfer coefficient characterizing diffusional exchange of solutes
between the mobile and immobile liquid phases. Notice that (10a) is identical to (5) for therectangular geometry-based model. The mass transfer coefficient is of the general form

where /3 is a geometry-dependent shape factor, and a is the characteristic length of the aggregate(e.g., the radius of a spherical or solid cylindrical aggregate, or half the width of a rectangularaggregate). Equation (1 lb) holds for a hollow cylindrical macropore for which (,=6/a where a now
represents the radius of the macropore and b the outer radius of the cylindrical soil mantle
surrounding the macropore. The value of j3 ranges from 3 for rectangular slabs to 15 for sphericalaggregates [van Genuchten and Dalton, 1990].

First-Order Models. Rather than using geometry-based transport models, many of thepreferential flow features can also be accounted for using models which assume first-orderexchange of solutes by diffusion between the macropore (mobile) and micropore (immobile)liquid regions. Table 3 summarizes the governing equations for such an approach, Models of
this type have been popularly used in petroleum and chemical engineering studies; they lumpthe effects of aggregate size, geometry, and the diffusion coefficient in the quasi-empirical masstransfer coefficient, CX. Because of their simplicity, first-order mobile-immobile models can be
included in one- and multi-dimensional numerical models with only minimal changes in thecomputer codes. The applicability of first-order rate models to transport processes at the scaleof laboratory soil columns has been demonstrated in a large number of studies in the soil
physics literature (Nielsen et al., 1986; van Genuchten et al., 1990).
Extension to Variably-Saturated Flow. One limitation of mobile-immobile transport models
such as those described in Tables 2 and 3 is the assumption that water flow is limited to themacropore region. This assumption is inconsistent with experimental observations whichindicate that water in the soil matrix is generally also mobile. Different types of models havebeen proposed to simulate transport in variably-saturated structured media. As an example,
Table 4 summarizes the dual-porosity model developed by Gerke and van Genuchtcn (1993).
This model assumes that the Richards equation for transient water flow and the convectiondispersion equation for solute transport can be applied to each of the two pore systems.
Similarly as for the first-order mobile-immobile transport models, water and solute masstransfer between the two pore systems is described with first-order rate equations. The modelin Table 4 contains two water retention functions, one for the matrix and one for the fracture
pore system, but three hydraulic conductivities functions: K,(/I,) for the fracture network K,_(k)for the matrix, and K,(h) for the fracture/matrix interface. The I@,) function is determinedby the structure of the fracture pore system, i.e., the size, geometry, continuity and wallroughness of the fractures, and possibly the presence of fracture fillings. Similarly, K,(/I,) is
determined by the hydraulic properties of single matrix blocks, and the degree of hydrauliccontact between adjoining matrix blocks during unsaturated flow. Finally, K,(h) is the effective
hydraulic conductivity function to be used in Eq. (14) for describing the exchange of water
between the two pore systems.

Table 4. Governing equations for variably-saturated flow and solute transport
in a dual-porosity medium [Gerke and van Genuchten, 1993]
Flow equations for the fracture (subscript f) and matrix (subscript m) pore systems are, respectively,C,? = ;(K,!$ -K,) - $

(l2a)

%

(12b)

a

rw

%

Cm= = -#mz -Km) + 1 -w,
where r, describes the rate of exchange of water between the fracture_- and matrix regions:(13)r, = cx,(h,-hmj
in which (Y, is a first-order mass transfer coefficient for water:P
ay = -Kaya2wwhere 0 and a are as defined in Table 3, Ka is the hydraulic conductivity of the fracture/matrix
interface, and y, (=0.4) is a dimensionless scaling factor. Solute transport for the fractures andmatrix are given by, respectively:

where r, is the solute mass transfer term given byr, = a(~,-c_) +

I

rW+p

rw9, cm/e
,

rw20
rw<o

(16)

in which Q is the same as used in the first-order mobile-immobile model (Table 3). The first termon the right-hand side of (16) specifies the diffusion contribution to r,, while the second term gives
the convective contribution. The above variably-saturated dual-porosity transport model reduces tothe first-order model (Table 3) for conditions of steady-state flow in the fracture (macropore) region
and no flow in the matrix pore system (q_=r, =0).

Several important features of preferential flow are illustrated here by using the model in Table
4 to calculate the infiltration of water at a constant rate of 50 cm/day into a 40-cm deep
structured soil profile having an initially uniform pressure head of -1000 cm. Water is allowed
to infiltrate exclusively into the fracture pore system, thus assuming that the matrix pore systemat the soil surface is sealed. The hydraulic properties of the fracture and matrix pore systems(Fig. 2) are indicative of relatively coarse- and fine-textured soils, respectively. The simulations
assume a macroporosity of 5% (w~=O.OS), and rectangular aggregates (/3=3) having a width of
2 cm (a = 1 cm). The hydraulic parameters for &(/a) were assumed to be the same as those for
K&J, except for the saturated hydraulic conductivity which was decreased by a factor of 100.
Figure 3 shows simulated pressure head and watercontent distributions during infiltration. The
results indicate a rapid increase in the pressure head of the fracture pore system, but arelatively slow response of the matrix (Fig. 3a). The resulting pressure head gradient between
the two pore systems causes water to flow from the fracture into the matrix pore system (Fig.
3b), thus increasing the water content of the matrix (Fig. 3c). Significant pressure headdifferences between the two pore systems are still present when the infiltration front in thefracture system reaches the bottom of the soil profile after approximately 0.08 days (Fig. 3a).
Notice that the water transfer rate, I’,,,, is highest close to the infiltration front, and gradually
decreases toward the soil surface (Fig. 3b). The shapes of the r,-curves reflect the combined
effects on I’, of the pressure head difference between the two pore regions (which decreases
in time) and the value of Ku (which increases in time) at any point behind the wetting front.

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Volumetric Water Content, 0

-lo+

-lb

-10'

-1d

-10’

-1u

-

Pressure Head, /I (cm)Fig. 2. Water retention (a) and hydraulic conductivity (b) functions of a dual-porosity
medium involving the fracture network (l), the matrix pore system (2), the composite
medium (3), and the conductivity of the fracture/matrix interface (4) (after Gerke and van
Genuchten, 1993).

Figure 4 shows the simulated concentration distributions. Results are for the infiltration ofsolute-free water into a structured medium having a relative initial concentration of 1. The
calculations show that the solute concentrations in the fracture pore system initially decrease
rapidly as solute-free water infiltrates (Fig. 4a). Water with relatively low solute concentration
subsequently flows from the fracture into the matrix pore system (see also Figure 3b). At thesame time, however, solutes begin to diffuse back from the matrix into the fracture pore systembecause of the large concentration gradients which develop between the two pore systems (Fig.
4a). The net solute transfer rate, P,, eventually becomes negative, indicating a net transferfrom the matrix into the fracture pore system (Fig. 4b). The solute mass in the matrix poresystem (19,cJ initially decreases only slightly (t=0.0l days in Figure 4c), but starts to decrease
more rapidly at later times (tr0.04 days). The results in Figure 4 illustrate the extremelytransient and complicated nature of transport in a structured medium involving verticalconvective transport and dispersion; and horizontal mass transfer by convection and diffusion.
Simulations such as those shown in Figures 3 and 4 may be used to explain previously observed
effects of several parameters on solute leaching during transient flow, including soil surface
boundary condition (Bond and Wierenga, 1990) water application rate (White et al., 1986;
McLay et al., 1991), and initial condition (Kluitenberg and Horton, 1990).
The potential value of process-based preferential flow simulations is further illustrated in Figure
5 which shows the sensitivity of the infiltration process to changes in the hydraulic conductivity
K. of the fracture/matrix interface. Results obtained with a relatively large saturated
conductivity, K_ , of 1 cm/day (equal to the matrix conductivity) closely approximate thelimiting case of pressure head equilibrium (Fig. 5a) with little or no preferential flow. The
equilibrium moisture front reaches a depth of only 5 cm after 0.02 days. The water transfer
rates shown in Figure 5b are So high that the two pore systems quickly approach equilibrium
when K,,a = 1 cm/day. However, for the smallest Ks,a (0.001 cm/day) water percolated rapidly
downward through the fracture pore system to a depth of 35 cm during the same time period(t=0.02 days or 29 minutes). This last situation represents an extreme case of preferential flow
with significant pressure head differences between the two pore systems (Fig. 5a).
The results in Figure 5 are important since they indicate that equilibrium between the fracture
and matrix pore systems should be expected when the hydraulic conductivity, K,,*, of the matrix/
fracture interface is equal to the conductivity of the soil matrix (assuming a fracture spacing of
2 cm). For preferential flow to initiate, K,,a must be much less than KS,: of the matrix. Thisconclusion is consistent with experimental studies which suggest that a soi1 aggregate may havea higher local bulk density (and hence lower conductivity) near its surface than in the aggregate
center, presumably because of the deposition of organic matter, fine-texture mineral particles,
or various oxides and hydroxides on the aggregate exteriors or macropore walls. For example,
Wilding and Hallmark (1984) noted that ped argillans have been demonstrated to markedly
reduce rates of diffusion and mass flow from the ped surface to the soil matrix. Cutans,consisting of coatings with modified physical, chemical or biological properties, may havepreferred orientations parallel to the aggregate surfaces; they are common in many clayey
horizons as a consequence of both natural and tillage-induced deformation. Unsaturated
fractured rock formations may exhibit similar features, i.e., fracture skins (Moench, 1984) orother types of coatings (Wang and Narasimhan, 1985; Pruess and Wang, 1987; Thoma et al.,1992) made up of fine clay particles, calcite, zeolites, or silicates, which may reduce thehydraulic conductivity. Preferential flow within the macropores or fractures themselves also

F04rh.‘,““,““,““,‘,,
1~025Pressure Head, h (cm)

Water Transfer Rate, I’, (l/day)

0.30

055

0.40

0.45

,

Volumetric Water Content, 0

Fig. 5. Simulated distributions versus depth of (a) the pressure head, 12, (b) the watertransfer rate, r,,,, and (c) the total volumetric water content, 0, for different values of theinterface hydraulic conductivity, K_ (t= 0.02 days, a = 1 cm).
contributes to a lower effective K&L). Such a situation would restrict water and solute exchange

between the two pore systems (notably infiltration into the matrix) to only a small portion of
the total interface area (Hoogmoed and Bouma, 1980), even in capillary-size pores (Omoti andWild, 1979). Hydrophobic fracture surfaces can similarly restrict the exchange between the twopore systems.
Application of the model in Table 4 requires several hydraulic and other parameters which are
not easily measured. Estimates for the K,- and Km -functions (Fig. 2) may be obtained byassuming that K, is primarily the conductivity function in the wet range, while K,,, is the
conductivity in the dry range (Othmer et al., 1991; Durner, 1992; Peters and Klavettcr, I988).
Obtaining accurate estimates for the hydraulic properties of the fracture pore system from thecomposite curves requires that the hydraulic functions be very well defined in the wet range.
This problem is indirectly demonstrated by Figure 2a which was obtained by assuming that thefracture pore system comprises 5% of the porous medium. Notice that the retention functionof the matrix differs only minimally from that of the composite medium. Hence, it may be very
difficult in practice to estimate separate soil water retention curves of the fracture and matrix
pore systems using bulk soil measurements which generally contain some noise. By contrast,
it appears more promising to assess the contributions of macropores from an accurately
measured bulk hydraulic conductivity function near saturation (e.g., Smettem and Kirby. 1990).

Finally, we note that the dual-porosity model in Table 4 assumes applicability of the Richards
equation, and hence of Darcy’s law. This assumption may not be strictly correct for the fracture
pore system. However, given the uncertainties in all of the physical and chemical processes
related to preferential flow, the real issue is not necessarily the validity of Darcy’s law as such,
but whether Darcy’s law - even if formally invalid - can provide a useful description of thepreferential flow process. Alternative descriptions of the flow regime in fractures, such as
Manning’s equation for turbulent overland flow, kinematic wave theory, or simple gravity-flow
models, may be too elaborate for routine use. Moreover, some of these approaches do not haveprovisions for flow to occur from the micropores back into the fractures, e.g., at or near thebottom boundary of a coarse-textured soil horizon overlaying a fine-textured horizon.
Stochastic Approaches
Dual-porosity type preferential flow models represent attempts to deal with pore structure
heterogeneity at spatial scales somewhere intermediate between laboratory-scale measurements
and the larger field scale. As such they are useful for predicting predominantly verticaltransport in structured but areally homogeneous field soils. Unfortunately, few soils are areallyhomogeneous. This is illustrated in Figure 6 which shows measured bromide concentrations
some 400 days after application of an areally uniform instantaneous bromide pulse to thesurface of a field soil in Switzerland. The unequal distribution of the tracer, especiallyhorizontally along the transect, raises important questions about how to design effective fieldsampling programs, which field instrumental methods to use for sampling, and how to simulate
the heterogeneous field-scale transport process.
The enormous variability of the subsurface environment, and the imprecision with whichparameters and processes can be measured, has led to the adoption of stochastic models andgeostatistical procedures to assist in the prediction and monitoring of contaminant transport in

2.5

5.0

7.5

IO.0

12.5

15.0

Distance Along Transect, x (m)Fig. 6. Observed gravimetric bromide concentrations (mg per kg dry soil) 399 days after
application of a bromide tracer pulse to the soil surface (after Schulin et al., 1987).

the unsaturated zone. A large number of stochastic approaches are currently available (Dagan,
1989; Jury and Roth, 1990), including scaling theories, Monte Carlo methods, and stochasticcontinuum models. A common assumption of most stochastic transport models is thatparameters are treated as random variables with values assigned according to a given
probability distribution. In practice, the stochastic approach is often used with several
simplifying assumptions, including (1) the stationarity hypothesis which assumes that a randomparameter has the same probability density function (pdf) at every point in the field, and (2)the ergodicity hypothesis which states that ensemble averages can be replaced by spatial
averages, and that spatial replicates can be used to construct the appropriate pdfs for thetransport parameters.
Current scaling theories applied to field-scale flow and transport problems have evolved fromthe early work of Miller and Miller (1956) on microscopic geometric similitude. The approachconsiders different regions of a heterogeneous field soil to be similar if their microscopic
geometric structures are scale magnifications of each other. Transport parameters at any point
within a given field soil are related to the parameters at an arbitrary reference point throughlength scale ratios, or scaling factors. Hydraulic conductivity and soil water retention
parameters at a particular location in the field are then calculated from those of the reference
soil by means of prescribed equations (Nielsen et al., 1983; Sposito and Jury, 1985). Recent
work (Jury et al., 1987; Vogel et al., 1991) suggests that two or more scaling factors may beneeded for soils which are not strictly similar. Most previous applications of scaling theory tofield problems assumed that the scaling factor is a random variable characterized by a certain
probability density function. The method has been a central part of the stochastic flow andtransport models of Bresler and Dagan (1983) and Dagan and Bresler (1983).
Monte Carlo simulations assume that the flow and transport parameters are random variables
with values assigned from a joint pdf. The water flow or solute transport equations arerepeatedly solved using coefficient values from the assumed pdf until a large number of possible
outcomes has been generated. These outcomes are then used to calculate sample means andvariances of the underlying stochastic transfer process. The method may be used todemonstrate, among other things, that macrodispersion in field soils is affected primarily bypore-water velocity variations, and much less by local dispersion phenomena (AmoozegardFard et al. 1982; Persaud et al., 1985).
Stochastic continuum models were initially used primarily in groundwater studies (e.g., Gelhar
et al., 1979; Gelhar and Axness, 1983). In these models all random variables are represented
by the sum of their mean values plus random fluctuations which, when substituted into theconvection-dispersion equation (Eq. 4) lead to a new mean transport model with additionalterms. The modified model may be evaluated by deriving first-order approximations for thefluctuations and solved by means of Fourier transforms. The approach leads to, among otherthings, a macro-scale dispersion coefficient whose value is reached asymptotically as distance
and/or time increase. Spatial correlations of solute velocity variations characterized by its
autocorrelation function, have been shown to play important roles in the derivation of theasymptotic convection-dispersion equation (Sposito et al., 1986).
A different continuum approach was followed by Simmons (1982) who neglected the dispersion
coefficient D in Eq. (4), and developed a formal theoretical approach using the pore-watervelocity and the travel time as random variables. Jury (1982) initially also neglected D in his
development of the transfer function model (TFM) of solute transport. The TFM involves theuse of a probability density function, f&), of travel times from the soil surface down to somereference depth L. The travel time probability density function for many transport experiments
is given by a lognormal distribution. The flux concentration in the profile is represented witha convolution integral off&) and the imposed flux concentration at the soil surface. Transfer
function models are expected to find increasingly wider applications in subsurface solutetransport as the underlying theory is being strengthened by the incorporation of a variety ofphysical, chemical and biological processes (Jury and Roth, 1990; Sardin et al., 1991).
Several other statistical approaches exist, including the use of continuous Markov processes(e.g., Knighton and Wagenet, 1987), fractal-mathematics (Wheatcraft and Cushman, 1991),
random walk methods (Kinzelbach, 1988), and a variety of procedures based on moment
analysis (Cvetkovic, 1991). More work in these areas of research can be expected in the near
future. While useful, the real challenge is to eventually include in the stochastic approach thosephysical, chemical and microbiological processes and soil properties known to affectmacroscopic transport behavior. In other words, the challenge is to combine stochastic behavior
with process-based knowledge that so painstakingly has been obtained over the past several
hundred or more years. It would be a major error to ignore this knowledge and to resort toa purely empirical-statistical description of field-scale transport processes. Moreover, several
farming practices or edaphic factors sometimes combine to produce an apparent stochastic
situation which, upon closer inspection, may be described or interpreted more advantageously
in terms of deterministic processes or soil properties. An excellent example is given by Rhoades(1994) on field-scale variations in soil salinity measured with standard methods as well as using
an automated, tractor-driven measurement system.
Parameter Estimation
The past several decades has seen a tremendous increase in the development and use ofcomputer models simulating various aspects of water and solute movement in the unsaturated
zone. The rapid proliferation of computer models in research and management will likely
continue as computer costs keep decreasing and the need for more realistic predictions
increases. Unfortunately, the simulation of field-scale processes requires considerable effort inquantifying a large number of spatially and temporally variable soil hydraulic and solutetransport parameters. It is likely that the completeness of experimental data, and the accuracyof the estimated model parameters, eventually will become critical factors determining theusefulness of site-specific simulations. Many of our current methods for measuring relevant
unsaturated flow and transport parameters are largely those that were introduced several
decades ago (Klute, 1986; Dane and Molz, 1991). New methods and technologies ofmeasurement are critically needed to keep pace with our ability to simulate increasingly
complex laboratory and field systems. New or improved instrumental methods must be
developed to more accurately and cost-effectively measure pertinent physical and chemical soil
properties, including especially the unsaturated hydraulic functions. At the same time, given
suitable instrumentation, challenges are to improve current methods for sampling design andoptimal analysis of laboratory and field experiments.
Instrumental Techniques. A number of potentially very powerful geophysical and related
approaches are now available for rapid non-invasive as well as in-situ measurements within theunsaturated zone. Among the non-invasive methods are (a) electromagnetic methods (EM) for
characterizing soil salinity and water content, (b) ground penetrating radar (GPR) methods for
detecting subsoil layers (e.g., those causing preferential flow), (c) nuclear magnetic resonance(NMR) techniques for studying soil porosity and permeability, and (d) X-ray computed
tomography (CT) for studying, among other things, the three-dimensional geometry andcontinuity of macropores. Similarly, time-domain reflectometry (TDR) methods are increasingly
being used for in-situ measurements of water content and/or salinity distributions, both in thelaboratory as well as in the field. The same is true for remote sensing techniques to obtainspatial and temporal distributions of soil moisture. More sophisticated equipment formeasuring the soil hydraulic conductivity is also becoming available, including a variety of new
or improved borehole and disk permeameters.
Inverse Methods. Inverse procedures provide potentially powerful methodologies to moreaccurately estimate a large number of water flow, solute transport, and other parameters. Early
studies using inverse procedures focused primarily on solute transport problems by coupling
appropriate parameter optimization methods with analytical solutions of the convectiondispersion equations (e.g., van Genuchten, 1981). Starting with the studies of Zachmann et al.
(1981) and Dane and Hruska (1983), the method was later extended to the analysis of
unsaturated soil hydraulic functions. Computer models applicable to one-step and multi-steplaboratory outflow measurements are given by Kool et al. (1985), van Dam et al., (1990) and
Eching and Hopmans (1993). While initially applied primarily-to laboratory type experiments,
inverse methods are equally well applicable to field data, or some appropriate combination offield and laboratory data. As such, inverse methods are extremely flexible and could involve
nearly any combination of observed water flow data (water contents, pressure heads, fluxes),
solute transport data (concentrations, fluxes), and heat transfer data (temperature, heat fluxes),
to yield simultaneously estimated hydraulic parameters (hysteretic water retention andunsaturated conductivity parameters), solute transport parameters (pore-water velocities,
retardation factors, dispersion coefficients, degradation coefficient, nonequilibrium rate
coefficients), and/or heat transfer parameters (soil thermal conductivity, volumetric heat
capacity). Formally, there is also no reason to limit the analysis to only one- or two-dimensions.

Application of inverse methods to observed field data appears especially promising, e.g.,
procedures that estimate soil hydraulic parameters from infiltration experiments (Russo et al.,1991) rather than using the much more time-consuming gravity-drainage (instantaneous profile)approach (Shouse et al., 1992). Another advantage of inverse procedures, if formulated within
the context of a parameter optimization problem, is that a detailed error analysis of theestimated parameter can be incorporated more easily into the procedure. lnverse procedures
may prove to be very appropriate for estimating regional-scale effective soil hydraulicparameters, either by appropriately manipulating in-situ measurement of the hydraulicproperties as shown by Feddes et al. (1993), or perhaps by using remotely-sensed measurementsof the soil surface water content.
Concluding Remarks

This paper gives a brief overview of some of the issues facing soil physicists and others inaccurately predicting, measuring, and managing water and chemical transport processes in theunsaturated zone. The rapidly increasing problems of soil and water pollution require an
integrated, multi-disciplinary research effort involving contributions within the soil sciences from

soil physicists, soil biologists, soil morphologists and soil geneticists, as well as expertise fromoutside soil science in hydrology, geology, micrometeorology, the atmospheric sciences, andespecially mathematics and computer science. The implementation of a multidisciplinary
approach is predicated first and foremost by a level of awareness that such an approach isneeded in our own individual research effort, rather than necessarily being a member of a large
mission-oriented multidisciplinary research group. As such, I do believe that the internationalsoil science community needs to restructure some of its activities. For example, most membersof the Soil Science Society of America (SSSA) and the International Soil Science Society (ISSS)
operate and publish their research products within strongly disciplinary divisions orcommissions, with a concomitant lack of interaction across disciplinary boundaries.
We also must do a better job in documenting our research products. A large number of models
have been developed for predicting water and contaminant transport in the subsurface, as well
as for geostatistically analyzing spatially variable soil properties. These models must be
documented and made available more effectively to our peers, as well as to users in extension
or the private sector. An international clearinghouse (van der Heijde, 1993) for computermodels has been established for this purpose. A similar clearinghouse is sorely needed for
databases on flow and transport experiments carried out by soil physicists and others. Designingand conducting field experiments on water flow and solute transport, especially in naturally
heterogeneous field soils, can be extremely costly and time-consuming in terms of manpower,
equipment, and analytical expenses. It is imperative that such experiments be properly
documented and its results be made available to the scientific community at large. The sameis true for the documentation of measured unsaturated soil hydraulic-data sets (van Genuchten
and Leij, 1992). Some of these data sets are the results of thousands of hours of work. The
International Soil Science Society could play an important leadership role in such efforts to
document research results.
Finally, the need for a more broadly defined soil physics/soil science framework has manyeducational implications. We, as soil physicists, must become increasingly familiar with theprinciples of computer science, mathematics, fluid dynamics, hydrogeology, andmicrometeorology, and perhaps also with such topics as radiation physics, ecosystems
management, and global-scale hydrology. This training should not only occur at the level ofgraduate education, but also at the undergraduate level or even earlier (see also discussion in
Sposito and Reginato, 1992).
Literature Cited
1.

2.
3.
4.
5.

Amoozegard-Fard, A., D. R. Nielsen, and A. W. Warrick. 1982. Soil solute concentration
distributions for spatially varying pore water velocities and apparent diffusion coefficients.
Soil Sci. Soc. Am. J. 46(1):3-9.
Beven, K. J., and P. Germann. 1982. Macropores and water flow in soils. Water Resour.
Res. 18:1311-1325.
Bond, W. J., and P. J. Wierenga. 1990. Immobile water during solute transport inunsaturated sand columns. Water Resour. Res. 26:2475-2481.
Bresler, E., and G. Dagan. 1983. Unsaturated flow in spatially variable fields. 3. Solutetransport models and their application to two fields. Water Resour. Res. 19(2):429-435.
Brusseau, M. L., R. E. Jessup, and P. S. C. Rao. 1992. Modeling solute transport

6.
7.
8.
9.
1 0 .
11.

12.

13.
14.

15.

16.
17.
18.
19.
20.
21.
22.
23.
24.

influenced by multiprocess nonequilibrium and transformation reactions. Water Resour.
Res. 28(1):175-182.
Cvetkovic, V. D. 1991. Mass arrival of reactive solute in single fractures. Water Resour.
Res. 27(2):177-183.
Dagan, G. 1989. Flow and Transport in Porous Formations. Springer Verlag, Berlin.
Dagan, G., and E. Bresler. 1983. Unsaturated flow in spatially variable fields. 1.
Derivation of models of infiltration and redistribution. Water Resour. Res. 19(2):413-420.
Dane, J. H., and S. Hruska. 1983. In-situ determination of soil hydraulic properties
during drainage. Soil Sci. Soc. Am. J. 47:619-624.
Dane, J. H., and F. J. Molz. 1991. Physical measurements in subsurface hydrology. Rev.
of Geophys., Vol. 29, Supplement, pp. 268-279.
Dumer, W. 1992. Predicting the unsaturated hydraulic conductivity using multi-porosity
water retention functions. In: M. Th. van Genuchten, F. J. Leij and L. J. Lund (eds.),
Proc. Int. Workshop on Indirect Methods for Estimating the Hydraulic Properties of
Unsaturated Soils, pp. 185-202, Univ. of California, Riverside, CA.
Eching, S. O., and J. W. Hopmans. 1993. Inverse Solution of Unsaturated Soil Hydraulic
Functions from Transient Outflow and Soil Water Pressure Data. Land, Air and Water
Resources Paper No. 100021, Dept. of Land, Air and Water Resources, Univ. of
California, Davis.
Evans, D. D., and T. J. Nicholson (eds.). 1987. Flow and Transport Through Unsaturated
Fractured Rock. Am. Geophys. Union, Washington, D.C.
Feddes, R. A, G. H. de Rooij, J. C. van Dam, P. Kabat, P. Droogers, and J. N. M.
Stricker. 1993. Estimation of regional effective soil hydraulic parameters by inverse
modeling. In: D. Russo and G. Dagan (eds.), Water Flow and Solute Transport in Soils,
pp. 211-230, Springer Verlag, New York, NY.
Gelhar, L. W., and C. L. Axness. 1983. Three-dimensional stochastic analysis ofmacrodispersion in aquifers. Water Resour. Res. 19(1):161-180.
Gelhar, L. W., A. L. Gutjahr, and R. N. Naff. 1979. Stochastic analysis of macrodispersion in aquifers. Water Resour. Res. 15(6):1387-1397.
Gerke H. H., and M. Th. van Genuchten. 1993. A dual-porosity model for simulating
the preferential movement of water and solutes in structured porous media. Water
Resour. Res. 29(2):305-319.
Gish, T. J., and A. Shirmohammadi (eds.). 1991. Preferential Flow, Proc. of the National
Symp. December 16-17, 1991, Chicago, Illinois. Am. Soc. Agric. Eng., St. Joseph, MI.
Grisak, G. E., J. F. Pickens, and J. A. Cherry. 1980. Solute transport through fractured
media, 2. Column study of fractured till. Water Resour. Res. 16:731-739.
Giiven, O., F. J. Molz, and J. G. Melville. 1990. Modeling migration of contaminants inthe subsurface environment. In: J. Saxena (ed.), Hazard Assessment of Chemicals, pp. 203282, Hemisphere Publ. Co., New York, NY.
Hillel, D. 1990. Fundamentals of Soil Physics. Academic Press, New York, NY.
Hillel, D. 1993. Unstable flow: A potentially significant mechanism of water and solutetransport to groundwater. In: D. Russo and G. Dagan (cds.), Water Flow and Solute
Transport in Soils, pp. 123-135, Springer Verlag, New York, NY.
Hoogmoed, W. B., and J. Bouma. 1980. A simulation model for predicting infiltration
into cracked clay soil. Soil Sci. Soc. Am. J. 44:458-461.
Jury, W. A. 1982. Simulation of solute transport using a transfer function model, Water
Resour. Res. 18(2):363-368.

25.

Jury, W. A., and K. Roth. 1990. Transfer Functions and Solute Movement in Soil; Theory
and Applications. Birkhauser Verlag, Basel.
26. Jury, W. A., D. Russo, and G. Sposito. 1987. Spatial variability of water and solutetransport in unsaturated soil, II. Scaling models of water transport, Hilgardia. 55:33-55.
27. Kaluarachchi, J. J., and J. C. Parker. 1989. An efficient finite element method formodeling multiphase flow. Water Resour. Res. 25(1):43-54.
28. Kinzelbach, W. 1988. The random walk method in pollutant transport simulation. In:
E. Custodio, A. Gurgui, and J. P. Lobo Ferreira (eds.), Groundwater Flow and Quality
Modeling. pp. 227-246, D. Reidel, Hingham, MA.
29. Kluitenberg, G. J., and R. Horton. 1990. Effect of solute application method on
preferential transport of solutes in soil. Geoderma, 46:283-297.
30. Klute, A. (ed.). 1986. Methods of Soil Analysis, part 1. Physical and Mineralogical
Methods. Agronomy 9, 2nd ed., Am. Soc. Agron., Madison, WI.
31. Knighton, R. E., and R. J. Wagenet. 1987. Simulation of solute transport using acontinuous time Markov process, I. Theory and steady-state application. Water Resour.
Res. 23(10):1911-1916.
32. Kool, J. B., J. C. Parker, and M. Th. van Genuchten. 1985. ONESTEP: A nonlinear
parameter estimation program for evaluating soil hydraulic properties from one-stepoutflow experiments. Virginia Agric. Exp. Stat. Bull. 85-3, Blacksburg, VA.
33. MacKay, D. M., P. V. Roberts, and J. A. Cherry. 1985. Transport of organic
contaminants in groundwater. Environ. Sci. Tech. 19:384-392.
34. Mangold, D. C., and C.-F. Tsang. 1991. A summary of subsurface hydrological andhydrochemical models. Rev. of Geophys. 29(1):51-79.
35. McCarty, P. L., B. E. Rittman, and E. J. Bouwer. 1984. Microbial processes affecting
chemical transformations in groundwater. In: G. Bitton and C. P. Cerba (eds.),
Groundwater Pollution Microbiology, pp. 90-115, John Wiley, New York, NY.
36. McLay, C. D. A., K. C. Cameron, and R. G. McLaren. 1991. Effect of time of applicationand continuity of rainfall on leaching of surface-applied nutrients., Aust. J. Soil Res.
29( l):l-9.
37. Miller, E. E., and R. D. Miller. 1956. Physical theory of capillary flow phenomena.
J. Appl. Phys. 27:324-332.
38. Moench, A F. 1984. Double-porosity models for a fissured groundwater reservoir withfracture skin. Water Resour. Res., 20:831-846.
39. Molz, F. J. 1981. Models of water transport in the soil-plant system. Water Resour. Res.
17(5):1245-1260.
40. Nielsen, D. R., P. M. Tillotson, and S. R. Vieira. 1983. Analyzing field-measured soilwater properties. Agric. Water Manage. 6(1):93-109.
41. Nielsen, D. R., M. Th. van Genuchten, and J. W. Biggar. 1986. Water flow and solutetransport processes in the unsaturated zone. Water Resour. Res. 22(9):89S-108S.
42. Omoti, U., and A. Wild. 1979. Use of fluorescent dyes to mark the pathways of solute
movement through soils under leaching conditions, 1. Laboratory experiments. Soil Sci.
128(4):28-37.
43. Othmer, H., B. Diekkriiger, and M. Kutilek. 1991. Bimodal porosity and unsaturated
hydraulic conductivity. Soil Sci. 152(3):139-150.
44. Persaud, N., J. V. Giraldez, and A. C. Chang. 1985. Monte-Carlo simulation of noninteracting solute transport in a spatially heterogeneous soil. Soil Sci. Soc. Am. J. 49(3):562-568.

45.

Peters, R. R., and E. A. Klavetter. 1988. A continuum model for water movement in an
unsaturated fractured rock mass. Water Resour. Res. 24:416-430.
46. Pinder, G. F., and L. M. Abriola. 1986. On the simulation of nonaqueous phase organic
compounds in the subsurface. Water Resour. Res. 22(9):109S-119s.
47. Pruess, K., and J. S. Y. Wang. 1987. Numerical modeling of isothermal and nonisothermal flow in unsaturated-saturated fractured rock - A review. In: D. D. Evans, andT. J. Nicholson (eds.), Flow and Transport Through Unsaturated Fractured Rock. pp. 1 l22, Am. Geophys. Union, Washington, D.C.
48. Rhoades, J. M. 1994. Salinity assessment: Recent Advances and Findings. Trans. 15th
Int. Congress of Soil Science, Int. Soil Sci. Soc., Acapulco, Mexico.
49. Richards, L. A. 1931. Capillary conduction of liquids in porous mediums. Physics 1:318333.
50. Russo D., E. Bresler, U. Shani, and J. C. Parker. 1991. Analysis of infiltration events inrelation to determining soil hydraulic properties by inverse problem methodology. Water
Resour. Res. 27(6):1361-1373.
51. Sardin, M., D. Schweich, F. J. Leij, and M. Th. van Genuchten. 1991. Modeling thenonequilibrium transport of linearly interacting solutes in porous media: A review. Water
Resour. Res. 27(9):2287-2307.
52. Schulin, R., M. Th. van Genuchten, H. Fluhler, and P. Ferlin. 1987. An experimental
study of solute transport in a stony field soil. Water Resour. Res. 23(9):1785-1791.
53. Shouse, P. J., J. B. Sisson, T. R. Ellsworth, and J. A. Jobes. 1992. Estimating in situ
unsaturated hydraulic properties of vertically heterogeneous soils. Soil Sci. Soc. Am. J.
56(6):1673-1679.
54. Simmons, C. S. 1982. A stochastic-convective ensemble -method for representing
dispersive transport in groundwater. Report CS-2258, Electric Power Research Institute,
Palo Alto, CA.
55. Simdnek, J., and D. Suarez. 1994. Modeling unsaturated water flow, and CC& and solutetransport with major ion chemistry. Trans. 15th Int. Congress of Soil Science, Int. Soil Sci.
Soc., Acapulco, Mexico.
56. Smettem, K. R. J., and C. Kirby. 1990. Measuring the hydraulic properties of stable
aggregated soil. J. Hydrol. 117:1-13.
57. Sposito, G., and W. A. Jury. 1985. Inspectional analysis in the theory of water flow
through soil. Soil Sci. Soc. Am. J. 49(4):791-797.
58. Sposito, G., and R. J. Reginato. 1992. Opportunities in Basic Soil Science Research. Soil
Sci. Soc. Am., Inc., Madison, WI.
59. Sposito, G., W. A. Jury, and V. K. Gupta. 1986. Fundamental problems in the stochastic
convection-dispersion model of solute transport in aquifers and field soils. Water Resour.
Res. 22( 1):77-88.
60. Stevenson, F. J. (ed.). 1982. Nitrogen in Agricultural Soils. Agronomy No. 22, Am. Soc.
Agron., Madison, WI.
61. Sudicky, E. A., and P. S. Huyakorn. 1991. Contaminant migration in imperfectly known
heterogeneous groundwater systems. Rev. of Gcophys., Supplement 29:238-25 1.
62. Thoma, S. G., D. P. Gallegos, and D. M. Smith. 1992. Impact of fracture coatings on
fracture/matrix flow interactions in unsaturated, porous media. Water Resour. Res. 28:
1357-1367.
63. van Dam, J. C., J. N. M. Stricker, and P. Droogcrs. 1990. From One-Step to Multi-Step
Determination of Soil Hydraulic Functions by Outflow Measurements. Dept. of Water

Resources, Rep. No. 7, Agricultural University, Wageningen, The Netherlands.
64. van der Heijde, P. K. M. 1993. Identification and Compilation of Unsaturated/Vadose
Zone Models. Robert S. Kerr Environ. Res. Laboratory, Office of Research andDevelopment, U. S. Environ. Protection Agency, ADA, OK.
65. van Genuchten, M. Th. 1981. Non-Equilibrium Transport Parameters from Miscible
Displacement Experiments. Research Report No. 119, U. S. Salinity Laboratory, USDA,
ARS, Riverside, CA.
66. van Genuchten, M. Th. 1985. A general approach for modeling solute transport instructured soils. Proc. Int. Congress, Hydrageology of Rocks of Low Permeability, Jan.
7-11, Tucson, Arizona. Memoires Int. Assoc. Hydrogeol. 17(2):513-526.
67. van Genuchten, M. Th., and F. N. Dalton. 1986. Models for simulating salt movementin aggregated field soils. Geoderma 38:165-183.
68. van Genuchten, M. Th., and F. J. Leij. 1992. On estimating the hydraulic properties of
unsaturated soils. In: M. Th. van Genuchten, F.J. Leij, and L. J. Lund (eds.), Proc. Int.
Workshop on Indirect Methods for Estimating the Hydraulic Properties of Unsaturated
Soils, pp. 1-14, University of California, Riverside.
69. van Genuchten, M. Th., D. E. Rolston, and P. F. Germann (eds.). 1990. Transport of
Water and Solutes in Macropores. Special Issue, Geoderma 46(1-3):1-297.
70. Vogel, T., M. Cislerova, and J. W. Hopmans. 1991. Porous media with linearly variable
hydraulic properties. Water Resour. Res. 27(10):2735-2741.
71. Wagenet, R. J. 1983. Principles of salt movement in soil._ In: D. W. Nelson, D. E. Elrick,and K. K. Tanji (eds.), Chemical Mobility and Reactivity in Soil Systems, pp. 123-140,
Publ. No. 11, Soil Sci. Soc. of Am., Madison, WI.
72. Wang, J. S. Y. 1991. Flow and transport in fractured rocks. Rev. of Geophys., Vol. 29,
Supplement, pp. 254-262.
73. Wang, J. S. Y., and T. N. Narasimhan. 1985. Hydrologic mechanisms governing fluid flow
in a partially saturated, fractured, porous medium. Water Resour. Res. 21:1861-1874.
74. Wheatcraft, S. W., and J. H. Cushman. 1991. Hierarchical approaches to transport inheterogeneous porous media. Rev. of Geophys., Vol. 29, Supplement, pp. 261-267.
75. White, R. E., J. S. Dyson, Z. Gerstl, and B. Yaron. 1986. Leaching of herbicides throughundisturbed cores of a structured clay soil. Soil Sci. Soc. Am. J. 50:277-283.
76. Widdowson, M. A., F. J. Molz, and L. D. Benefield. 1988. A numerical transport modelfor oxygen- and nitrate-based respiration linked to substrate and nutrient availability inporous media. Water Resour. Res. 24(9):1553-1565.
77. Wilding, L. P., and C. T. Hallmark. 1984. Development of structural and microfabricproperties in shrinking and swelling clays. In: J. Bouma and P. A. C. Raats (eds.), Proc.
of the ISSS Symposium on Water and Solute Movement in Heavy Clay Soils, ILRI Publ.
37, pp. 1-18, Int. Inst. for Land Reclamation and Improvement, Wageningcn, The
Netherlands.
78. Wood, E. F. 1991. Global scale hydrology: Advances in land surface modeling. Rev. ofGeophys., Vol. 29, Supplement, pp. 193-201.
79. Zachmann, D. W., P. C. Duchatcau, and A. Klute. Simultaneous approximation of watercapacity and soil hydraulic conductivity by parameter identification. Soil Sci. 134: 157-163.