Chapter 2: Basic Relationships and the Soil Phases Under Equilibrium Conditions

2.1      Where we are going

2.1.1    Soil as a Three-Phase System

2.1.2    Volume and Mass Relationships

2.1.3    Liquid Content

2.2       Characterization of the Solid Phase

2.2.1    Texture and systems of classification

2.2.2    Sieving

2.2.3    Sedimentation

2.2.4    Surface Area of Soil Particles

2.2.5    Soil Structure and Bulk Density

Aggregation, bio-pores, packing geometries

2.3       Characterization of the Liquid Phase

2.3.1    Soil Water Content Measurement Methods

Gravimetric, dielectric, neutron probe, cosmic ray backscatter, remote sensing via radar and gravitation (SMAP, SMOS, GRACE); direct, indirect, and surrogate methods, errors.

2.3.2    Applications of Soil Water Content Information – The Water Balance across scales of soil profile, field, catchment, continent, and globe and seconds to decades.

2.3.3    Field Capacity, Wilting Point and Plant-Available Soil Water

2.4       The Energy State of Soil Water

2.4.1    Total Soil Water Potential

2.4.2    Gravitational Potential

2.4.3    Matric Potential

2.4.3.1 Water Properties Affecting the Matric Potential

2.4.4    Pressure Potential

2.4.5    Solute Potential

2.4.6    Water Potentials Under Equilibrium Conditions (Hanks, 1992 – pp. 31-37)

2.4.6.1.1          Constructing Potential Diagrams under Equilibrium Conditions

2.4.6.1.2          Water Potentials Under Nonequilibrium Conditions

2.5       Bringing gas, solids and liquids together

2.5.1    Capillarity and Surface Tension

2.5.1.1 What is surface tension

2.5.1.2 Deriving the LaPlace Equation

2.5.1.3 Temperature dependence of surface tension

2.5.1.4 Surface absorption (Tuller and Or Equation) and the myth of dry soil

2.5.1.5 Contact angle

Idealized, impacts of surface texture and defects, hysteresis

2.5.2    Soil Water Content-Energy Relationships

2.5.2.1 Soil Water Characteristic  (SWC)

2.5.2.2 Measurement of SWC Curves

2.5.2.3 Fitting Parametric SWC Expressions to Measured Data

2.5.2.4 Hysteresis in the Soil Water Characteristic Relation

2.5.2.4.1          The independent domain model for hysteresis

2.5.2.4.2          Analytical expressions for hysteresis

2.5.3    Measurement of Soil Water Potential Components

2.5.3.1 Tensiometer for Soil Matric Potential Measurement

2.5.3.2 Heat Dissipation in Rigid Porous Matrix for Measuring Matric Potential

2.5.3.3 The Psychrometric equation and Psychrometric methods for Measuring Water Potential

2.5.3.4 Piezometer for Measuring Hydrostatic Pressure Potential

2.5.3.5 Measurement of Soil Water Osmotic Potential

2.6       SUMMARY

2.7       Review Questions

2.8       REFERENCES

2.9       PROBLEMS

2.10     Appendices:

2.10.1  The thermodynamics of liquid and solid (Hamaker) interfaces

2.10.2  Advanced methods in soil moisture measurement (GPR, ERT, EMI, NMR, COSMOS

2.10.3  The independent domain model for hysteresis

2.1 Where we are going

We are concerned with describing the retention and movement of fluids in natural porous media. This objective fits neatly into a single sentence, however the challenge of attaining this goal remains to be met to the degree required to solve many contemporary problems in environmental sciences and engineering. Nonetheless, there are many powerful tools at our disposal. In this chapter we will cover the basic physical properties upon which much of the remainder of the text rests. If all goes well you will get a “feel” for both the physical processes which control flow through porous media and the conceptual and mathematical approaches we employ to characterize the processes. Having a well balanced understanding from both perspectives is essential; you must be able to guide your application of mathematical solutions based on a clear and accurate intuitive concept of what you are trying to describe.How do you describe a chunk of moist porous medium? What are some basic features which would allow you to abstract that chunk of material in terms of a set of mathematical parameters which will assist in describing: the way in which it might dry; how rain falling on it would soak in (and how much would run off); and what would happen to it if you spilled some gasoline on it? We require a system of measurable quantities for a porous medium that may be used in addressing the broadest possible set of questions about how fluids will interact with it under a variety of conditions. Porous media is characterized by a set of parameters selected to provide as concisely as possible the greatest insight into the response of the media to a range of physical processes. In this section we identify the fundamental parameters employed in describing both the media and the fluid of interest, and introduce the fundamental equations of motion for fluids in porous media.

2.1.1 Soil as a Three-Phase System

Before digging into the details of hydraulic properties and flow processes, we need to provide some background, definitions and quantitative measures. We divide natural porous media into three constituent phases: solid, liquid, and gaseous, with each of these phases consisting of an admixture of compounds.

The gaseous phase is generally dominated by the major constituents of the atmosphere; N₂, O₂, CO₂, H₂O vapor etc. Clearly, as microbes and plant roots respire, a soil’s CO₂ and CH₄ content may become significantly elevated, while in industrially contaminated sites there may be significant levels of organic vapors. By and large we keep track of the state of the gaseous phase in the soil just as we do in the atmosphere; we consider the partial pressures of each of the constituents, and keep track of the partitioning of each compound. Primary parameters of interest, then, include the gas diffusivity, and the partial pressure versus temperature relationship, and the Henry’s law coefficient (Henry’s law states that at a given temperature, for low concentrations, there is a linear relationship between the partial pressure of the gas and the concentration of the material in the liquid phase). It is noteworthy that transport of gaseous phases are dominated by molecular diffusion which generates about 2 cm/day of spreading, while molecular diffusion in the liquid phase is 100 times less effective in contaminant mass movement.† The movement of gaseous phase constituents is also a function of the changes in atmospheric pressure, which generates a pumping effect, as noted by Buckingham (1904) and Weeks (1979). Of course, gaseous phase movement can also be driven by liquid phase movement, as is the case of a saturated wetting front pushing gas down through a profile during flood irrigation. Gas phase transport has received renewed interest with the advent of soil vapor extraction (SVE) as a remediation technique for volatile organics. In this procedure, gas is pumped out of the vadose zone through a well, taking the volatile fraction of contamination with it.

Of the possible liquid phase constituents, water is by far of greatest interest. In the vadose zone it is valid to assume liquids are incompressible, since by definition liquids will only experience pressures which are negative, and generally with a magnitude of less than 10 atmospheres. In contrast, below the water table, groundwater commonly experiences extreme pressures (e.g. in cases of confined and very deep aquifers). Otherwise liquid phase properties are most often parameterized in terms of the interaction with solid and gaseous phases, such as contact angle, capillary pressure, volatility, etc., as is discussed later.

Example

What is the percent volumetric change in water in the vadose zone with a change in pressure from zero to –10 bars (10 bars is 10⁶ N/m²)? Recall that the isothermal compressibility of water, b_w, is about 4.8 x 10^-10 m²/N. The percentage change in volume of water is 100 x b_w x the change in pressure = 0.048%. Since the material is mostly gas filled when the pressure is this low, this very small change in volume will be undetectable.

The solid phase may be as complicated or simple as you like. Typically composed of largely inert minerals (rock, gravel, sand, silt and clay), solid constituents may include roots, carbon, and other highly reactive substances. Note that in this discussion we have slighted the semi-solid or soft tissue “goo” such as microbes, plants and worms that are prevalent in the vadose zone. These components are very important to understanding the chemical behavior of the unsaturated zone, and are discussed in chapter 4.

To get yourself in the proper frame of mind, pretend you are an ant. What humans consider to be uniform porous medium now looks lumpy and irregular. Although this may seem all too obvious, it is crucial to consider this problem in detail. To do this we will have to introduce the basic parameters used to conceptualize a porous media. First, let’s identify the constituents which we will concern ourselves with: the mass and volume of each phase; gas, liquid and solid. As identified in Figure 2.1a, the volume and mass of phase i are identified as V_i and M_i respectively.

2.1.2 Volume and Mass Relationships

We can now define the dry bulk density of the medium, which is the dry mass per unit volume of the medium, is denoted ρ_vb.

(1)  

Here the subscripts “V_t” denotes the total volume over which this density is calculated, “b” denotes bulk and “s” denotes the solid phase (without fluid). Typical units for density are grams per cubic centimeter, kilograms per cubic meter, or, heaven forbid, pounds per cubic foot. The upper limit on ρ_vb is the density of the pure solid mineral, often called the particle density, which we denote ρ_s.

(2)  

Notice that the value of ρ_s is independent of the value of V_s (since the value of density is taken to be a material property independent of the volume selected), whereas as we shall soon see, ρ_vb is not quite as cooperative. The solid phase density is often denoted by the specific gravity of the material, denoted by s, which is simply the ratio of the solid density to that of water ρ_w:

(3)  

Typical values of s are 2.65 for quartz (often used as the default value in lieu of further information), 2.54 for feldspar, 2.72 for calcite and 5.0 for pyrite.

Complementing these definitions, is that of porosity, or void fraction, which we denote by n_v.

(4)  

n_v indicates the fraction of the volume V which is not filled by solids. The bulk density of a medium, ρ_vb, depends primarily on two characteristics of the medium: its particle size distribution and the medium’s structure. Figure 2.1 illustrates that the size of particles alone does not dictate density: A set of media which are made up of single sizes of particles (and which are packed in the same arrangement) will have identical densities. Hence, the particle size distribution, and to a lesser degree the packing, strongly affect ρ_vb. The material structure of the medium is also significant in determining ρ_vb, particularly in the case of clays and very sharply angular aggregates.

So far so good; but what about this annoying V subscript? Can’t we just forget about the sample volume? Not really, and it is worthwhile seeing why. Let’s just run a quick “experiment” to make this clear. Imagine that you are in the middle of a large field with a tool which measures ρ_vb for any prescribed V. Start with the instrument set with V smaller than any individual grain in the soil (Figure 2.2). You will be sampling the density of either pure solid, pure void, or maybe some arbitrary mixture of the two. Now, holding the instrument steady, enlarge V continuously, all the while recording the values of ρ_vb. Figure 2.2 shows two possible realizations of such an experiment.

Although the measured density becomes “better behaved” as the control volume enlarges from the grain scale, there is no scale over which density is entirely constant. As the control volume increases, ever larger structural features of the setting are encompassed with differing properties, and thus the measured bulk density changes. The existence of such multi scale heterogeneity is widely accepted, and Dagan (1986) provides a very nice treatise on the topic. Keep this in mind: the value of all physical parameters in the vadose zone are a function of both position and sample volume.

We have two choices then: 1) throw up our arms, and figure that natural systems are hopelessly complex, or, 2) make some reasonable simplifying assumptions. For the sake of progress, we will opt for the latter option (however tempting the former). We will assume that there is some volume, much larger than the grains of the porous medium, yet smaller than the distance between dissimilar regions, which provides a representative sample of our porous medium. This volume is referred to as the Representative Elementary Volume (REV) for the system (Bear, 1972). Each extensive property of the medium is defined at a point by its value calculated taking V to be one REV about that point. We do get around to a more complete analysis of this problem in chapter 6 where we explicitly address scale dependent variability.

With that settled, we can proceed to complete our laundry list of definitions of physical properties of porous media (dropping the v subscripts since this is now implicitly understood):

Total Bulk Density : ρ_t (solid and liquid phase)

(5)  

where M_t is the total mass (solids and liquids).

Void Ratio: e (ratio of pore volume to solid volume)

(6)  

where V_a is the air volume. Note that using our earlier definitions, we may relate the void ratio to the porosity:

(7)  

2.1.3 Liquid Content

Central to most investigations of vadose zone processes is the quantification of liquid content (typically water, often referred to as the moisture content). This is expressed in either a mass or volumetric basis, as shown in 8 and 9, with the volumetric being more commonly employed, either as a water content or degree of saturation.

a. Mass basis (a.k.a., gravametric)

(8)  

b. Volume basis. θ is called the volumetric liquid content and S the degree of saturation:

(9)  

(10)  

where:

θ_o = “residual” or “irreducible” liquid content.

θ_s = saturated liquid content.

Equivalent Depth of Soil Water (Soil Water Storage): D_e = θD, where D is the soil depth increment having water content θ. D_e provides a means for expressing water content in a given soil depth increment as soil water storage (Length), and for relating aboveground water dimensions (Length; as in rainfall, irrigation, or Et) to belowground dimensions (soil water content, L³ water/L³ soil). This is very useful in water balance calculations, for example.

The typical units for moisture content are in mass per mass (for gravimetric water content, e.g., g/g) or volume per volume (for volumetric water content, e.g., cm³/cm³). Watch out though, because you will see water content in other units akin to the depth of storage; how about inches of water per foot of depth? Typically the units will indicate whether the values are on a volume or mass basis.

We need to take some time to iron out a few necessary fussy points regarding θ_s and θ_o. Lets start with the easy one: what could be more straightforward than the concept of saturation? Just fill all the pores with water, right? Well, although this is simple in principle, it is rarely achieved in the vadose zone. Due to the complexity of the pore space, it is common for some pores to remain gas filled even if there is ponded water on the surface due. This is due to trapping in dead end pores and pores that become completely surrounded by water filled pores. For this reason, several terms have been proposed to describe the state of natural soil being as wet as it will get under normal conditions. This is referred to as a “satuated”, “satiated”, or “field saturated” condition. This can be used in the definition of degree of saturation if clearly specified. Typically, the field saturated conditions render 90% of the pore space filled with fluid (Bower, 1966).

Now, for this business of residual moisture content. Unless soil is heated to over 700 ^oC , water held in hydrogen bonds will remain in place (Nagata et al., 1974, Figure 2.XX). This is clearly an extreme which will not be reached in typical natural settings, but the point is that the residual moisture content a particular medium will achieve is a function of the drying process. There are three drying processes which are widely considered in the definition of θ_o: gravity drainage, plant uptake, and atmospheric surface drying. The residual moisture content achieved under gravity drainage is referred to as “field capacity,” θ_fc, which is defined as the condition which arises after a soil has been allowed to drain to a stable condition under gravity driven conditions (assuming an infinitely deep water table, for a detailed discussion, see Assouline, S., and D. Or, 2014, The concept of field capacity revisited: Defining intrinsic static and dynamic criteria for soil internal drainage dynamics, Water Resour. Res., 50, 4787–4802, doi:10.1002/2014WR015475). From an operational perspective, this is often assumed to be the conditions which occurs when the soil water is at a pressure of -1/3 bar, although clearly such a threshold is very much a function of the soil texture, ranging from -1/30 bar for coarse sands to -1 bar for clayey soils (Assouline et al., 2014). With more time, water eventually stops draining when the films of water break between soil particles leaving isolated pockets of water.

The second commonly taken value for residual moisture content is referred to as the “permanent wilting point,” θ_pwp. You might have guessed that this is the moisture content which is observed when plants can no longer retrieve water from soil (always lower than θ_fc). Since plants vary greatly in their ability to extract water, this is operationally defined as the moisture content when the water in the medium is a pressure of -15 bar. Since water no longer flows along continuous streams at moisture contents less than θ_fc this measurement relies on vapor transport of water between pores, and so is necessarily slow in achieving equilibrium. The final method of drying is due to direct evaporation to the atmosphere. This process can bring the upper soil surface to pressures of less than -100 Bar. In summary, so long as you define the conditions over which a flow process will occur, there will be a real and measurable residual moisture content, but if you are unsure of the context where the parameter will be used, it is more reasonable to provide a range of possible values, or provide the extreme values just discussed. In fact, many soil scientists believe that θ_o should be simply taken as 0, regardless of the problem under study, as this is the only physically defensible value for arbitrary conditions. Table 2.1a and Figure 2.X provide some typical ranges for various parameters.

A primary practical importance of the field capacity and wilting point concepts is for determination of a range of plant-available soil water (PASW; Fig.2.X). Soil water storage available for plant use is generally calculated as being between these two limits (θ_FC -θ_WP), as water contents higher than θ_FC are generally not sustained for long times. It takes only one to a few days for the attainment of field capacity after wetting if conditions permit. Plant available soil water storage is an important factor in the determination of irrigation amounts for a cropped field or other soil-plant system. For practical purposes irrigation amounts in excess of a soil’s field capacity are lost to deep percolation and thus should be avoided.

A useful ‘rule of thumb’ in absence of other information is to estimate θ_FC as θ_S/2, and θ_WP as θ_FC/2; i.e., a soil exhibiting this property will have lost one-half it’s saturated water content at field capacity, and another one-half the remaining water by wilting point.

Knowledge of soil water content is also important to geotechnical applications where the attainment of high bulk density by compaction minimizes future settlement and increases the soil strength. For each level of compactive effort there is an “optimal” water content at which maximum compaction is attained (this occurs at an approximate degree of saturation of 0.8). Hence, knowledge of the soil water content can assist a geotechnical engineer in determining how much water should be applied to a layer of earth fill to achieve optimal compaction. Such information is also useful for farmers who want to avoid excessive compaction of their fields by farm implements (see Hillel 1980, Chapter 14 for details).

Table 2.1a Typical values of physical properties for several soils.

2.4       The Energy State of Soil Water

Imagine a column of saturated soil that is allowed to drain: water will flow downward draining the upper soil, but why?  What drives this motion? Water content alone is not sufficient to characterize the water status of soils. An additional critical attribute is the energy state of soil water. Soil water contains different forms and quantities of energy, as does any body of matter. Some familiar forms of energy are kinetic, potential, thermal, electrostatic, and chemical. The last aspect of characterizing the state of water in porous media requires the introduction of the energetic status of the liquid in the porous medium. While all of these quantities represent energy, we will  discover that most liquid movement in porous media results from gradient in only the mechanical potentials (as opposed to electrostatic, chemical or thermal potential). In principle this is as simple as water flowing downhill in a stream: water is flowing from higher elevation (thus higher gravitational potential) to lower elevation under constant pressure (the gauge pressure of water at the air-water surface is zero), thus moving from areas of high gravitational potential to areas of lower potential with constant pressure potential.  But given the differing roles of the different potentials, we must be deliberate in identifying these terms and their role in soil-water processes.

Let’s start with the potential that you experience every time you move: Kinetic energy (E_k) is acquired by virtue of motion and is proportional to the velocity squared . Illustrating how different the primary processes are in natural porous media compared to daily life, because the movement of water in soils is generally quite slow (e.g., of the order of meters per year), its kinetic energy is almost always negligible. We call these slow-moving flows “creeping flow” to indicate that the forces of acceleration associated with changing kinetic energy are much smaller than the forces generated by viscosity and gradients in pressure.  These flows are associated with Reynolds numbers less than 1.  For every general rule there are special case exceptions, and in soils this can here occur during extreme events: flow along macropores can have velocities with kinetic energy that are appreciable. We will discuss this in detail in the context of dynamics of fluid in porous media, but for the moment we simply posit that potential energy is the primary driver of flow outside of these rather extreme cases.

With respect to potential energy we need to discriminate between terms which are mechanical from those which are electrostatic, chemical or thermal. The point here is that although hot water has more energy than cold water, this energetic difference does not cause the hot water to flow toward the cold water if they are in physical contact (think of a lake in summer with warm water on top and cool water below) . Similarly, a pure water has higher osmotic potential than does a saline water, but when you have salty water below fresh water you will not see the salty water flow in bulk toward the fresh. So how can these chemical or thermal potentials give rise to movement of water?  To translate chemical potential into mechanical potential requires there be an interface which selectively constrains the motion of selected components of the system (i.e., an charge barrier, permeable membrane, or gas-liquid interfaces). Some examples are helpful here. Consider a semi-permeable membranes that allows water to pass but excluding solutes.  This can result in water delivery to the side with higher solute concentration, resulting in bulk water movement.  At a clay surface with electrostatic bound cations population, when exposed to water, these cations will hydrate, pulling in bulk water into the clay.  Systems with such selective regions can develop huge mechanical pressures (e.g., turgor in cells or the pressure drop across reverse-osmosis desalinization systems). Such interfaces occur in the natural world at roots, cell walls, and air-water interfaces (Figure 2.X123). While very important to life, these processes are somewhat esoteric in the case of the movement of bulk water in soils.

To predict fluid movement in the vadose zone we are generally concerned with keeping track of the total mechanical potential (TMP which we will indicate with the letter H) of the fluid phase. The total mechanical potential is made up of the sum of the gravitational potential (due to elevation) and the pressure potential (the force per unit  area exerted by the fluid on its boundary – Figure 2.X124). As we introduced above, fluid motion results from a gradient in the TMP.  So we are concerned with gradients in the TMP.  Still we need to define reference (or datum) in order to provide a frame of reference  from which the TMP might vary.  For pressure we typically define the atmospheric pressure as the datum (by assigning it a value of 0). This is the pressure of the water table.  In selecting a datum for gravity, we choose an elevation which is convenient for each problem. In the lab this may be the elevation of the surface a pool of water, and in the field it is often taken to be the soil surface. We emphasize that the specific choice of datum is immaterial: fluid moves relative to local changes (gradients) in potential, so any constant added or subtracted from the total potential falls out when we take the derivative required to compute the flux.

So how are we going to add elevation to pressure to compute the TMP?  We need to employ common units, and the most natural approach is to use the units of elevation: length. We typically denote the elevation with the variable z which becomes more positive in the upward direction. The pressure potential, however, most naturally has units of force per area(see Table 2.2), and can be calculated within a hydrostatic liquid by multiplying the distance to the water table (say, h) by the liquid density and the acceleration of gravity: . Here we have the trick which makes flow calculations work out: we can express the pressure in terms of the equivalent height of a column of water which would generate this pressure. This is the pressure “head,” which we denote h, and has units exactly compatible with elevation. We denote the TMP as H and is equal to the sum of h and z.

To get a feel for how this total mechanical potential works, let’s consider the simplest possible case: a tank of water at rest (Figure 2.X124). Since there is no flow in the tank, the condition which is called “hydrostatic,” we can be sure that the total mechanical potential is the same everywhere (remember water would flow towards lower potential, but we posited that the water was at rest). As mentioned, we may define our datum for pressure and elevation arbitrarily: here we chose the top of the water surface as the zero pressure (denoted with the letter h when in units of head) datum and the bottom of the tank as the zero elevation datum. Clearly elevation increases with height in the tank, while pressure decreases such that the total mechanical, denoted with the letter H when measured in terms of head, potential is constant, consistent with the lack of flow in the tank. Clearly one may convert between units of pressure and equivalent water column height (Table 2.2). A word of caution is in order regarding the assumed density of the fluid which fills the pore space: what about cases of saline solutions, or non-aqueous phase liquids (NAPLs; hydrocarbons etc.)? We must be consistent in our units of pressure measurement, and in this case we will stick to height a column of pure water would be lifted (at 20 deg C) regardless of the fluid under consideration.

Table 2.2 Conversions between common units of pressure. To convert the units given in the left column to that listed on the upper row, multiply by the conversion factor listed.

People often wonder about pressures more negative than -1 atmosphere: isn’t this somehow non-physical? We know from our kitchens that boiling occurs when the vapor pressure of the water reaches the atmospheric conditions, so why doesn’t the pore water boil? Don’t feel bad if you came up with these questions; in contemporary literature some authors still wondered if water pressures more negative than -1 ATM were physically possible (Gray and Hassanizadeh, 1991; Miller, 1994; Gray and Hassanizadeh, 1994). However, experiments have shown that liquid water can sustain negative pressures of up to 150 bars before vaporizing at room temperature since the vapor pressure of water is actually much less than one atmosphere at this temperature (e.g., Zheng et al., 1991) supporting Miller’s position. This conclusion is confirmed in trees over 10 m in height, where the pore pressures more negative than -1 bar are essential to the delivery of water to their leaves.  Stable water at such negative pressures only occur in very small bodies of water (e.g., xylem of plants), but at such great negative pressures water will only be found in very small regions, as is discussed below. Thus a few points are clear: (a) negative pressures less than -1 bar arise commonly in porous media due to the forces listed in Table 2.1, and (b) it is not unreasonable to consider the fluid-dynamic behavior of water at pressures more negative than -1 bar.

Example 2-X:  Gravitational Potential

Figure XXX. Self-study of the elevation potential of a soil column to understand the role of selection of datum.  All distance values shown are in mm.

Case 1:

Point A is 100 mm above the datum and point B is 200 mm below the datum. Find the difference in gravitational potential: between the points.

Solution:

and , so .

Case 2:

For the same two points A and B we now change the elevation of the reference to 100 mm above point A.

Solution:

and , so .

Notice we are not here dealing with the possible difference in pressure potential at points A and B, just their elevation potential.  The key point here is that the choice of datum does not change the value of the difference in potential between two points.

When do we employ total potential, versus mechanical potential, in the estimation of fluxes from high to low potential? As we discussed above, the osmotic potential can cause movement of water when there is an interface in the system which constrains the flux of selected constituents. For example, if there were two side-by-side cups of water solutions in a sealed box which could exchange water only in the vapor phase, then the air-water interfaces would act as constraints on the non-volatile solutes, while letting water be freely exchanged.  If left long enough the water would move in the vapor phase until the water came to have the same total potential in the two cups. For instance, if one was pure water, and the other had a dissolved salt, then the pure water would tend to lose water in the vapor phase which would then condense into the saline solution, diluting the salt, and bringing the total potential of the two vessels closer to equal. Another important example is that of plant turgor: in contact with pure water, plants can obtain high cell pressures by having high dissolved concentrations in their tissue (low potential), so that the equilibrium total potential is obtained when the internal pressure in the tissue is high. Similarly, if the salt concentration in soil is too high, plants cannot create an energetic gradient to drive the saline water into the plant, and thus can wilt.

You will often hear of water described to be under “tension,” “matric potential” or “suction.” These terms all refer to the same quantity, the strength of attraction of water to a medium.This term is not rigorously tied to a particular physical mechanism, though some authors employ this term to identify the energy of interaction between water and solid surfaces (recall Figure 2.X illustrating the high energy of this interaction by way of the temperatures exceeding 700 C to remove all of this water – for a classic treatment of this see Brown and Montgomery, 1913, “Dehydration of Clays,” Bureau of Standards Technologic Paper no. 21, Washington, D.C.: https://nvlpubs.nist.gov/nistpubs/nbstechnologic/nbstechnologicpaperT21.pdf). While this is a perfectly fair energy term to quantify, it is important to recognize that, in general, solid surfaces will be at some intermediate state of water absorption, so any theoretical maximum value is rarely relevant. Generally speaking the energy per unit area of this interface is far in excess of the air-water surface tension, which theoretically should lead water to spread across all such surfaces.

Hydrostatics in Porous Media

Having defined the parameterization of the porous media, and the units of the liquid phase potential, we can now discuss the interactions between soil and fluids. In principle there could be any number of fluids in these spaces, say m, each of them filling some partial volume of the pore space; V₁, V₂, …Vm. The only constraint on this system is that

(19)  

where we recall that n is the porosity. Fluids of primary interest are air and water; however other gases and non-aqueous phase liquids (e.g., gasoline and diesel in both liquid and gaseous forms) are also of significant interest. The focus of this text are on water and air, but the formulations are sufficiently general to allow direct extension to other fluids.

Why doesn’t all the water just fall out of the porous media? There are at least four groups of forces contributing to the retention of water in soil, which are enumerated here in order of decreasing strength. Each of these forces has a connection to the structure of water, which has a polar nature, slightly negatively charged on the side where oxygen projects, and positively charged on the side where hydrogen is prominent.  This means that water molecules have a highly dipolar electrostatic structure.  Water molecules then develop a hydrogen bond between neighbors between their oxygen and hydrogen.  Among other things, this leads to high bonding energy between water molecules, which leads to high surface tension, high latent heat of vaporization, and high solubility of ions.  To summarize:

1. The hydrogen of water is attracted to the negative surface charge of mineral surfaces (Van der Waals attraction).
2. The periodic structure of the clay surfaces gives rise to an electrostatic dipole which results in an attractive force to the water dipole.
3. Osmotic force, caused by ionic concentration near charged surfaces, hold water.
4. Surface tension at water/air/solid interfaces maintains macroscopic units of water in pore spaces.

These forces are discussed in detail in Sposito (1989) and Kirkham and Powers (1972) for readers who would like to learn more about these fascinating and fundamental relations. Attention in this text is focused on transport phenomena in the vadose zone. The first three forces listed in Table 2.1 act only at very short range, and essentially immobilize water. In contrast, surface tension produces effects which affect water in bulk, and is therefore more influential in transport.

Surface Tension

Surface tension is a measure of the energy required to generate an interface between two phases.  Here we will be most concerned with the air-water interface, as shown in Figure 2.X2, but the concept applies to solid-gas interfaces, and liquid-solid interfaces.  An intuitive understanding of surface tension can be obtained from a simple thought experiment. Imagine a block of water in a container that can be split in two (Figure 2.X3). For convenience think of each water molecule as being organized in a cubic lattice structure, with six bonds. Now we quickly split this block of water into two halves. The molecules on the new air/water surfaces are bound to only five of their neighbors. It took energy to break these bonds, which will be manifest as a free surface energy. Since the water surface has a constant number of molecules on its surface per unit area, the energy required to create these surfaces is directly related to the surface area created. This would suggest that surface tension should have units of energy per unit area, or force per unit length.

It is easy to measure surface tension using the sliding wire set-up shown in Figure 2.X4. If the force measured on the scale is F, and the width of the film is L, then the surface tension, σ, is simply

(20)  

Wait a minute! How did that factor of 2 sneak into Equation 20? Simple: since the film is open on two sides, we are creating two air/water interfaces, and thus the effective length of a single-sided free surface is 2L. Note that you could also keep the entire problem in terms of energy per unit area in keeping with the thought experiment of the previous paragraph: carefully record the product of force required to stretch the film multiplied by the distance over which that force is exerted (to obtain the energy input) and divide by the area (l x L). Table 2.3 lists some typical values of surface tension for fluid pairs.

Table 2.3 Some values of surface tensions for various interfaces (taken from the CRC Handbook of Chem. and Physics 68th ed., 1987).

Classical Derivation of Laplace’s Equation

We are concerned with how surface tension manifests itself in a porous medium. In particular, we would like to know the static fluid pressures which result from surface tension acting on fluid surfaces, which are curved within a porous medium. Consider the infinitesimal curved fluid surface which is shown in Figure 2.5. For the sake of generality, we allow for curvature along two axes of radius r₁ and r₂. To determine the pressure p acting on the surface we simply add up the forces acting on the interface in the vertical direction. Since the system is at equilibrium, these forces must add to zero. The forces acting upward (which we will take as positive) result from p₂ and r₂

(21)  

where σ is the surface tension of the interface. Now since this is a very small patch, the angle δθ₂ is very small, and we may use the small angle approximation for sine:

(22)  

By symmetry, the expression for F_down is identical to Equation 21, substituting appropriate subscripts and signs. The net force on the surface is then given by

(23)  

which, canceling δθ₁δθ₂ from all terms, simplifies to

(24)  

or

(25)  

which is Laplace’s Equation for pressure across an interface at a saddle point (an anticlastic surface). Notice that if both radii were of the same sign and magnitude (i.e., the surface was spherical) we have r₁ = – r₂ = R in Figure 2.5, and Equation 25 would reduce to

(26)  

which is the familiar form of Laplace’s equation for pressure across a spherical fluid interface. Of course this is the same expression one obtains for the pressure behind the gas/liquid interface in a capillary tube of radius R. For the general case where , we define the reduced radius of curvature, R, as

(27)  

or

(28)  

where we have used a sign convention for r₁ and r₂ such that they are taken to be positive if measured in the wetting fluid and negative if measured in the non-wetting fluid (typically air).  Note that an element of water in a porous media often displays complex curved interfaces.  Even the simplest case of two spheres in contact with water between has two characteristic radii, as shown in Figure 2.X

Introduction to the of Fluid Interfaces

The derivation of Laplace’s equation provided in section 2.3.2.1 is completely acceptable, but in this case there is good reason to look at this problem from another perspective, that of thermodynamics, to gain further insight into the workings of fluid interfaces. In the following sections we shall re-derive Laplace’s equation, using a somewhat less intuitive and more formal approach. This will be a warm-up, for then we will push thermodynamics a bit further to obtain some very handy additional results.

The Thermodynamics of Flat Interfaces

Consider a small, partially filled container of water with a flat water-air interface. The system has water in three phases: liquid, gaseous and as a taut interfacial surface, denoted by subscripts ‘l’, ‘g’, and ‘a’, respectively. We will denote the surface area as Σ and the surface tension at the air/water interface as σ. Since there is no barrier to water moving between phases, the chemical potential of the water in each phase is the same, i.e.,

(29)  

The subscript ∞ is used to indicate values of intensive parameters which are constant in all three phases (recall intensive parameters are those which do not depend upon the amount of material in question, in contrast to extensive parameters which do depend on the amount, or ‘extent’ of the material). We will consider a system at constant temperature, T. We will use N_i to denote the number of atoms, S_i the entropy, E_i the internal energy, P_i the pressure and V_i the volume in phase i.

Considering the gas and liquid phases, we have three fundamental differential forms

(30)  

(31)  

(32)  

The total energy and entropy for the system are simply the sum of the components, i.e.

(33)  

(34)  

Since the surface is two-dimensional, we take the number of atoms in the surface to be zero in comparison to the atoms in the three-dimensional volumes of gas and liquid. So we may write

(35)  

Thus the fundamental differential form for the system is given by

(36)  

Employing the Gibbs-Duhem relationship (i.e., for an exact differential, the differentiation may be shifted from the extensive to intensive variables maintaining equality) we find

(37)  

or

(38)  

which is the equation of state for the surface phase (analogous to the equation of state for an ideal gas ), and relates the temperature dependence of surface tension to the magnitude of the entropy of the surface.

Laplace’s Equation from a sphere in Space

Now let’s consider the effect of a curved air-water interface. Now P_g and P_l are not equal, while μ_g and μ_l remain equal. Considering a droplet floating in gas filled space of constant volume, the fundamental differential form for the system is given by

(39)  

Considering a very small transfer of volume, δV, between the two phases, we see the chemical potential terms would be equal and opposite, that the total change in energy in the system is unchanged (we are doing no work on the system), and the entropy would be constant to first order. Since the volume of the system is held constant from 39 we find

(40)  

or defining P_d as the pressure difference across the interface of the droplet, we have

(41)  

We can calculate the differential on the right-hand side by noting that for a sphere and , so , hence

(42)  

which is Laplace’s equation for the pressure across a curved interface where the two characteristic radii are equal (see 26).

Although you probably have had it with this equation by now, it is only fair to show the simple way to obtain this result. Consider the pressure balance across the middle of the droplet. We know the surface tension of the water about the center of the droplet must equal the pressure exerted across the area of the droplet by the liquid. The area of the droplet at its midpoint is πr² at pressure P_d, while the length of surface applying this pressure is 2πr at tension σ. Equating these we obtain

(43)  

as expected.

Across Curved Interfaces

The curved interface also affects the vapor pressure of the liquid species. Suppose we have our spherical water droplet in a fixed volume at equilibrium with the vapor phase. The chemical potential of the water in these two phases will be equal, and will remain equal through any reversible process, i.e.,

(44)  

and

(45)  

Writing the fundamental differential forms for the gas and liquid phases for a system with one mole of each species (following equations 30 and 31) we have

(46)  

(47)  

As before, the Gibbs-Duhem relation of the fundamental differential forms follows the principle that for an exact differential form, equality is maintained if differentiation for each term is exchanged between intensive and extensive variables, hence 46 and ?? may be written

(48)  

(49)  

Dividing 48 by N_g and 49 by N_l we get the equations in a per mole basis. By further assuming that temperature and surface tension remain constant, we obtain

(50)  

(51)  

where v in lower case indicates the volume per mole. Now employing ?? we find

(52)  

which may be written

(53)  

Using Laplace’s equation (e.g., 43) this may be re-written as

(54)  

or

(55)  

We note that v_l is approximately four orders of magnitude less than v_g, so we will suppose that can be accurately approximated by , and we obtain

(56)  

Now assuming the vapor phase acts like an ideal gas, , whence 56 becomes

(57)  

This may be integrated from a flat interface () to that with radius to obtain

(58)  

where P_∞ is the vapor pressure of water at temperature T (P_g/P_∞ is the relative humidity of the soil air system). Taking R̄ as the specific gas constant for water (i.e., R̄ = R/v_l), and noting that the left-hand side of 58 is just P_d, this may also be written

(59)  

where P_d is the pressure potential of the water. This very useful relation is known as the psychrometric equation, which provides the relationship between the humidity of the air in equilibrium with an air-water interface and the potential of the water.  This allows the determination of very high negative potentials through measurement of the vapor pressure (or relative humidity) of water in porous media. For instance, at a pressure of -15 Bar (the permanent wilting point of many plants), P_g/P_∞ is 0.99. Measurement of P_g/P_∞ may be made in a number of ways. Most commonly a small thermocouple is cooled by application of a reverse current while its temperature is read with a second thermocouple. As the temperature of the first thermocouple attains the dew point of the vapor, its temperature decline is sharply reduced due to the energy of condensation of water on the device. Knowing the dew point, it is straightforward to obtain the relative humidity with good precision (Rawlins and Campbell, 1986).

Temperature Dependence of Surface Tension

It is often overlooked that all the measurements we take regarding water/media interactions are strongly temperature dependent. In many of the equations derived above we have treated surface tension as a constant property of a fluid. This is certainly not the case. The temperature dependence of air-water surface tension is illustrated below in Figure 2.6 where we see that for this pair, surface tension decreases at approximately one percent per 4℃, which could give rise to discrepancies of greater than 5% between common laboratory and field conditions, or even through diurnal fluctuations.

The Characteristic Curves

Now for the real excitement: putting the solid and fluid together! To do this, we will undertake yet another simple thought experiment. Take a funnel, and seal a flat porous plate in it to make a Buchner funnel, as shown in Figure 2.7. Select the porous plate so that its pores, of radius r, are much smaller than the pores in the soil sample. Now attach a tube to the bottom of the funnel, and fill the region below the porous plate to the end of the tube with water. This tube is fed from a beaker on a precision balance, which keeps track of the amount of water in the soil sample. By wetting the porous plate, we know from Equation 26 that air cannot enter the porous plate until a suction greater than P_entry

(60)  

is applied to the plate. Since the porous plate has pore radii smaller than the pore in the soil sample, we know that the plate will remain saturated as we draw the water from the sample. We now place a thin piece of dry soil on the porous funnel surface. The experiment will take place over three stages: (1) wetting, (2) drying, and (3) re-wetting. For stage (1), we incrementally bring the end of the tube up, finally attaining the elevation of the soil. At each lift we measure the amount of water drawn up by the soil, with pressures ranging from P_entry to 0. Each measurement is taken allowing the system to come to a steady state. Measuring elevation, H, as positive upward, the pressure applied to the water in this soil will be given by

(61)  

where ρ_w is the density of water and h is the pressure head of the water in the sample, which we will assume is much larger than the thickness of the sample l. In stage (2) the end of the tube is incrementally lowered to apply a negative head to the water in the soil. Finally, in stage (3) we again lift the beaker to an elevation slightly below the height of the soil sample keeping track of the amount of water taken up by the sample. A typical plot illustrating how the soil sample’s volumetric water content (equation \ref{eq:volumetric_liquid}) changes in response to pressure changes for this type of experiment is shown in Figure 2.8, with a few more lifts and drops thrown in for good measure. The vertical axis in Figure 2.8 represents the pressure in the sample, as reflected by the hanging height of water, h. The horizontal axis shows the soil’s volumetric moisture content.

The equilibrium relationship between pressure and moisture content is called the soil water characteristic (SWC). The SWC is an important soil property related to the distribution of pore space (sizes, interconnectedness), which is strongly affected by texture and structure, as well as related factors including organic matter content. The SWC is a primary hydraulic property required for modeling water flow, for irrigation management, and for many additional applications related to managing or predicting water behavior in the porous system. A SWC is a highly nonlinear function and is relatively difficult to obtain accurately. Because the pressure potential extends over several orders of magnitudes for the range of water contents encountered in the field, the potential is often plotted on a logarithmic scale (Figure 2.X).

This simple experiment illustrates most of the physics which control the retention and movement of fluids through porous media, so it is worth considering these results carefully. First consider the wetting curve alone, spanning the points labeled (1), (2) and (3) in Figure (2.8). We see that more water is taken up by the soil (i.e., the moisture content increases) as the beaker comes closer to the elevation of the soil (i.e., as the pressure of the feed water increases).

One can gain a basic understanding of the characteristic S-shape of this curve by relating it to the soil particle and pore size distributions. Consider a soil made up of particles with size distributed between d_min and d_max as shown in Figure 2.9. By passing the sample through successive sieves, from very fine to coarse, we may evaluate the percentage of the grains which have size greater than any given size, as plotted in Figure (2.9a) (see section 1.3.2). Without stretching your imagination too far, it is plausible that the pore size distribution in the original soil is quite similar to the particle size distribution. Now the ordinate goes from mass of particles, to volume of pores. The plot of pore size distribution which is obtained is quite similar to that of particle size, as shown in Figure (2.9b).

This is often conceptualized in terms of the soil being made of many differing pore sizes classes which could be represented as equivalent parallel capillary tubes of many radii (Figure 2.X).

Figure 2.X Capillary tube conceptual models for the pores in soils. (a) illustrates the notion that the pores of larger diameter draw water to lesser heights, while fine pores draw water higher. Thus, a coarse soil will draw in water with less pressure than will a fine soil. (b) considers a single vertical soil sample at equilibrium with water provided to the lower boundary, with a variety of pore sizes represented by capillary tubes. This shows that all pores will be filled at 0 pressure, then as the the pressure potential becomes more negative with increasing elevation (recall that at equilibrium the total potential is constant in the sample), only ever smaller pores can retain water.

(a)	(b)

Using Equation 26, we can equate pore size to the capillary filling pressure of each pore, and re-plot the pore size distribution as a filling pressure plot; Figure (2.9c). Notice the volume of pores which are filled in Figure (2.9c) is simply the degree of saturation. At last, in Figure (2.9d), we have obtained the wetting curve from the particle size distribution. The idea here is to point out the connection between the wetting characteristic curve and the physical properties of the medium. A similar argument is often used when trying to determine how a medium will drain using the particle size distribution (e.g., Arya and Paris, 1981).

Now let’s follow the process of draining, going from points 3 to 4 and 5 in Figure (2.8). Why isn’t this curve following right along the wetting curve? To understand this we shall consider the filling and emptying of a single pore, as shown in Figure (2.10). The idea represented in this figure again follows directly from Equation 26, and were introduced by Haines in 1930.

The wetting process for this pore is illustrated in stages 1 through 3 of Figure 2.8 and (b) of Figure 2.10, while the draining process is illustrated in stages 4 through 7 of Figure 2.8 and (c) and (d) of Figure 2.10. For any given pressure, there are two states that this pore may be in, depending on whether it is in the process of filling or draining. In this sense the saturation of the pore is a function of both the pressure and the history of the pore. For this reason, the wetting and drying curves for soil are referred to as “hysteretic.”

Figure 2.10 shows that the filling process is held up until the pressure at the soil is sufficient to fill the largest pores in the medium. On the other hand, once the medium is filled, it cannot drain significantly until the pressure becomes sufficiently negative to pull the water out of the pore necks that lead into the medium. Clearly between the pressures of -2σ/r₂ and -2σ/r₁ the moisture content is a function of not only the pore pressure, but what the history of wetting has been.

Based on the ideas contained in Figure (2.9), h_w can be identified as the pressure at which the largest group of pore bodies fill. From Figure (2.9) and Equation 26, this indicates that the largest pores in the medium are of radius r_max = 2σ/h_w. Similarly, the pressure h_a is dictated by the size of the typical pore throats, which have a radius of approximately 2σ/h_a.

Why doesn’t the soil drain completely, and have point #5 in Figure (2.8) return to the dry state? There are two phenomena which account for this behavior, with relative magnitude depending on the chemical and physical properties of the medium. First, any water which entered into chemical bonds with charged surfaces in the medium cannot be removed through simple fluid tension: the bonds holding them are as strong as surface tension, and so they cannot be ripped away by the bulk fluid. Furthermore, even in an inert porous medium, fluid can be held in the very small radius regions at particle contacts. These isolated donut shaped pockets of water, referred to as “pendular,” remain in place, no longer hydraulically connected to the fluid reservoir (as shown in Figure 2.10 (d)).

Contact Angle; Sources and Effects

The main hysteretic mechanism illustrated in Figure 2.10 is referred to as the “ink bottle effect,” due to the similarity of the behavior of bottles of ink with very narrow spouts. Another source of hysteretic behavior is hysteresis in the contact angle. Although contact angle was mentioned before in passing, it is now time to take a closer look at this interesting corner of fluid behavior.

Figure 2.11 illustrates a cross-section of a straight-line contact between a fluid and a solid. What should the angle of contact between the solid and fluid be, and why? Well, since the system is at equilibrium, we should be able to balance the forces at the point of contact. Considering only the component of forces acting in the horizontal direction, we find

(62)  

where F_sg is the solid-gas surface force per unit length, F_sl is the solid-liquid force, F_lg is the liquid-gas surface force and γ is the contact angle. Per unit length, F_ij = σ_ij, so this may be put in terms of the relative surface tensions

(63)  

which is the Youngs-Laplace equation; the relationship between the contact angle and the three surface tensions. Solving [2.56] for the contact angle we find

(64)  

Care needs to be taken when applying these results, as there are some physical limits on possible values of γ: the contact angle is bounded by . So if the operand of cos^-1 is greater than 1, then g will be 0 ^o, while if the value is less than -1, the value will be 180^o. Since it is often true that for water, the contact angle for water going into geologic material is often taken to be 0 ^o.

Figure 2.XX: Illustration of wetting and non-wetting surfaces. (a) Wettable silt soil surface (γ ~ 0^o) (b) Treated water-repellant silt soil (γ ~ 70^o) (Bachmann et al., 2000) ]

Now how does this contact angle relate to hysteresis? As shown in Figure 2.12, the contact angle often exhibits a different value in the advancing and receding cases. This is the reason a drop of water on a flat plate will not start to move as soon as you tilt the plate: more energy is required to remove the water from the trailing edge of the plate than is given up by the sum of the gravitational potential plus the energy released wetting the plate. This is sometimes referred to as the rain-drop effect, perhaps referring to how a raindrop moves down a window pane. The relevance of the rain-drop effect to capillary hysteresis is simply an extension of the observation regarding the plate and drop: a medium will retain water more vigorously than it will absorb water (Bartell and Wooley, 1933).

That this rain-drop effect exists is well known to any child who has bemoaned a rainy day at the window (without it all the droplets would run off the pane), but what is the physical basis? Consider the effect which might result from a scratch on a solid surface as shown in Figure 2.13. At the microscopic level the Youngs-Laplace equation is adhered to, but from a macroscopic point of view, the drip cannot advance until the apparent contact angle is quite large. Similarly, upon retreat the macroscopic contact angle will be much smaller than the true microscopic magnitude (Figure 2.13). The same pair of effects would result from a smudge of oil or other contamination which caused a local reduction in the microscopic contact angle.

A fine cylindrical tube (capillary tube) is often taken as the simplest model for capillary attraction. Given its prominence in describing capillary processes, and because it is a useful conceptual model, we now discuss this system. Figure 2.14 shows two possible capillary tubes. The tube on the left has a zero-contact angle between the tube and the liquid, while the tube on the right has a non-zero contact angle.

So why does the water go up into the tubes? Because energy is liberated when this liquid coats solid surfaces. The ultimate height of the capillary surface is that which balances the energy required to pull the liquid up against gravity and the energy liberated when the liquid coats the solid surfaces. Thus, the derivative of the gravitational energy and surface tension terms are equal. The gravitational (E_G) and surface (E_S) energies of the system are

(65)  

and

(66)  

taking the reservoir surface as a datum, for liquid of density ρ. A stable elevation demands that the rate of change in these equal and opposite with change in height, so we find

or

(67)  

Which is quite akin to our good friend the Laplace equation, but only the same if . The tricky part here, which is often unappreciated, is that the gas-liquid interfacial energy does not show up in this result, nor does the contact angle. Many people write 67 in terms of the gas-liquid interfacial energy and the contact angle using ?? but this really obscures the important physical nature of the process, as pointed out by Selker and Schroth (1998).

We should look more closely at the case of the zero contact angle, which contains two special cases: (1) and (2) . In case (1), the result of 67 is applicable. But in the case (2), when the surface energy of the gas-liquid interface is relatively small, the liquid will spontaneously spread in a thin film over the entire inner surface of the capillary tube. Now the capillary rise of the liquid in the tube is simply a matter of the pressure which is generated by a curved film of water that is tangential to the walls of the tube. In this case the radius of the curved surface is the radius of the tube wall and the pressure is given by Laplace’s equation for a spherical surface 42).

The Independent Domain Model of Hysteresis

If you confine yourself to the conceptual model presented in Figure (2.10), you will always be on either the main wetting or main draining curve, but Figure (2.8) shows that if one stops draining at some point short of θ_r, θ{Ψ} will depart from both the main wetting and main draining characteristic curves. This can be understood and accounted for through a simple extension of Figure (2.10). We noted that the key features of filling and draining of a single pore had to do with that pore’s neck and body radii, which we will refer to as r_n and r_b respectively. With this in mind, consider categorizing a medium by neck and body pore radius. As before, we will immediately translate these radii into equivalent capillary pressures via equation 26. So the idea is that we can characterize any medium as the aggregation of independent pores, each of which has a characteristic filling pressure (h_f) dictated by the body radius and emptying pressure (h_e) controlled by the neck radius. [Although we appeal to equation 26 as a motivation, we can simply measure the filling and emptying pressures, so this is really just a crutch for visualizing the system.] Finally, we assume each pore is hydraulically connected to the bulk medium so that if a pressure is established at one of the medium’s boundaries, all the pores will experience that pressure. Each pore responds to changes in pressure independently (without regard to the state of the surrounding pores). This perspective was brought to vadose hydrology by A. Poulovassilis (1962).

The essence of this assumption is that the filling and draining of each pore is determined strictly by that pore’s geometry, regardless of the connection of that pore to surrounding pores (thus the term “independent”). This assumption is certainly fine for the wet end of the characteristic curve, when water fills most pores, but when the medium is quite dry it is possible the pores will become isolated hydraulically by being surrounded by empty pores. In this case the independent domain assumption breaks down for the isolated pore since it cannot drain regardless of the level of pressure exerted on the boundaries of the sample. We may expect then that the independent domain assumption will not lead to good predictions of residual moisture content, which must be carefully defined based on the process of interest, as discussed in section 2.2.

Pores will have a range of magnitudes for he and hf which can be described by a joint probability density function, as plotted in Figure 2.16. Note that the distribution functions for h_e and h_f are not independent, since pores with very small necks are more likely to have similarly small bodies. The total volume under the curve in Figure 2.16 is 1, corresponding to the fact that all pores will have some combination of the two characteristic radii, which is stated mathematically in equation 68.

(68)  

Pore necks are, by definition, smaller than pore bodies, so the absolute value of he is necessarily larger than h_f for a given pore. So we may draw a line at 45° on the h_e-h_f plane and know that all pores will lay below that line, as indicated in Figures 2.16 and 2.17.

To understand the density function plot more intuitively, it is useful to consider the physical interpretation of several of the points indicated in Figure 2.17. You see that the end-points of the pressure scales are defined by the largest and smallest bodies and necks. The appeal to Laplace’s equation here is primarily heuristic: the exact radii identified are not what you would find under a microscope, but provide a useful conceptual model. If we consider some horizontal line with a value h1 going from the right boundary (h_emax) to the 45° line, this delineates all of the pores which fill at pressure h1. Since the filling pressure is determined by the body radius, we can consider these to be the pores with some body radius r_b such that h_f=2σ/r₁. Similarly, the vertical line shown in Figure 2.17 at pressure he defines all the pores which empty at pressure he. Realizing that the pore neck diameter determined the draining pressure, we can again associate all of these pore as having neck diameters of rn such that h_e=2σ/r_n. The intersection of these two lines corresponds to the pores which fill at h_f and empty at he, and the height of f at this point indicates the fraction of pores which have this particular set of characteristic radii.

With this understanding of how to organize filling and emptying of pores, it is straightforward to obtain the various characteristic curves (once you have determined the joint density function). Applying this method boils down to figuring out the range over which to integrate the density function. As an example, consider the series of characteristic curves shown in Figure 2.8 within the context of the independent domain model. Starting from point 1 where the medium is dry, all the pores are empty. So now we will go from h = -∞ to h = 0. Through this sequence of pressures we want to add up all the pores which fill regardless of their emptying pressure. Thus we sweep downward through hf adding up all the pores with emptying pressures between h to -∞, integrating the density function as we go. This first wetting process is illustrated in Figures 2.18a and 2.18b.

Having saturated the sample, to obtain the main draining curve we follow the same procedure. This time the pressure starts at 0 and becomes more negative. Since all the pores are initially filled, we need to integrate along vertical lines all the pores which empty at a given pressure regardless of the pressure at which they filled (Figure 2.18c).

Things get interesting when we stop prior to emptying all the pores, and turn around and start to re-fill the medium. Pores which are already filled with water cannot be refilled, which is the only feature which has changed from the main filling procedure. We still sweep down with a horizontal line, but now we only add the pores which fill between pressures h₁ and h_f (Figure 2.18d). Figure 2.18e continues this sequence, now following the draining process. So once the joint density function is known, it is simply a matter of determining which area of the domain to integrate over to determine the moisture content for any sequence of pressure changes. This procedure is quite amenable to numerical implementation.

This procedure can be stated mathematically as integrating over the domain of filled pores. We set up a few of the integrals which would be evaluated in Figure 2.8 as examples. In the main wetting process (Figure 2.18a), we simply carry out the integral

(69)  

Supposing we take the system to saturation, and now start draining the medium to obtain the main draining curve. Having filled all the pores to start with, we must calculate how much has been drained, and subtract that from 1:

(70)  

Suppose now that we stop draining at a pressure h₁ and start refilling (Figure 2.18d). For convenience, let’s define the moisture content at the completion of the draining as

(71)  

In refilling, only the pores which were empty should be counted. So the integral will be restricted to the range of emptying pressures from h_f to h₆:

(72)  

The only trick here is to keep careful track of the range of integration, with some bookkeeping: you have to keep track of the turning points of the soil’s history. By turning points, we mean the occasions when the soil either stopped wetting and began to dry, or vise-versa (i.e., there was a change in sign in the time derivative of the moisture content). A notation for recording the turning points was presented by Enderby (1955), and adopted by Everett (1955) and Mualem (1974)

(73)  

Figure 2.19. Experimental primary drying scanning curves (solid lines) and theoretical points obtained using the independent domain model (small circles). Also included is the joint density function for this medium. Reproduced with permission from Poulovassilis (1962).

In this notation, subscripts denote the order of pressures, and relative position indicates whether the transition was wetting or drying. In the case shown, the medium wetted from h₀ to h₁, dried to h₂, and then re-wetted to the present pressure h.

In 73 we have written the relationship between pressure and moisture content using curly brackets {}. Later in this text, and in general in the literature, this relationship is written as a function with regular parentheses (), but this really is inappropriate, as θ is not a function of h but is a functional of h. That is to say, there is not a one-to-one mapping between θ and h without consideration of the antecedent conditions. Of course we can relate θ and h from a known initial state and through a known sequence of either θ or h as stated in 73, without requiring any knowledge of the other independent variables (most notably time). The ability of this model to reproduce hysteretic scanning curves is demonstrated in Figure 2.19 which is reproduced from Poulovassilis (1962).

So we have set up the conceptual framework for this method of modeling hysteresis, but how do we obtain the joint density function practically? You can just carry out a terrific number of experiments where you map out the entire domain of possible filling and draining pressures to obtain f(h_e,h_f) by brute force. While this would appear to be impractical, with numerically controlled valves and scales, this is not infeasible using an automated pressure cell. However, in 1973 Mualem introduced a simplification of this model which greatly eased the task of determining the density function.

Mualem (1973) noted that the joint density function f(h_e,h_f) could be well approximated by the product of two univariate density functions

(74)  

where g() and l() are probability density functions that depend only on h_e and h_f, respectively. This constrains, for instance, the filling pressure distributions to be the same, up to a constant multiplier, along any draining pressure line, as illustrated in Figure 2.20. Using Mualem’s formulation, we only need the main filling and emptying curves to obtain g() and l(), and from these results all the hysteretic results can be derived. Mualem showed that this formulation maintained its ability to reproduce scanning curves (Figure 2.21).

Parlange (1976) went on to show that an even more constraining similarity model based on data from the main draining curve alone is sufficient to reproduce the full family of scanning curves. This was demonstrated to be capable of fitting experimental data by Hogarth et al. in 1988, as shown in Figure 2.22, reproduced from their work.

We have gone into a fair bit of detail in how to handle hysterisis, but have not mentioned when and why it must be included when describing vadose zone processes. First off, to understand the general character of many processes, dealing specifically with hysterisis is not critical. For instance in irrigated agriculture, as long as you know the moisture content at field capacity and the permanent wilting point, you may not worry about the exact pressure the soil achieves at intermediate water contents. On the other hand, hysterisis has been shown to be central to the nature of fingered flow processes (see Chapter 3), as well as the movement of NAPLs (e.g., Lenhard and Parker, 1987; Lenhard et al., 1991). Beyond these very specific settings the literature does not provide general guidance on when and where hysterisis is most important. Beyond the examples given, we really don’t know when hysterisis will play a central role.

If hysterisis has been shown to be important to several important processes, and in general we can’t be sure hysterisis will not be influential, why haven’t more people included consideration of hysterisis in the description of vadose zone processes? We think this is simply a matter of practical constraints. None of the analytical solutions to flow, and precious few of the numerical codes for modeling flow in the vadose zone consider the possibility of hysterisis. Until we have convenient numerical codes which include hysteretic effects, the possible effects of hysterisis will continue to be ignored in quantitative analyses. Of course nature doesn’t care if hysterisis is hard to model or easy to model: it happens! Failing access to a useful model which includes hysterisis, a responsible approach is to qualitatively assess the potential impact of hysterisis on a problem. Given the rapid increase in availability and sophistication of numerical models, we believe that hysterisis will be included in many popular models in the coming decade, and is in fact available in research models at this time (White et al., 1995).

Geometric Relations of Mualem Similarity:

Emptying pressure: and for any h

Filling Pressure: and for any h

Figure 2.21. Mualem’s demonstration of the ability of the independent domain model with the similarity assumption to reproduce experimental scanning curves found by Poulovassilis (1970) (Reproduced from Mualem, 1974, with permission).

Figure 2.22. Primary draining curves as found experimentally and reproduced using the Parlange similarity model form of the independent domain model (Reproduced from Hogarth et al., 1988).

[Insert “Measurement of water potential” here]

Measurement of Water Content

The ability to accurately measure water content in time and space is central to the study of vadose zone hydrology. Remarkably, the ability to measure this critical parameter continues to present many challenges, with significant developments continuing to come forward.

The benchmark for water content measurements is the direct determination by oven drying (Gardner, 1986). In this method a sample of soil is gathered from the field and weighed. It is then placed in a 100 – 110℃ oven until it’s weight has stabilized. The mass basis water content may then be computed using 8. Though simple to describe, even this method is subject to significant error and interpretation. For instance, if stones are avoided, either as an unplanned feature of the method of collection or by design, the value of water content obtained will be systematically higher than the true field value. Furthermore, much of the water bound to clay remains in place at these temperatures. Thus you are not obtaining the total water content, but that which can be liberated under these specific conditions. Further discussion of these issues is presented by Gardner (1986). Perhaps the most problematic aspects of this method are that it is destructive and cannot be completed in the field in a timely manner which renders the method inappropriate for many applications.

There are a variety of field methods, all of which use indirect measurements of bulk soil physical properties to infer a moisture content. We will discuss three methods that illustrate the key physical properties that can be related to water content. These are the neutron thermalization method (neutron probe), the gamma ray attenuation method (gamma probe) and the dielectric method (time domain reflectometry).

Neutron Probe

The neutron probe has long been the most accurate method of measuring water content in the field. The method relies on the thermalization, or slowing, of neutrons. Relatively high energy neutrons are directed into the medium from a radioactive source, with a detector provided to observe the resulting number of slow neutrons that are back-scattered. These neutrons then collide with the various atoms in the soil. The loss of energy which occurs upon each collision is far greater when the collision is with an object of mass similar to the neutron. For instance, in a single collision the loss of energy from a neutron to hydrogen is about 4,000 time that of a collision with oxygen (Weinberg and Wigner, 1958). So the presence of hydrogen is closely correlated with the production of thermalized neutrons.

Since water is the largely hydrogen, the rate of thermalized neutron production can be calibrated to the surrounding volumetric water content. Site calibration is required to account for the local effects of the mineral and organic components and structurally bound water in the soil which will also contribute to neutron thermalization. The calibration equation is typically written in the form

(75)  

where θ is the volumetric water content, I_o is the incident neutron flux, I is the thermalized neutron count, and a and b are calibration parameters that depend upon the neutron source strength, soil properties, and detector characteristics. Since the soil characteristics will vary significantly in the features which will effect this calibration (e.g., organic matter), proper calibration must be carried out for each measurement location. This presents some logistical difficulties. First, the measurement must be taken at least two moisture contents to be able to obtain the two parameters. Therefore one must either come to the field twice, under differing moisture regimes, or one must carry out the calibration in the lab. Since the calibration requires a volume of soil at least as large as the instrument’s sphere of influence (about 15 cm in wet soil and 75 cm in dry soil), and in light of the sensitivity to bulk density, it is impractical to obtain sufficient samples to calibrate the probe using disturbed laboratory samples. This leaves only in situ calibration as being a viable option. This calibration is typically carried using the samples obtained from the original probe installation, then returning to the field when a second measurement can be made when the soil has quite different moisture content. In fact, while two points are sufficient in theory, a reliable calibration generally requires many more measurements, typically more like ten samples.

A great deal of work has been carried out to obtain estimates of the error associate with neutron probe measurements, calibration of measurements made close to the soil surface, and other logistical issues. A good starting place to learn about these practical concerns is Gardner (1986). Use of the neutron probe method has diminished in the past decade due to the logistical issues involved with handling and transporting radioactive sources. Yet due to the precision which may be obtained, the large sphere of influence of measurements, and low cost of operation, this method remains a mainstay of vadose moisture measurement.

Gamma Attenuation

The gamma attenuation method of liquid content measurement relies on the adsorption of gamma rays. The measurement is made by locating the sample under study between a gamma source and detector. The attenuation of the beam follows Beer’s law

(76)  

where the subscript λ denotes the gamma frequency of interest, I_oλ is the detected gamma intensity without sample, I_λ is the detected gamma intensity after passing through the sample, α_iλ is the attenuation coefficient for material i, and x_i is the path length of gamma travel in material i. If all of the features of the sample are constant except for liquid content, the attenuation of non liquid components can be identified with a single multiplicative constant C. For the fluid content, x_i = Tθ_i where T is the beam path length in the sample.

Each gamma energy will be attenuated at differing rates for a given liquid. By using multiple frequencies of gamma radiation it is possible to track the content of multiple variable constituents. For this reason, the dual frequency gamma method is the most common method of documenting the joint movement of water/oil systems (e.g., Schiegg, 1990), and may be used to document swelling effects (Gardner et al., 1972).

Due to the need to place the sample between source and detector, the gamma method is not often convenient to apply in the field, but is widely used in laboratory environments. Furthermore, given the use of intense radioactive sources, the gamma method is limited to highly controlled settings. For an excellent overview of this method, see Gardner (1986).

Time Domain Reflectometry

The measurement of water content using time domain reflectometry (TDR), although long standing (e.g., Fellner-Feldegg, 1969), has rapidly grown in popularity with the advent of affordable TDR systems. The TDR method can be configured with multiplexing to allow nearly continuous automated measurement at multiple locations with a precision sufficient for most applications. Furthermore, TDR requires less calibration than the neutron probe method, and eliminates the use of radioactive sources. Finally, TDR can be used to measure bulk electrical conductivity in parallel with water content, which can be very useful in tracer and salinity studies. Disadvantages of the method are a relatively small sphere of influence (<10 cm radius), and the high cost of probes (on the order of $100 per measurement location). Due to the significant benefits of TDR and developmental nature of this technology, we will provide a slightly more in-depth discussion of this method.

TDR depends upon the measurement of travel time of an electromagnetic wave down a wave guide which is either inserted or laid on the soil surface (Selker et al., 1993; Maheshwarla et al., 1995). Travel time increases with the square root of increase in the real component of the dielectric surrounding the probe. Since the mineral portion of porous media typically have constant dielectric (between 2-5) at frequencies > 50 MHz (Hoekstra and Delaney, 1974; Campbell, 1990), and liquid water in such media has a dielectric of 70-80 (temperature dependent), variation in signal travel time along the probes is attributed to variations in water content. Typical accuracy of TDR measurements better than +/-2% volumetric (Roth et al., 1990).

With the travel times of signals along TDR probes being in the order of nanoseconds, the precise measurement of this travel time has only been economically feasible in the last 15 years with the advent of inexpensive dedicated TDR oscilloscopes. The specific implementation of the TDR method has been intimately tied to the characteristics of the measurement tools. In particular, the Tektronix 1502 series TDR meters which were designed for field evaluation of signal line integrity have been the primary instrument employed in TDR moisture measurement. A key determinant of spatial resolution in TDR is the rise time of the device; approximately 2×10^-10 s for the Tektronix 1502 series. This dictates the minimum probe length that can be used. The change in travel time along a probe, D_t, for a change in travel velocity from v₁ to v₂ is given by

(77)  

where L is the probe length, with the factor 2 coming from the fact that two-way travel time is observed. Solving for L, we may calculate the minimum probe length for a given ability to measure D_t, and a given tolerance in velocity change. For example, when considering a dry soil-water system with a dielectric of 10, one can ask what minimum probe length allows determination of a change 2 percent in volumetric moisture content. This would create a change in dielectric of about 1.5, with v = 3×10⁸ m/s x K^1/2, the minimum desirable probe length is approximately 12 cm. Measurement accuracy can be increased, or the probe length decreased, when using a faster TDR system. Kelly et al. (1995) determined that probes as short as 0.02 m could be used with a TDR system which had a 25 ps rise time. The physical limitations of the frequency dependence of the dielectric of water become significant at these time scales. The period of molecular resonance of water is of the order of 20 ps. The real dielectric of water drops dramatically for signals with periods in this range, while the complex component (which is an energy loss term) climbs (Merabet and Bose, 1988). The combined loss of propagation delay and energy for single components with periods < 40 ps obviates any advantage of operating at frequencies in excess of about 10 Ghz. In addition to considerations of length, issues of spacing, conductor diameter, and number of wave guides are all significant in determining the sampling volume and sensitivity of the device (Knight, 1992; Whalley, 1993).

An additional limitation experienced in the application of TDR has been signal loss due to bulk conductivity of the porous media, which essentially shorts out the signal sent down the TDR probe (Arcone, 1986). The loss of signal strength due to the electrical conductivity of the media can be exploited to provide estimates of soil solution salinity, as demonstrated by Dalton et al. (1984) and many others. Signal loss difficulties are also encountered using Ground Penetrating Radar (GPR) methods in fine-textured media, as will be described later in this chapter. The signal loss problem can in part be mitigated through the use of minimum length insulated probes (e.g., Kelly et al., 1995), but it requires individual calibration curves for each probe geometry. In addition, high conductivity conditions reduces signal energy due to eddy currents established within the media, thereby limiting the range of applicability of the insulated probes to solutions with concentration of smaller than one molar total salinity, about 25 times greater than can be accommodated by un-insulated probes.

TDR methods are well established in vadose zone monitoring, yet there are various areas where significant enhancements can be achieved. One recent advancement is the development of multilevel shorting-diode TDR probe (Cook et al., 1992), which allows the monitoring of a linear moisture profile with the installation of a single probe. Another opportunity for enhancement of TDR methods lies in the interpretation of wave forms. At present, most TDR signals are interpreted with an empirical relationship between moisture content and travel time as developed by Topp et al. (1980). More recently, authors have derived physically-based methods to interpret the signals (Herkelrath et al., 1991; Heimovaara, 1993; Whalley, 1993), thereby allowing a more clear distinction between dielectric effects on the signal by the solid and liquid phase of the soil matrix. At this time, explicit electromagnetic models for signal propagation are being developed which predict the exact shape of the reflected wave form. These models can be used in a forward modeling mode, but the ultimate objective is to estimate phase contributions by inverse modeling using the full TDR trace as the basis of estimation of moisture content and soil conductivity, which will make more complete use of the wave form obtained from a TDR system.

Measurement of moisture content using TDR will continue to grow in the research community, but rapid adoption in the commercial sector has been hampered by the relatively high cost of TDR systems and the requirement of careful signal interpretation. New products have reduced these concerns, with increasingly robust signal interpretation software and more affordable meters. The TDR method has yet to be packaged into a system which is broadly accessible to individuals in the commercial sector. Nevertheless, TDR is continuing to grow in popularity in the measurement of moisture content within the vadose zone. Competing methods of moisture content measurement through sensing of the soil dielectric are under development or marketed (e.g., Campbell Scientific). As shown by Hoekstra and Delaney (1974) and Campbell (1990), these measurements must be made at frequencies of greater than 50 MHz to obtain a calibration which is independent of soil type.

Measurement of Hydraulic Conductivity

Hydraulic conductivity is one of the most important parameters in characterizing the transport of fluids and solutes through the vadose zone, but presents significant practical difficulty in measurement. In vadose applications, knowledge of the saturated conductivity is often of limited use, as these conditions are rare in most locations. Yet reliable methods of measuring unsaturated conductivity are very few. One of the most common approaches to obtaining the unsaturated conductivity relationship is to measure both the saturated conductivity and the water retention curve, and estimate the unsaturated conductivity based on a conceptual model relating permeability to water content (van Genuchten, 1980). This process is facilitated through application of computer programs written expressly for this purpose, such as RETC (van Genuchten et al., 1991).

The only direct method of measuring unsaturated conductivity in the field is through a pressure/moisture content profile method. The most widely used approach is the “instantaneous profile method” (IPM). The IPM requires that an area of land be subjected to ponded infiltration until the entire profile of interest is saturated. The site is instrumented with devices to measure both the moisture content and pressure with depth and time. At t = 0 the ponding is ceased, and the site is covered with an impermeable barrier to avoid evaporative losses. Pressure and water content measurements are made in time, with measurements extending to periods of weeks and months if low permeability values are required. Using the water content measurements one can calculate the flux at various depths by way of the total time rate of change in water content in the profile above the height of interest. By observing the local pressure gradients, as measured with tensiometers, one then has values for q and dh/dz in Darcy’s law, and it is straightforward to solve for K. A very nice added benefit of this method is that the joint pressure/water content measurements provide an in situ determination of the water retention relationship. Although the IPM is a direct measurement, due to the large number of measurements that need to be made over long periods of time, it is difficult to carry out, and often provide results which are less than satisfactory (e.g., Selker, 1994).

The remaining methods of measurement rely on an analytical solution for the flow from some simple geometric boundary condition. The solution must predict the time rate of infiltration in terms of the desired parameters. The non-linearity in Richards equation generally leads to the adoption of a solution with very restrictive assumptions in order to allow development of a general solution. The two most widely used assumptions are (1) that the infiltration process may be approximated by a sharp wetting front followed by a region of nearly saturated media (Green and Ampt infiltration: see Chapter 3) and (2) that the hydraulic conductivity function may be approximated by the Gardner (1958) relationship

(78)  

where α is a parameter approximately equal to -1/h_we, and thus is large for coarse textured soils (on the order of 10 m^-1 for sands) and small for fine textured soils (on the order of 0.1 m^-1 for clays).

The best control over the boundary at which infiltration occurs is obtained by carrying out surface measurements. Ring infiltration methods cover a broad set of these tests. In general these method involve establishing a prescribed pressure over a circular area and observing the resulting infiltration. The scale of these measurements ranges from 100 cm² (Ankeny et al., 1988) to 4 m² (ASTM, 1990). An approach that make use of the Green and Ampt concepts was introduced by Bouwer, and is discussed in section 3.1.1.3. Here we will introduce the ponded double ring method, and some newer tension infiltrometer approaches.

By far and a way the most common hydraulic test carried out for to estimate soil hydraulic properties is the double ring infiltrometer (Figure 2.40). The rings are inserted deep enough to preclude leakage from the outer ring, and to have the tops of the rings level with each other. A constant water level is quickly established in both rings to the same level, and the infiltration of water may be measured by watching the drop in water level on a floating ruler, or by using a Marriott bottle constant head water source.

The data may be analyzed using any of several infiltration models. Below we list the Horton Equation, the two-term Philip model, the Green and Ampt model and the Brutsaert model.

Horton Equation (1940): In this empirical infiltration model the rate of infiltration, i, is given by

(79)  

where i_f is the infiltration rate after long time, i_o is the initial infiltration rate and β is and empirical soil parameter. Integrating this with time yields the cumulative infiltration

(80)  

The model is widely used, but since the parameters obtained are not directly associated with physically measurable soil properties (e.g., conductivity), it is of limited use except for predicting infiltration.

The two term Philip model (1957) is a simplification of the infinite series Philip introduced for infiltration, and is widely used. The result given by

(81)  

or for cumulative infiltration

(82)  

has the advantage of being closely related to measurable soil physical parameters (i.e., S is the soil sorptivity and A is typically taken to be an estimate of field saturated conductivity).

The Green and Ampt Model (1911) for constant head infiltration is derived in Chapter 3, and is also useful in interpreting infiltration results. The Green and Ampt result for constant head infiltration is given by

(83)  

where L is the depth of wetting front penetration, n is the porosity of the soil, d is the depth of ponding, and h_f is the water entry pressure for the soil. The cumulative infiltration is simply I = nL. To use this equation you must find the values of K_sat and h_f which give the best fit to the data.

The Brutsaert Model. Brutsaert (1977) revisited the infinite series solution of Philip (1957) and developed a very convenient result that provides an excellent approximation to the full series solution for infiltration, and is written explicitly in the physical variables of interest. For the rate of infiltration Brutsaert (1977) found

(84)  

where S is the sorptivity and 1/3<β<1 is a soil parameter related to the distribution of pore sizes. For soils with very wide pore size distributions you should take β close to 1, for most other soils β = 2/3 is quite good, while for uniform pore size soils, β = 1/3 is suggested. The cumulative infiltration using the Brutsaert model is simply

(85)  

Fitting this equation to infiltration data allows ready determination of both K_sat and S.

At negative pressures, features of textural interfaces and fissures cause drastic reduction in unsaturated permeability, while at near-zero pressure macroporosity can greatly increase permeability. Thus, the ability to measure hydraulic properties of vadose materials in-situ at realistic negative pressures has been recognized as an area in need of attention.

Controlled infiltration under tension to determine unsaturated water retention and transport properties has been considered for decades, however, only in the past few years have experimental and analytical tools been improved to make the technique practical in the field. Clothier and White (1981) introduced the “sorptivity tube” which allowed for fixed pressure (<0) application of water to a soil using a Marriott bottle supply with a sintered glass porous interface. This design was tested and improved, culminating in the fabric interface device presented by Perroux and White (1988). This device provided ready adjustment of interface pressure, could be made with a wide range of contact areas, and was robust and simple enough to be used conveniently in the field. Developed in parallel, Ankeny et al. (1988) fitted the tension infiltrometer with two pressure transducers which allowed automated acquisition of both the pressure at the soil interface and the flux from the device. Further, this device was fitted with three adjustable bubble tubes that allowed for rapid evaluation of flux from a series of interface pressures. A key development in interpretation of the results of these devices was identification of the utility of the result of Wooding (1968) to disk infiltrometer (Ankeny et al., 1991). Wooding found that using a Gardner exponential conductivity function 78 the flux into a circular wetted area of radius r, at steady state, was

(86)  

where Q (measured) is the flux into the soil over a device of radius r (measured), h (measured) is the pressure at which the infiltration boundary, K_s (unknown) is the saturated hydraulic conductivity, and α (unknown) is as defined in ??. By carrying out infiltration tests at two or more pressures, or with two or more radius disks, one may solve for both of the unknowns in 86.

One significant shortcoming of Wooding’s approach is the notion of a well-defined steady-state condition, and the time required to achieve this condition. As shown by Smettem et al. (1995), these effects may take hours to become small, and the outflow might appear to be steady far before the true steady state has been achieved. Haverkamp et al. (1994) present an approximate solution for the time-dependent rate of infiltration for one and three dimensional cases as

(87)  

(88)  

where S_o is the sorptivity, β is a shape parameter (0<β<1), γ is a parameter with value between 0.6 and 0.8 which Smettem et al. (1994) found to be well approximated by 0.75, K_o & K_n, and θ_o & θ_n are the conductivity and moisture content at the antecedent and instrument pressures respectively. As shown by Smettem et al. (1995), a double ring tension infiltrometer may be made in which the outflow of the inner ring may be assumed to be one-dimensional, and that of the combined inner and outer rings three-dimensional.

Smettem et al. (1995) show that a good estimate of sorptivity may be obtained using 88 by plotting I_3d/t^1/2 against time. We see that 88 then reduces to

(89)  

for some constant A. By fitting a straight line to 89, the t=0 value of this relationship provides an estimate of S_o. Smettem et al. found this value to be consistent with values obtained from one-dimensional infiltration tests. This result may be used in conjunction with 86 or 88 for K_s.

These measurements are typically made at a series of tensions and often with two sizes of circular sources. The measurements require access to a horizontal soil surface, thereby limiting application of this technique to the first few meters of the vadose zone, although a down-hole tension infiltrometer designs has also been developed (Sisson and Honeycutt, 1994).

An alternative method of interpreting tension infiltrometer results is to carry out the estimation of sorptivity and hydraulic conductivity functions by inverting a numerical simulation of the process (Simunek et al., 1995 and 1996). Since this method makes full use of the early time data, it holds the promise of significantly increasing the utilization of site-specific hydraulic properties by reducing the experimental time required to obtain reliable values. The accuracy and practicality of these methods has also been explored recently, comparing a number of methods both for field and laboratory application (e.g., Logsdon and Jaynes, 1993; Gupta et al., 1993; Stolte et al., 1994; Cook and Broeren, 1994). It is clear that each of the present implementations of the tension infiltrometer method have difficulty under some situations, but the outlook for improved methods of measurement appears quite bright in the coming years using this technique.

To determine the hydraulic properties of vadose materials at depth, down-hole methods are essential. This is an area of active development, with some very useful techniques available. The method most widely employed is to auger a hole of constant radius to the depth at which the measurement is required. A constant head of water is then established in the hole, while the rate of infiltration is observed. Several tools have been developed for this test, and generally consist of a Marriott bottle source of water with an extension tube to reach the required depth. Such tools can easily be made without access to specialized equipment, as shown in Figure 2.41.

The analysis of constant head well permeameter can be carried out using a number of formulations. The original work directed to analysis of such data was presented by Glover (Zangar, 1953). As presented by Elrick and Reynolds (1993), Glover found that at steady state the rate of infiltration from such a permeameter could be approximated by

(90)  

where K_fs is the “field saturated” hydraulic conductivity and C is a dimensionless shape coefficient reflecting the geometry of the wetted area. This was estimated by Glover to be

(91)  

This solution assumes complete saturation in the wetted region, and does not include gravity driven flow (Elrick and Reynolds, 1993). For ratios of H/a < 2, this solution has been shown to be reasonably accurate, while for H/a > 10, it under estimates the correct value of C by about 50%, leading to an error in K_fs of about 40%. To account for the effects of gravity and unsaturated flow effects, Elrick et al. (1989) provided the solution

(92)  

Elrick and Reynolds (1992a and 1992b) suggest that α* may be estimated based on field conditions using the values given in Table 2.5. Alternately, tests using multiple depths of ponding may be carried out, providing for multiple equations from which the values of K_fs and α* may be estimated using regression analysis. Difficulties with this approach have been reported (e.g. Elrick and Reynolds, 1992b) due to heterogeneity in soils which will yield inconsistent values of the parameters depending upon the depth of ponding. This can lead to non-physical values of estimated parameters (e.g., negative conductivity), and thus should only be employed in soils without significant vertical variability in properties.

Table 2.5 Suggested values for α* as provided by Elrick and Reynolds (1992).

Index

A
admixture

anticlastic

atmospheric pumping

B

Beer’s law

Brutsaert Infiltration model

Buchner funnel

Buckingham

Bukner funnel

bulk density

C

calcite

capillary tube

characteristic curves

Characteristic Curves

chemical potential

conductivity, down hole measurement

constant head permeameter

constituents

contact angle

conversions

creeping flow

curved fluid surface

D

Darcian Length

Darcian velocity

Darcy’s law

datum

dead end pore

density

dew point

dielectric

diffusion

diffusivity

dipole

DNAPL

double ring infiltrometer

down-hole tension infiltrometer

Droplet

dynamic viscosity

E

entropy

F

feldspar

field capacity

field saturated

Fokker-Planck equation

free product

functional

fundamental differential form

G

gamma attenuation

gamma method, dual frequency

Gardner conductivity

gas movement

gas permeability

gaseous

geometric similarity

Gibbs-Duhem relationship

GPR

gravitational potential

Green and Ampt

ground penetrating radar

gypsum block

H

Hagen Poiseuille

Hagen-Poiseuille Equation

Haines jump

Henry’s law

Horton Equation

hydraulic conductivity, measurement

Hydrostatics, 21

hysteresis, gypsum blocks

hysteretic

hysterisis

I

immiscible

independent domain assumptions

inertial forces

infiltration models

ink bottle effect

instantaneous profile method

internal energy

intrinsic permeability

irreducible

K

kinematic viscosity

Klinkenberg effect

L

Liquid Content

LNAPL

LNAPLs in wells

M

macropore

matric potential

mean free path length

microbe

Miller similarity

N

NAPL

NAPLs

Navier Stokes Equation

negative pressure

neutron probe

neutron probe, errors

Newtonian fluid

non aqueous phase liquids

non-spreading NAPL

non-wetting fluid

O

oil spill

Osmotic

osmotic potential

oven drying

P

particle size distribution

pendular

permanent wilting point

phases

piezometer

pore size distribution

porosity

potential energy

pressure head

pressure potential

pressure potential, measurement

psychrometer

psychrometric equation

pyrite

R

rain-drop effect

reduced radius

relative humidity

Representative Elementary Volume

residual

REV

Reynolds number

Richards equation

S

satiated

satuated

saturated liquid content

scaled conductivity

scaled pressure

scaled sorptivity

scaled time

scanning curve

similar media

similarity

similarity assumptions

soil diffusivity

sorptivity tube

specific flux

specific gravity

spreading NAPL

suction

surface tension

Surface Tension

surface tensions

T

TDR

TDR, insulated probes

TDR, shorting diode

tensiometer

tension

tension infiltrometer

tension infiltrometer, numerical simulation

thermalization

thermocouple

Thermodynamics

time domain reflectometery

Time domain reflectometery

tortuosity

Total Bulk Density

total head

turbulence

two-term Philip model

U

units

V

Van der Waals

vapor extraction

vapor pressure

Vapor Pressure

viscous dissipation

viscous forces

void fraction

Void Ratio

volatile organics

W

water content, measurement

wetability

Wooding solution

Y

Youngs-Laplace equation

References

Brown and Montgomery, 1913, “Dehydration of Clays,” Bureau of Standards Technologic Paper no. 21, Washington, D.C.: https://nvlpubs.nist.gov/nistpubs/nbstechnologic/nbstechnologicpaperT21.pdf

Kirkham and Powers (1972)

Revil, A. and A. Jardani, The Self potential Method: theory and applications in environmental geosciences, Cambridge University Press, ISBN 978-1-107-01927-0, 2013.

Sposito (1989)

†† A few words about the terms in the K tensor.  So the idea here is to be able to represent the fact that the resistance to flow might be different per direction due to orientation of the particles at a scale below that of the REV.  So we will “take apart” the conductivity to explicitly represent the relationship between the gradient in potential in a particular direction, and the resulting flow that might happen in that same direction, or along a different axis.  Writing Darcy’s law

To get a feel for the terms in the K-tensor, look at the terms in the x-coordinate of flux.  We see that gradients in each of the x, y, and z directions can give rise to flux in the x direction.  In general we see that the term K_ab is the proportionality constant relating gradients in potential in the b-direction to flux that results in the a-direction.  How could such a thing happen?  If the medium is aniostropic (i.e. permeability is a function of direction) then gradients that are not aligned with the directionality of permeability will give rise to flux not parallel with potential gradients.  An easy analogy is movement on a slide in a playground.  The direction of energy gradient due to gravity is vertically downward, but your resulting motion includes a lateral component, along the direction of the slide.  This really is similar to the motion of water in a clay with oriented particles:  the water will largely travel parallel to the sheet-structures of the clay rather than across the clay sheets because the resistance is so much lower, even if the gradient is not exactly in the same direction as the orientation of the clay sheets. If the coordinate system is aligned such that the coordinate system and primary directions of anisotropy are coincident the off diagonal terms in the permeability tensor will be zero (i.e., the K_ab terms where a ≠ b), leaving permeabilities in the three coordinate directions.  If these are all equal, then the tensor collapses to a scalar.  So the only reason to take the permeability as a tensor is to capture the effects of anisotropy.

† Some of you may ask what is this up-side-down triangle all about?  The ∇ symbol stands for the “dell” operator.  It is just mathematical short hand for calculating the three-directional derivative instead of the one directional derivative.  It is defined as

where i, j, and k are the unit vectors in the x, y, and z directions.  The result of “operating” on a function with ∇ is to obtain a new vector which is the slope of the function.   In practice, for any scalar valued function of x,y,z, ∇F tells us the magnitude and direction of the slope of F: ∇F points directly towards the steepest direction of uphill with a length proportional to the slope of the hill.   Later you will also see this symbol in the notation ∇·F.  The dot just tells us to take the dell and calculate the dot product of that and the function F (which needs to be a vector for this to make sense).  This term “dell-dot-F” then tells us the divergence of the F function.  If F were three directional flux content, this would be the amount of water leaving the point (x,y,z).  Notice that this is a scalar result!  So ∇F needs a scalar function F and gives a vector slope, while ∇·F uses a vector function F and gives a scalar result.