Physical Properties and Principles

2.1 Darcy’s Law

The birth of groundwater hydrology as a quantitative science can be traced to the year 1856. It was that year that French hydraulic engineer named Henry Darcy published his report on the water supply of the city of Dijon, France. In the report Darcy described a laboratory experiment that he had carried out to analyze the flow of water through sands. The result of his experiment can be generalized into the empirical law that now bears his name.

Consider an experimental apparatus like that shown in Figure 2.1. A circular cylinder of cross section A is filled with sand, stoppered at each end, and outfitted with inflow and outflow tubes and a pair of manometers. Water is introduced into the cylinder and allowed to flow through it until such time as all the pores are filled with water and the inflow rate Q is equal to the outflow rate. If we set an arbitrary datum at elevation z = 0, the elevations of the manometer intakes are z₁ and z₂ and the elevations of the fluid levels are h₁ and h₂. The distance between the manometer intakes is Δl.

**Figure 2.1** Experimental apparatus for the illustration of Darcy’s law.

We will define v, the specific discharge through the cylinder, as

$v = \frac{Q}{A}$ (2.1)

If the dimensions of the Q are [L³/T] and those of A are [L²], v has the dimensions of velocity [L/T].

The experiments carried out by Darcy showed that v is directly proportional to h₁ – h₂ when Δl is held constant, and inversely proportional to Δl when h₁ – h₂ is held constant. If we define Δh = h₂ – h₁ (an arbitrary sign convention that will stand us in good stead in later developments), we have $v \propto \Delta h$ and $v \propto 1/\Delta l$ . Darcy’s law can now be written as

$v = -K\frac{\Delta h}{\Delta l}$ (2.2)

Or in differential form,

$v = -K\frac{dh}{dl}$ (2.3)

In Eq. (2.3), h is called the hydraulic head and dh/dl is the hydraulic gradient. K is a constant of proportionality. It must be a property of the soil in the cylinder, for were we to hold the hydraulic gradient constant, the specific discharge would surely be larger for some soils than for others. In other words, if dh/dl is held constant, v $\propto$ K. The parameter K is known as the hydraulic conductivity. It has high values for sand and gravel and low values for clay and most rocks. Since Δh and Δl both have units of length [L], a quick dimensional analysis of Eq. (2.2) shows that K has the dimensions of a velocity [L/T]. In Section 2.3, we will show that K is a function not only of the media, but also of the fluid flowing through it.

An alternative form of Darcy’s law can be obtained by substituting Eq. (2.1) in Eq. (2.3) to yield

$Q = -K\frac{dh}{dl}A$ (2.4)

This is sometimes compacted even further into the form

$Q = KiA$ (2.5)

where i is the hydraulic gradient.

Darcy’s law is valid for groundwater flow in any direction in space. With regard to Figure 2.1 and Eq. (2.3), if the hydraulic gradient dh/dl and the hydraulic conductivity K are held constant, v is independent of the angle θ. This is true even for θ values greater than 90° when the flow is being forced up through the cylinder against gravity.

We have noted that the specific discharge v has the dimensions of a velocity, or flux. For this reason it is sometimes known as the Darcy velocity or Darcy flux. The specific discharge is a macroscopic concept and it is easily measured. It must be clearly differentiated from the microscopic velocities associated with the actual paths of individual particles of water as they wind their way through the grains of sand (Figure 2.2). The microscopic velocities are real, but they are probably impossible to measure. In the remainder of the chapter we will work exclusively with concepts of flow on a macroscopic scale. Despite its dimensions we will not refer to v as a velocity; rather we will utilize the more correct term, specific discharge.

**Figure 2.2** Macroscopic and microscopic concepts of groundwater flow.

This last paragraph may appear innocuous, but it announces a decision of fundamental importance. When we decide to analyze groundwater flow with the Darcian approach, it means, in effect, that we are going to replace the actual ensemble of sand grains (or clay particles or rock fragments) that make up the porous medium by a representative continuum for which we can define macroscopic parameters, such as the hydraulic conductivity, and utilize macroscopic laws, such as Darcy’s law, to provide macroscopically averaged descriptions of the microscopic behavior. This is a conceptually simple and logical step to take, but it rests on some knotty theoretical foundations. Bear (1972), in his advanced text on porous-media flow, discusses these foundations in detail. In Section 2.12, we will further explore the interrelationships between the microscopic and macroscopic descriptions of groundwater flow.

Darcy’s law is an empirical law. It rests only on experimental evidence. Many attempts have been made to derive Darcy’s law from more fundamental physical laws, and Bear (1972) also reviews these studies in some detail. The most successful approaches attempt to apply the Navier-Stokes equations, which are widely known in the study of fluid mechanics, to the flow of water through the pore channels of idealized conceptual models of porous media. Hubbert (1956) and Irmay (1958) were apparently the earliest to attempt this exercise.

This text will provide ample evidence of the fundamental importance of Darcy’s law in the analysis of groundwater flow, but it is worth noting here that it is equally important in many other applications of porous-media flow. It describes the flow of soil moisture and is used by soil physicists, agricultural engineers, and soil mechanics specialists. It describes the flow of oil and gas in deep geological formations and is used by petroleum reservoir analysts. It is used in the design of filters by chemical engineers and in the design of porous ceramics by materials scientists. It has even been used by bioscientists to describe the flow of bodily fluids across porous membranes in the body.

Darcy’s law is a powerful empirical law and its components deserve our more careful attention. The next two sections provide a more detailed look at the physical significance of the hydraulic head h and the hydraulic conductivity K.

2.2 Hydraulic Head and Fluid Potential

The analysis of a physical process that involves flow usually requires the recognition of a potential gradient. For example, it is known that heat flows through solids from higher temperatures toward lower and that electrical current flows through electrical circuits from higher voltages toward lower. For these processes, the temperature and the voltage are potential quantities, and the rates of flow of heat and electricity are proportional to these potential gradients. Our task is to determine the potential gradient that controls the flow of water through porous media.

Fortunately, this question has been carefully considered by Hubbert in his classical treatise on groundwater flow (Hubbert, 1940). In the first part of this section we will review his concepts and derivations.

Hubbert’s Analysis of the Fluid Potential

Hubbert (1940) defines potential as “a physical quantity, capable of measurement at every point in a flow system, whose properties are such that flow always occurs from regions in which the quantity has higher values to those in which it has lower, regardless of the direction in space” (p. 794). In the Darcy experiment (Figure 2.1) the hydraulic head h, indicated by the water levels in the manometers, would appear to satisfy the definition, but as Hubbert points out, “to adopt it empirically without further investigation would be like reading the length of the mercury column of a thermometer without knowing that temperature was the physical quantity being indicated” (p. 795).

Two obvious possibilities for the potential quantity are elevation and fluid pressure. If the Darcy apparatus (Figure 2.1) were set up with the cylinder vertical (θ = 0), flow would certainly occur down through the cylinder (from high elevation to low) in response to gravity. On the other hand, if the cylinder were placed in a horizontal position (θ = 90°) so that gravity played no role, flow could presumably be induced by increasing the pressure at one end and decreasing it at the other. Individually, neither elevation nor pressure are adequate potentials, but we certainly have reason to expect them to be components of the total potential quantity.

It will come as no surprise to those who have been exposed to potential concepts in elementary physics or fluid mechanics that the best way to search out our quarry is to examine the energy relationships during the flow process. In fact, the classical definition of potential as it is usually presented by mathematicians and physicists is in terms of the work done during the flow process; and the work done in moving a unit mass of fluid between any two points in a flow system is a measure of the energy loss of the unit mass.

Fluid flow through porous media is a mechanical process. The fluid forces driving the fluid forward must overcome the frictional forces set up between the moving fluid and the grains of the porous medium. The flow is therefore accompanied by an irreversible transformation of mechanical energy to thermal energy through the mechanism of frictional resistance. The direction of flow in space must therefore be away from regions in which the mechanical energy per unit mass of fluid is higher and toward regions in which it is lower. In that the mechanical energy per unit mass at any point in a flow system can be defined as the work required to move a unit mass of fluid from an arbitrarily chosen standard state to the point in question, it is clear that we have uncovered a physical quantity that satisfies both Hubbert’s definition of a potential (in terms of the direction of flow) and the classical definition (in terms of the work done). The fluid potential for flow through porous media is therefore the mechanical energy per unit mass of fluid.

It now remains to relate this quantity to the elevation and pressure terms that we anticipated earlier. Consider an arbitrary standard state (Figure 2.3) at elevation z = 0 and pressure p = p₀ where p₀ is atmospheric. A unit mass of fluid of density ρ₀ will occupy a volume V₀, where V₀ = 1/ρ₀. We wish to calculate the work required to lift the unit mass of fluid from the standard state to some point P in the flow system which is at elevation z and where the fluid pressure is p. Here, a unit mass of the fluid may have density ρ and will occupy a volume V = 1/ρ. In addition, we will consider the fluid to have velocity v = 0 at the standard state and velocity v at the point P.

**Figure 2.3** Data for calculation of mechanical energy of unit mass of fluid.

There are three components to the work calculation. First, there is the work required to lift the mass from elevation z = 0 to elevation z:

$w_1 = mgz$ (2.6)

Second, there is the work required to accelerate the fluid from velocity v = 0 to velocity v:

$w_2 = \frac{mv^2}{2}$ (2.7)

Third, there is the work done on the fluid in raising the fluid pressure from p = p₀ to p:

$w_3 = m \int^p_{p_0}\frac{V}{m}dp = m \int^p_{p_0}\frac{dp}{\rho}$ (2.8)

If the fluid were to flow from point P to a point at the standard state, Eq. (2.6) represents the loss in potential energy, Eq. (2.7) is the loss in kinetic energy, and Eq. (2.8) is the loss in elastic energy, or p – V work.

The fluid potential Φ (the mechanical energy per unit mass) is the sum of w₁, w₂, and w₃. For a unit mass of fluid, m = 1 in Eqs. (2.6), (2.7), and (2.8), and we have

$\Phi = gz + \frac{v^2}{2}+\int^p_{p_0}\frac{dp}{\rho}$ (2.9)

For porous-media flow, velocities are extremely low, so the second term can almost always be neglected. For incompressible fluids (fluids with a constant density, so that ρ is not a function of p), Eq. (2.9) can be simplified further to give

$\Phi = gz + \frac{p-p_0}{\rho}$ (2.10)

Our earlier premonitions as to the likely components of the fluid potential are now seen to be correct. The first term of Eq. (2.10) involves the elevation z and the second term involves the fluid pressure p.

But how do these terms relate to the hydraulic head h? Let us return to the Darcy manometer (Figure 2.4). At P, the fluid pressure p is given by

$p = \rho g\psi + p_0$ (2.11)

where ψ is the height of the liquid column above P and p₀ is atmospheric pressure or pressure at the standard state. It is clear from Figure 2.4 and Eq. 2.11 that

$p = \rho g(h-z) + p_0$ (2.12)

**Figure 2.4** Hydraulic head h, pressure head ψ, and elevation head z for a laboratory manometer.

Substituting Eq. (2.12) in Eq. (2.10) yields

$\Phi = gz + \frac{[\rho g(h-z)+p_0]-p_0}{\rho}$ (2.13)

Or cancelling terms,

$\Phi = gh$ (2.14)

Our long exercise has led us to the simplest of conclusions. The fluid potential Φ at any point P in a porous medium is simply the hydraulic head at the point multiplied by the acceleration due to gravity. Since g is very nearly constant in the vicinity of the earth’s surface, Φ and h are almost perfectly correlated. To know one is to know the other. The hydraulic head h is therefore just as suitable a potential as is Φ. To recall Hubbert’s definition: it is a physical quantity, it is capable of measurement, and flow always occurs from regions where h has higher values to regions where it has lower. In fact reference to Eq. (2.14) shows that, if Φ is energy per unit mass, then h is energy per unit weight.

It is common in groundwater hydrology to set the atmospheric pressure p₀ equal to zero and work in gage pressures (i.e., pressures above atmospheric). In this case Eqs. (2.10) and (2.14) become

$\Phi = gz + \frac{p}{\rho} = gh$ (2.15)

Dividing through by g, we obtain

$h = z + \frac{p}{\rho g}$ (2.16)

Putting Eq. (2.11) in terms of gage pressures yields

$p = \rho g \psi$ (2.17)

and Eq. (2.16) becomes

$h = z + \psi$ (2.18)

The hydraulic head h is thus seen to be the sum of two components: the elevation of the point of measurement, or elevation head, z, and the pressure head ψ. This fundamental head relationship is basic to an understanding of groundwater flow. Figure 2.4 displays the relationship for the Darcy manometer, Figure 2.5 for a field measurement site.

Those who are familiar with elementary fluid mechanics may already have recognized Eq. (2.9) as the Bernoulli equation, the classical formulation of energy loss during fluid flow. Some authors (Todd, 1959; Domenico, 1972) use the Bernoulli equation as the starting point for their development of the concepts of fluid potential and hydraulic head.

**Figure 2.5** Hydraulic head h, pressure head ψ, and elevation head z for a field piezometer.

If we put Eq. (2.9) in terms of head and use a simplified notation, it becomes

$h_T = h+z + h_p + h_v$ (2.19)

where h_z is the elevation head, h_p the pressure head, and h_v the velocity head. In our earlier notation: h_z = z, h_p = ψ, and h_v = v²/2g. The term h_T is called the total head, and for the special case where h₀ = 0, it is equal to the hydraulic head h, and Eq. (2.18) holds.

Dimensions and Units

The dimensions of the head terms h, ψ, z are those of length [L]. They are usually expressed as “meters of water” or “feet of water.” The specification “of water” emphasizes that head measurements are dependent on fluid density through the relationship of Eq. (2.17). Given the same fluid pressure p at point P in Figure 2.5, the hydraulic head h and pressure head ψ would have different values if the fluid in the pores of the geological formation were oil rather than water. In this text, where we will always be dealing with water, the adjectival phrase will usually be dropped and we will record heads in meters.

As for the other terms introduced in this section; in the SI system, with its [M][L][T] base, pressure has dimensions [M/LT²] mass density has dimensions [M/L³], and the fluid potential, from its definition, is energy per unit mass with dimensions [L²/T²]. Table 2.1 clarifies the dimensions and common units for all the important parameters introduced thus far. Reference to Appendix I should resolve any confusion. In this text, we will use SI metric units as our basic system of units, but Table 2.1 includes the FPS equivalents. Table A1.3 in Appendix I provides conversion factors.

Note that Table 2.1 that the weight density of water, γ, defined by

$\gamma = \rho g$ (2.20)

is a more suitable parameter than the mass density ρ for the FPS system of units, which has force as one of its fundamental dimensions.

Table 2.1 Dimensions and Common Units for Basic Groundwater Parameters*

		Système International† SI		Foot-pound-second system,‡ FPS
Parameter	Symbol	Dimension	Units	Dimension	Units
Hydraulic head	h	[L]	m	[L]	ft
Pressure head	ψ	[L]	m	[L]	ft
Elevation head	z	[L]	m	[L]	ft
Fluid pressure	p	[M/LT²]	N/m² or Pa	[F/L²]	lb/ft²
Fluid potential	Φ	[L²/T²]	m²/s²	[L²/T²]	ft²/s²
Mass density	ρ	[M/L³]	kg/m³	–	–
Weight density	γ	–	–	[F/L³]	lb/ft³
Specific discharge	v	[L/T]	m/s	[L/T]	ft/s
Hydraulic conductivity	K	[L/T]	m/s	[L/T]	ft/s

*See also Tables A1.1, A1.2, and A1.3, Appendix I.
†Basic dimensions are length [L], mass [M], and time [T].
‡Basic dimensions are length [L], force [M], and time [T].

Piezometers and Piezometer Nests

The basic device for the measurement of hydraulic head is a tube or pipe in which the elevation of a water level can be determined. In the laboratory (Figure 2.4) the tube is a manometer; in the field (Figure 2.5) the pipe is called a piezometer. A piezometer must be sealed along its length. It must be open to water flow at the bottom and be open to the atmosphere at the top. The intake is usually a section of slotted pipe or a commercially available well point. In either case the intake must be designed to allow the inflow of water but not of the sand grains or clay particles that make up the geologic formation. It must be emphasized that the point of measurement in a piezometer is at its base, not at the level of the fluid surface. One can view the functioning of a piezometer much like that of a thermometer. It is simply the instrument, if such it can be called, used to determine the value of h at some point P in a groundwater reservoir. In recent years, the simple standpipe piezometer has been replaced in some applications by more complex designs utilizing pressure transducers, pneumatic devices, and electronic components.

Piezometers are usually installed in groups so that they can be used to determine directions of groundwater flow. In Figure 2.6(a) three piezometers tap a water-bearing geological formation. It is a worthwhile exercise to remove the instruments of measurement from the diagram [Figure 2.6(b)] and consider only the measured values. Flow is from higher h to lower, in this case from right to left.

**Figure 2.6** Determination of hydraulic gradients from piezometer installations.

If the distance between the piezometers were known, the hydraulic gradient dh/dl could be calculated; and if the hydraulic conductivity K of the geological formation were known, Darcy’s law could be used to calculate the specific discharge (or volume rate of flow through any cross-sectional area perpendicular to the flow direction).

Sometimes it is the vertical potential gradient that is of interest. In such cases a piezometer nest is utilized, with two or more piezometers installed side by side at the same location (or possibly in the same hole), each bottoming at a different depth and possibly in a different geological formation. Figure 2.6(c) and (d) shows a piezometer nest in a region of upward groundwater flow.

The distribution of hydraulic heads in a groundwater system is three-dimensional through space. The piezometer groupings shown in Figure 2.6 only prove the existence of components of flow in the directions indicated. If a large number of piezometers could be distributed throughout the three-dimensional hydrogeologic system, it would be possible to contour the positions of equal hydraulic head. In three dimensions the locus of such points forms an equipotential surface. In any two-dimensional cross section, be it horizontal, vertical or otherwise, the traces of the equipotential surfaces on the section are called equipotential lines. If the pattern of hydraulic heads is known in a cross section, flowlines can be constructed perpendicular to the equipotential lines (in the direction of the maximum potential gradient). The resulting set of intersecting equipotential lines and flowlines is known as a flow net. Chapter 5 will provide detailed instructions on the construction of flow nets, and Chapter 6 will prove their usefulness in the interpretation of regional groundwater flow.

Coupled Flow

There is now a large body of experimental and theoretical evidence to show that water can be induced to flow through porous media under the influence of gradients other than that of hydraulic head. For example, the presence of a temperature gradient can cause groundwater flow (as well as heat flow) even when hydraulic gradients do not exist (Gurr et al., 1952; Philip and de Vries, 1957). This component becomes important in the formation of frost wedges in soil (Hoekstra, 1966; Harlan, 1973).

An electrical gradient can create a flow of water from high voltage to low when earth currents are set up in a soil. The mechanism of flow involves an interaction between charged ions in the water and the electrical charge associated with clay minerals in the soil (Casagrande, 1952). The principle is used in soil mechanics in the electrokinetic approach to soil drainage (Terzaghi and Peck, 1967).

Chemical gradients can cause the flow of water (as well as the movement of chemical constituents through the water) from regions where water has higher salinity to regions where it has lower salinity, even in the absence of other gradients. The role of chemical gradients in the production of water How is relatively unimportant, but their direct influence on the movement of chemical constituents is of major importance in the analysis of groundwater contamination. These concepts will come to the fore in Chapters 3, 7, and 9.

If each of these gradients plays a role in producing flow, it follows that a more general flow law than Eq. (2.3) can be written in the form

$v = -L_1\frac{dh}{dl} - L_2\frac{dT}{dl} - L_3\frac{dc}{dl}$ (2.21)

where h is hydraulic head, T is temperature, and c is chemical concentration; L₁, L₂, L₃ are constants of proportionality. For the purposes of discussion, let us set dc/dl = 0. We are left with a situation where fluid flow is occurring in response to both a hydraulic head gradient and a temperature gradient:

$v = -L_1\frac{dh}{dl} - L_2\frac{dT}{dl}$ (2.22)

In general, L₁ dh/dl $\gg$ dT/dl.

If a temperature gradient can cause fluid flow as well as heat flow in a porous medium, it should come as no surprise to find that a hydraulic gradient can cause heat flow as well as fluid flow. This mutual interdependency is a reflection of the well-known thermodynamic concept of coupled flow. If we set dh/dl = i₁, and dT/dl = i₂ we can write a pair of equations patterned after Eq. (2.22):

$v_1 = -L_{11}i_1 - L_{12}i_2$ (2.23)

$v_2 = -L_{21}i_1 - L_{22}i_2$ (2.24)

where v₁ is the specific discharge of fluid through the medium and v₂ is the specific discharge of heat through the medium. The L’s are known as phenomenological coefficients. If L₁₂ = 0 in Eq. (2.23), we are left with Darcy’s law of groundwater flow and L₁₁ = 0 is the hydraulic conductivity. If L₂₁ = 0 in Eq. (2.24), we are left with Fourier’s law of heat flow and L₂₂ is the thermal conductivity.

It is possible to write a complete set of coupled equations. The set of equations would have the form of Eq. (2.23) but would involve all the gradients of Eq. (2.21) and perhaps others. The development of the theory of coupled flows in porous media was pioneered by Taylor and Cary (1964). Olsen (1969) has carried out significant experimental research. Bear (1972) provides a more detailed development of the concepts than can be attempted here. The thermodynamic description of the physics of porous media flow is conceptually powerful, but in practice there are very few data on the nature of the off-diagonal coefficients in the matrix of phenomenological coefficients L_ij. In this text we will assume that groundwater flow is fully described by Darcy’s law [Eq. (2.3)]; that the hydraulic head [Eq. (2.18)], with its elevation and pressure components, is a suitable representation of the total head; and that the hydraulic conductivity is the only important phenomenological coefficient in Eq. (2.21).

2.3 Hydraulic Conductivity and Permeability

As Hubbert (1956) has pointed out, the constant of proportionality in Darcy’s law, which has been christened the hydraulic conductivity, is a function not only of the porous medium but also of the fluid. Consider once again the experimental apparatus of Figure 2.1. If Δh and Δl are held constant for two runs using the same sand, but water is the fluid in the first run and molasses in the second, it would come as no surprise to find the specific discharge v much lower in the second run than in the first. In light of such an observation, it would be instructive to search for a parameter that can describe the conductive properties of a porous medium independently from the fluid flowing through it.

To this end experiments have been carried out with ideal porous media consisting of uniform glass beads of diameter d. When various fluids of density ρ and dynamic viscosity μ are run through the apparatus under a constant hydraulic gradient dh/dl, the following proportionality relationships are observed:

$v \propto d^2$

$v \propto \rho g$

$v \propto \frac{1}{\mu}$

Together with Darcy’s original observation that $v \propto - dh/dl$ , these three relationships lead to a new version of Darcy’s law:

$v = -\frac{Cd^2\rho g}{\mu}\frac{dh}{dl}$ (2.25)

The parameter C is yet another constant of proportionality. For real soils it must include the influence of other media properties that affect flow, apart from the mean grain diameter: for example, the distribution of grain sizes, the sphericity and roundness of the grains, and the nature of their packing.

Comparison of Eq. (2.25) with the original Darcy equation [Eq. (2.3)] shows that

$K = \frac{Cd^2\rho g}{\mu}$ (2.26)

In this equation, ρ and μ are functions of the fluid alone and Cd² is a function of the medium alone. If we define

$k = Cd^2$ (2.27)

then

$K = \frac{k\rho g}{\mu}$ (2.28)

The parameter k is known as the specific or intrinsic permeability. If K is always called hydraulic conductivity, it is safe to drop the adjectives and refer to k as simply the permeability. That is the convention that will be followed in this text, but it can lead to some confusion, especially when dealing with older texts and reports where the hydraulic conductivity K is sometimes called the coefficient of permeability.

Hubbert (1940) developed Eqs. (2.25) through (2.28) from fundamental principles by considering the relationships between driving and resisting forces on a microscopic scale during flow through porous media. The dimensional considerations inherent in his analysis provided us with the foresight to include the constant g in the proportionality relationship leading to Eq. (2.25). In this way C emerges as a dimensionless constant.

The permeability k is a function only of the medium and has dimensions [L²]. The term is widely used in the petroleum industry, where the existence of gas, oil, and water in multiphase flow systems makes the use of a fluid-free conductance parameter attractive. When measured in m² or cm², k is very small, so petroleum engineers have defined the darcy as a unit of permeability. If Eq. (2.28) is substituted in Eq. (2.3), Darcy’s law becomes

$v = \frac{-k\rho g}{\mu}\frac{dh}{dl}$ (2.29)

Referring to this equation, 1 darcy is defined as the permeability that will lead to a specific discharge of 1 cm/s for a fluid with a viscosity of 1 cp under a hydraulic gradient that makes the term ρg dh/dl equal to 1 atm/cm. One darcy is approximately equal to 10^-8 cm².

In the water well industry, the unit gal/day/ft² is widely used for hydraulic conductivity. lts relevance is clearest when Darcy’s law is couched in terms of Eq. (2.4):

$Q = -K\frac{dh}{dl}A$

The early definitions provided by the U.S. Geological Survey with regard to this unit differentiate between a laboratory coefficient and a field coefficient. However, a recent updating of these definitions (Lohman, 1972) has discarded this formal differentiation. It is sufficient to note that differences in the temperature of measurement between the field environment and the laboratory environment can influence hydraulic conductivity values through the viscosity term in Eq. (2.28). The effect is usually small, so correction factors are seldom introduced. It still makes good sense to report whether hydraulic conductivity measurements have been carried out in the laboratory or in the field, because the methods of measurement are very different and the interpretations placed on the values may be dependent on the type of measurement. However, this information is of practical rather than conceptual importance.

Table 2.2 indicates the range of values of hydraulic conductivity and permeability in five different systems of units for a wide range of geological materials. The table is based in part on the data summarized in Davis’ (1969) review. The primary conclusion that can be drawn from the data is that hydraulic conductivity varies over a very wide range. There are very few physical parameters that take on values over 13 orders of magnitude. In practical terms, this property implies that an order-of-magnitude knowledge of hydraulic conductivity can be very useful. Conversely, the third decimal place in a reported conductivity value probably has little significance.

Table 2.3 provides a set of conversion factors for the various common units of k and K. As an example of its use, note that a k value in cm² can be converted to one in ft² by multiplying by 1.08 × 10^–3. For the reverse conversion from ft² to cm², multiply by 9.29 × 10².

The various approaches to the measurement of hydraulic conductivity in the laboratory and in the field are described in Sections 8.4 through 8.6.

Table 2.2 Range of Values of Hydraulic Conductivity and Permeability

Table 2.3 Conversion Factors for Permeability and Hydraulic Conductivity Units

	Permeability, k*			Hydraulic conductivity, K
	cm²	ft²	darcy	m/s	ft/s	U.S. gal/day/ft²
cm²	1	1.08 × 10^–3	1.01 × 10⁸	9.80 × 10²	3.22 × 10²	1.85 × 10⁹
ft²	9.29 × 10²	1	9.42 × 10¹⁰	9.11 × 10⁵	2.99 × 10⁶	1.71 × 10¹²
darcy	9.87 × 10^–9	1.06 × 10^–11	1	9.66 × 10^–6	3.17 × 10^–5	1.82 × 10¹
m/s	1.02 × 10^–3	1.10 × 10^–6	1.04 × 10⁵	1	3.28	2.12 × 10⁶
ft/s	3.11 × 10^–4	3.35 × 10^–7	3.15 × 10⁴	3.05 × 10^–1	1	6.46 × 10⁵
U.S. gal/day/ft²	5.42 × 10^–10	5.83 × 10^–13	5.49 × 10^–2	4.72 × 10^–7	1.55 × 10^–6	1

*To obtain k in ft², multiply k in cm² by 1.08 × 10^–3.

2.4 Heterogeneity and Anisotropy of Hydraulic Conductivity

Hydraulic conductivity values usually show variations through space within a geologic formation. They may also show variations with the direction of measurement at any given point in a geologic formation. The first property is termed heterogeneity and the second anisotropy. The evidence that these properties are commonplace is to be found in the spread of measurements that arises in most field sampling programs. The geological reasoning that accounts for their prevalence lies in an understanding of the geologic processes that produce the various geological environments.

Homogeneity and Heterogeneity

If the hydraulic conductivity K is independent of position within a geologic formation, the formation is homogeneous. If the hydraulic conductivity K is dependent on position within a geologic formation, the formation is heterogeneous. If we set up an xyz coordinate system in a homogeneous formation K(x, y, z) = C, C being a constant; whereas in a heterogeneous formation, K(x, y, z) ≠ C.

There are probably as many types of heterogeneous configurations as there are geological environments, but it may be instructive to draw attention to three broad classes. Figure 2.7(a) is a vertical cross section that shows an example of layered heterogeneity, common in sedimentary rocks and unconsolidated lacustrine and marine deposits. Here, the individual beds making up the formation each have a homogeneous conductivity value K₁, K₂, . . . , but the entire system can be thought of as heterogeneous. Layered heterogeneity can result in K contrasts of almost the full 13-order range (Table 2.2), as, for example, in interlayered deposits of clay and sand. Equally large contrasts can arise in cases of discontinuous heterogeneity caused by the presence of faults or large-scale stratigraphic features. Perhaps the most ubiquitous discontinuous feature is the overburden-bedrock contact. Figure 2.7(b) is a map that shows a case of trending heterogeneity. Trends are possible in any type of geological formation, but they are particularly common in response to the sedimentation processes that create deltas, alluvial fans, and glacial outwash plains. The A, B, and C soil horizons often show vertical trends in hydraulic conductivity, as do rock types whose conductivity is primarily dependent on joint and fracture concentration reading heterogeneity in large consolidated or unconsolidated sedimentary formations can attain gradients of 2–3 orders of magnitude in a few miles.

Many hydrogeologists and petroleum geologists have used statistical distributions to provide a quantitative description of the degree of heterogeneity in a geological formation. There is now a large body of direct evidence to support the statement that the probability density function for hydraulic conductivity is log-normal. Warren and Price (1961) and Bennion and Griffiths (1966) found this to be the case in oilfield reservoir rocks, and Willardson and Hurst (1965) and Davis (1969) support the conclusion for unconsolidated water-bearing formations. A log-normal distribution for K is one for which a parameter Y, defined as Y = log K, shows a normal distribution. Freeze (1975) provides a table, based on the references above, that shows the standard deviation on Y (which is independent of the units of measurement) is usually in the range 0.5–1.5. This means that K values in most geological formations show internal heterogeneous variations of 1–2 orders of magnitude. Trending heterogeneity within a geological formation can be thought of as a trend in the mean value of the probability distribution. The same standard deviation may be evident in measurements at different positions in the formation, but the trending means lead to an increase in the overall observed range for the formation.

**Figure 2.7** Layered heterogeneity and trending heterogeneity.

Greenkorn and Kessler (1969) have provided a set of definitions of heterogeneity that are consistent with the statistical observations. In effect, they argue that if all geologic formations display spatial variations in K, then under the classical definitions, there is no such thing as a homogeneous formation. They redefine a homogeneous formation as one in which the probability density function of hydraulic conductivity is monomodal. That is, it shows variations in K, but maintains a constant mean K through space. A heterogeneous formation is defined as one in which the probability density function is multimodal. To describe a porous medium that satisfies the classical definition of homogeneity (K constant every. where. such as in experimental glass beads of diameter d) they use the term uniform. If we wish to adapt the classical definitions given at the start of this section to this more rational set of concepts, we can do so by adding the adjective “mean” and coaching the original definitions in terms of mean hydraulic conductivity.

Isotropy and Anisotropy

If the hydraulic conductivity K is independent of the direction of measurement at a point in a geologic formation, the formation is isotropic at that point. If the hydraulic conductivity K varies with the direction of measurement at a point in a geologic formation, the formation is anisotropic at that point.

Consider a two-dimensional vertical section through an anisotropic formation. If we let θ be the angle between the horizontal and the direction of measurement of a K value at some point in the formation, then K = K(θ). The directions in space corresponding to the θ angle at which K attains its maximum and minimum values are known as the principal directions of anisotropy. They are always perpendicular to one another. In three dimensions, if a plane is taken perpendicular to one of the principal directions, the other two principal directions are the directions of maximum and minimum K in that plane.

If an xyz coordinate system is set up in such a way that the coordinate directions coincide with the principal directions of anisotropy, the hydraulic conductivity values in the principal directions can be specified a K_x, K_y, K_z. At any point (x, y, z), an isotropic formation will have K_x = K_y = K_z whereas an anisotropic formation will have K_x ≠ K_y ≠ K_z. If K_x ≠ K_y ≠ K_z as is common in horizontally bedded sedimentary deposits, the formation is said to be transversely isotropic.

To fully describe the nature of the hydraulic conductivity in a geologic formation, it is necessary to use two adjectives, one dealing with heterogeneity and one with anisotropy. For example, for a homogeneous, isotropic system in two dimensions: K_x(x, z) = K_z(x, z) = C for all (x, z), where C is a constant. For a homogeneous, anisotropic system, K_x(x, z) = C₁ for all (x, z) and K_z(x, z) = C₂ for all (x, z) but C₁ ≠ C₂. Figure 2.8 attempts to further clarify the four possible combinations. The length of the arrow vectors is proportional to the K_x and K_z values at the two points (x₁, z₁) and (x₂, z₂).

The primary cause of anisotropy on a small scale is the orientation of clay minerals in sedimentary rocks and unconsolidated sediments. Core samples of clays and shales seldom show horizontal to vertical anisotropy greater than 10:1, and it is usually less than 3:1.

On a larger scale, it can be shown (Maasland, 1957; Marcus and Evenson, 1961) that there is a relation between layered heterogeneity and anisotropy. Consider the layered formation shown in Figure 2.9. Each layer is homogeneous and isotropic with hydraulic conductivity values K₁, K₂, . . . , K_n. We will show that the system as a whole acts like a single homogeneous, anisotropic layer. First, consider flow perpendicular to the layering.

**Figure 2.8** Four possible combinations of heterogeneity and anisotropy.

**Figure 2.9** Relation between layered heterogeneity and anisotropy.

The specific discharge v must be the same entering the system as it is leaving; in fact, it must be constant throughout the system. Let Δh₁ be the head loss across the first layer, Δh₂ across the second layer, and so on. The total head loss is then Δh = Δh₁ + Δh₂ + . . . + Δh_n, and from Darcy’s law,

$v = \frac{K_1\Delta h_1}{d_1} = \frac{K_2\Delta h_2}{d_2} = ... = \frac{K_n\Delta h_n}{d_n} = \frac{K_z\Delta h}{d}$ (2.30)

where K_z, is an equivalent vertical hydraulic conductivity for the system of layers. Solving the outside relationship of Eq. (2.30) for K_z, and using the inside relationships for Δh₁, Δh₂, . . . , we obtain

$K_z = \frac{vd}{\Delta h} = \frac{vd}{\Delta h_1 + \Delta h_2 + ... \Delta h_n} = \frac{vd}{vd_1/K_1+vd_2/K_2+...+vd_n/K_n}$

which leads to

$K_z = \frac{d}{\sum_{i=1}^n\frac{d_i}{K_i}}$ (2.31)

Now consider how parallel to the layering. Let Δh be the head loss over a horizontal distance l. The discharge Q through a unit thickness of the system is the sum of the discharges through the layers. The specific discharge v = Q/d is therefore given by

$v = \displaystyle\sum_{i=1}^n \frac{K_id_i}{d}\frac{\Delta h}{l} = K_x\frac{\Delta h}{l}$

where K_x, is an equivalent horizontal hydraulic conductivity. Simplification yields

$K_x = \sum_{i=1}^n\frac{K_id_i}{d}$ (2.32)

Equations (2.31) and (2.32) provide the K_x and K_z values for a single homogeneous but anisotropic formation that is hydraulically equivalent to the layered system of homogeneous, isotropic geologic formations of Figure 2.9. With some algebraic manipulation of these two equations it is possible to show that K_x > K_z for all possible sets of values of K₁, K₂, . . . , K_n. In fact, if we consider a set of cyclic couplets K₁, K₂, K₁, K₂, . . . with K₁ = 10⁴ and K₂ = 10², then K_x/K_z = 25. For K₁ = 10⁴ and K₂ = 1, K_x/K_z = 2500. In the field, it is not uncommon for layered heterogeneity to lead to regional anisotropy values on the order of 100:1 or even larger.

Snow (1969) showed that fractured rocks also behave anisotropic because of the directional variations in joint aperture and spacing. In this case, it is quite common for K_z > K_x.

Darcy’s Law in Three Dimensions

For three-dimensional flow, in a medium that may be anisotropic, it is necessary to generalize the one-dimensional form of Darcy’s law [Eq. (2.3)] presented earlier. In three dimensions the velocity v is a vector with components v_x, v_y, and v_z, and the simplest generalization would be

$v_x =-K_x\frac{\partial h}{\partial x}$

$v_y =-K_y\frac{\partial h}{\partial y}$ (2.33)

$v_z =-K_z\frac{\partial h}{\partial z}$

where K_x, K_y, and K_z, are the hydraulic conductivity values in the x, y, and z direction. Since h is now a function of x, y, and z, the derivatives must be partial.

In this text we will assume this simple generalization to be an adequate description of three-dimensional flow, but it is worth noting that a more generalized set of equations could be written in the form

$v_x = -K_{xx}\frac{\partial h}{\partial x} - K_{xy}\frac{\partial h}{\partial y} - K_{xz}\frac{\partial h}{\partial z}$

$v_y = -K_{yx}\frac{\partial h}{\partial x} - K_{yy}\frac{\partial h}{\partial y} - K_{yz}\frac{\partial h}{\partial z}$ (2.34)

$v_x = -K_{zx}\frac{\partial h}{\partial x} - K_{zy}\frac{\partial h}{\partial y} - K_{zz}\frac{\partial h}{\partial z}$

This set of equations exposes the fact that there are actually nine components of hydraulic conductivity in the most general case. If these components are put in matrix form, they form a second-rank symmetric tensor known as the hydraulic conductivity tensor (Bear, 1972). For the special case K_xy = K_yx = K_yz = K_zx = K_zy = 0, the nine components reduce to three and Eq. (2.33) is a suitable generalization of Darcy’s law. The necessary and sufficient condition that allows use of Eq. (2.33) rather than Eq. (2.34) is that the principal directions of anisotropy coincide with the x, y, and z coordinate axes. In most cases it is possible to choose a coordinate system that satisfies this requirement, but one can conceive of heterogeneous anisotropic systems in which the principal directions of anisotropy vary from one formation to another, and in such systems the choice of suitable axes would be impossible.

Hydraulic Conductivity Ellipsoid

Consider an arbitrary flowline in the xz plane in a homogeneous, anisotropic medium with principal hydraulic conductivities K_x and K_z, [Figure 2.10(a)].

**Figure 2.10** (a) Specific discharge *v_s* in an arbitrary direction of flow. (b) Hydraulic conductivity ellipse.

Along the flow line

$v_s = - K_s\frac{\partial h}{\partial s}$ (2.35)

where K_s is unknown, although it presumably lies in the range K_x – K_z. We can separate v_s into its components v_x and v_z, where

$v_x = -K_x\frac{\partial h}{\partial x} = v_s \cos \theta$
(2.36)

$v_z = -K_z\frac{\partial h}{\partial z} = v_s \sin \theta$

Now, since h = h(x, z),

$\frac{\partial h}{\partial s} = \frac{\partial h}{\partial x} \cdot \frac{\partial x}{\partial s} + \frac{\partial h}{\partial z} \cdot \frac{\partial z}{\partial s}$ (2.37)

Geometrically, ∂x/∂s = cos θ and ∂z/∂s = sin θ. Substituting these relationships together with Eqs. (2.35) and (2.36) in Eq. (2.37) and simplifying yields

$\frac{1}{K_s} = \frac{\cos^2 \theta}{K_x} + \frac{\sin^2 \theta}{K_z}$ (2.38)

This equation relates the principal conductivity components K_x and K_z to the resultant K_s in any angular direction θ. If we put Eq. (2.38) into rectangular coordinates by setting x = r cos θ and z = r sin θ, we get

$\frac{r^2}{K_s} = \frac{x^2}{K_x} + \frac{z^2}{K_z}$ (2.39)

which is the equation of an ellipse with major axes $\sqrt{K_x}$ and $\sqrt{K_z}$ [Figure 2.10(b)]. In three dimensions, it becomes an ellipsoid with major axes $\sqrt{K_x}$ , $\sqrt{K_y}$ and $\sqrt{K_z}$ , and it is known as the hydraulic conductivity ellipsoid. In Figure 2.10(b), the conductivity value K_s for any direction of flow in an anisotropic medium can be determined graphically if K_x and K_z are known.

In Section 5.1, the construction of flow nets in anisotropic media will be discussed, and it will be shown there that, in contrast to isotropic media, flowlines are not perpendicular to equipotential lines in anisotropic media.

2.5 Porosity and Void Ratio

If the total unit volume V_T of a soil or rock is divided into the volume of the solid portion V_s and the volume of the voids V_v, the porosity n is defined as n = V_v/V_T. It is usually reported as a decimal fraction or a percent.

Figure 2.11 shows the relation between various rock and soil textures and porosity. It is worth distinguishing between primary porosity, which is due to the soil or rock matrix [Figure 2.11(a), (b), (c), and (d)], and secondary porosity, which may be due to such phenomena as secondary solution [Figure 2.11(e)] or structurally controlled regional fracturing [Figure 2.11(f)].

**Figure 2.11** Relation between texture and porosity. (a) Well-sorted sedimentary deposit having low porosity; (b) poorly sorted sedimentary deposit having low porosity; (c) well-sorted sedimentary deposit; (d) well-sorted sedimentary deposit whose porosity has been diminished by the deposition of mineral matter in the interstices; (e) rock rendered porous by solution; (f) rock rendered porous by fracturing (after Meinzer, 1923).

Table 2.4, based in part on data summarized by Davis (1969), lists representative porosity ranges for various geologic materials. In general, rocks have lower porosities than soils; gravels, sands, and silts, which are made up of angular and rounded particles, have lower porosities than soils rich in platy clay minerals; an poorly sorted deposits [Figure 2.11(b)] have lower porosities than well-sorted deposits [Figure 2.11(a)].

Table 2.4 Range of Values of Porosity

	n(%)
Unconsolidated deposits
Gravel	25–40
Sand	25–50
Silt	35–50
Clay	40–70
Rocks
Fractured basalt	5–50
Karst limestone	5–50
Sandstone	5–30
Limestone, dolomite	0–20
Shale	0–10
Fractured crystalline rock	0–10
Dense crystalline rock	0–5

Source: Davis, 1969.

The porosity n can be an important controlling influence on hydraulic conductivity K. In sampling programs carried out within deposits of well-sorted sand or in fractured rock formations, samples with higher n generally also have higher K. Unfortunately, the relationship does not hold on a regional basis across the spectrum of possible rock and soil types. Clay-rich soils, for example, usually have higher porosities than sandy or gravelly soils but lower hydraulic conductivities. In section 8.7 techniques will be presented for the estimation of hydraulic conductivity from porosity and from grain-size analyses.

The porosity n is closely related to the void ratio e, which is widely used in soil mechanics. The void ratio is defined as e = V_v/V_s, and e is related to n by

$e = \frac{n}{1-n} \hspace{1cm} \text{or} \hspace{1cm} n = \frac{e}{1+e}$ (2.40)

Values of e usually fall in the range 0–3.

The measurement of porosity on soil samples in the laboratory will be treated in Section 8.4.

2.6 Unsaturated Flow and the Water Table

Up until this point, Darcy’s law and the concepts of hydraulic head and hydraulic conductivity have been developed with respect to a saturated porous medium, that is, one in which all the voids are filled with water. It is clear that some soils, especially those near the ground surface, are seldom saturated. Their voids are usually only partially filled with water, the remainder of the pore space being taken up by air. The flow of water under such conditions is termed unsaturated or partially saturated. Historically, the study of unsaturated flow has been the domain of soil physicists and agricultural engineers, but recently both soil scientists and groundwater hydrologists have recognized the need to pool their talents in the development of an integrated approach to the study of subsurface flow, both saturated and unsaturated.

Our emphasis in this section will be on the hydraulics of liquid-phase transport of water in the unsaturated zone. We will not discuss vapor-phase transport, nor will we consider soil water-plant interactions. These latter topics are of particular interest in the agricultural sciences and they play an important role in the interpretation of soil geochemistry. More detailed consideration of the physics and chemistry of moisture transfer in unsaturated soils can be found at an introductory level in Baver et al. (1972) and at a more advanced level in Kirkham and Powers (1972) and Childs (1969).

Moisture Content

If the total unit volume V_T of a soil or rock is divided into the volume of the solid portion V_s, the volume of the water V_w, and the volume of the air V_a, the volumetric moisture content θ is defined as θ = V_w/V_T. Like the porosity n, it is usually reported as a decimal fraction or a percent. For saturated flow, θ = n; for unsaturated flow, θ < n.

Water Table

The simplest hydrologic configuration of saturated and unsaturated conditions is that of an unsaturated zone at the surface and a saturated zone at depth (Figure 2.12(a)]. We commonly think of the water table as being the boundary between the two, yet we are aware that a saturated capillary fringe often exists above the water table. With this type of complication lurking in the background, we must take care to set up a consistent set of definitions for the various saturated-unsaturated concepts.

The water table is best defined as the surface on which the fluid pressure p in the pores of a porous medium is exactly atmospheric. The location of this surface is revealed by the level at which water stands in a shallow well open along its length and penetrating the surficial deposits just deeply enough to encounter standing water in the bottom. If p is measured in gage pressure, then on the water table, p = 0. This implies ψ = 0, and since h = ψ + z, the hydraulic head at any point on the water table must be equal to the elevation z of the water table at that point. On figures we will often indicate the position of the water table by means of a small inverted triangle, as in Figure 2.12(a).

**Figure 2.12** Groundwater conditions near the ground surface. (a) Saturated and unsaturated zones; (b) profiles of moisture content versus depth; (c) pressure-head and hydraulic-head relationships; insets: water retention under pressure heads less than (top) and greater than (bottom) atmospheric; (d) profile of pressure head versus depth; (e) profile of hydraulic head versus depth.

Negative Pressure Heads and Tensiometers

We have seen that ψ > 0 (as indicated by piezometer measurements) in the saturated zone and that ψ = 0 on the water table. It follows that ψ < 0 in the unsaturated zone. This reflects the fact that water in the unsaturated zone is held in the soil pores under surface-tension forces. A microscopic inspection would reveal a concave meniscus extending from grain to grain across each pore channel [as shown in the upper circular inset on Figure 2.12(c)]. The radius of curvature on each meniscus reflects the surface tension on that individual, microscopic air-water interface. In reference to this physical mechanism of water retention, soil physicists often call the pressure head ψ, when ψ < 0, the tension head or suction head. In this text, on the grounds that one concept deserves only one name, we will use the term pressure head to refer to both positive and negative ψ.

Regardless of the sign of ψ, the hydraulic head h is still equal to the algebraic sum of ψ and z. However, above the water table, where ψ < 0, piezometers are no longer a suitable instrument for the measurement of h. Instead, h must be obtained indirectly from measurements of ψ determined with tensiometers. Kirkham (1964) and S. J. Richards (1965) provide detailed descriptions of the design and use of these instruments. Very briefly, a tensiometer consists of a porous cup attached to an airtight, water-filled tube. The porous cup is inserted into the soil at the desired depth, where it comes into contact with the soil water and reaches hydraulic equilibrium. The equilibration process involves the passage of water through the porous cup from the tube into the soil. The vacuum created at the top of the airtight tube is a measure of the pressure head in the soil. It is usually measured by a vacuum gage attached to the tube above the ground surface, but it can be thought of as acting like the inverted manometer shown for point 1 in the soil profile of Figure 2.12(c). To obtain the hydraulic head h, the negative ψ value indicated by the vacuum gage on a tensiometer must be added algebraically to the elevation z of the point of measurement. In Figure 2.12(c) the instrument at point 1 is a tensiometer; the one at point 3 is a piezometer. The diagram is, of course, schematic. In practice, the tensiometer would be a tube with a gage and a porous cup at the base; the piezometer would be an open pipe with a well point at the base.

Characteristic Curves of the Unsaturated Hydraulic Parameters

There is a further complication to the analysis of flow in the unsaturated zone. Both the moisture content θ and the hydraulic conductivity K are functions of the pressure head ψ. On reflection, the first of these conditions should come as no great surprise. In that soil moisture is held between the soil grains under surface-tension forces that are reflected in the radius of curvature of each meniscus, we might expect that higher moisture contents would lead to larger radii of curvature, lower surface-tension forces, and lower tension heads (i.e., less-negative pressure heads). Further, it has been observed experimentally that the θ – ψ relationship is hysteretic; it has a different shape when soils are wetting than when they are drying. Figure 2.13(a) shows the hysteretic functional relationship between θ and ψ for a naturally occurring sand soil (after Liakopoulos, 1965a). If a sample of this soil were saturated at a pressure head greater than zero and the pressure was then lowered step by step until it reached levels much less than atmospheric $(\psi \ll \theta)$ the moisture contents at each step would follow the drying curve (or drainage curve) on Figure 2.13(a). If water were then added to the dry soil in small steps, the pressure heads would take the return route along the wetting curve (or imbibition curve). The internal lines are called scanning curves. They show the course that θ and ψ would follow if the soil were only partially wetted, then dried, or vice versa.

One would expect, on the basis of what has been presented thus far, that the moisture content θ would equal the porosity n for all ψ > 0. For coarse-grained soils this is the case, but for fine-grained soils this relationship holds over a slightly larger range ψ > ψ_a, where ψ_a is a small negative pressure head [Figure 2.13(a)] known as the air entry pressure head. The corresponding pressure p_a is called the air entry pressure or the bubbling pressure.

Figure 2.13(b) displays the hysteretic curves relating the hydraulic conductivity K to the pressure head ψ for the same soil. For ψ > ψ_a, K = K₀, where K₀ is now known as the saturated hydraulic conductivity. Since K = K(ψ) and θ = θ(ψ), it is also true that K = K(θ). The curves of Figure 2.13(b) reflect the fact that the hydraulic conductivity of an unsaturated soil increases with increasing moisture content. If we write Darcy’s law for unsaturated flow in the x direction in an isotropic soil as

$v_x = -K(\psi)\frac{\partial h}{\partial x}$ (2.41)

we see that the existence of the K(ψ) relationship implies that, given a constant hydraulic gradient, the specific discharge v increases with increasing moisture content.

**Figure 2.13** Characteristic curves relating hydraulic conductivity and moisture content to pressure head for a naturally occurring sand soil (after Liakopoulos, 1956a).

In actual fact, it would be impossible to hold the hydraulic gradient constant while increasing the moisture content. Since h = ψ + z and θ(ψ), the hydraulic head h is also affected by the moisture content. In other words, a hydraulic-head gradient infers a pressure-head gradient (except in pure gravity flow), and this in turn infers a moisture-content gradient. In Figure 2.12, the vertical profiles for these three variables are shown schematically for a hypothetical case of downward infiltration from the surface. Flew must be downward because the hydraulic heads displayed in Figure 2.12(e) decrease in that direction. The large positive values of h infer that $|z| \gg |\psi|$ . In other words, the z = 0 datum lies at some depth. For a real case, these three profiles would be quantitatively interlinked through θ(ψ) and K(ψ) curves for the soil at the site. For example, if the θ(ψ) curve were known for the soil and the θ(z) profile measured in the field, then the ψ(z) profile, and hence the h(z) profile, could be calculated.

The pair of curves θ(ψ) and K(ψ) shown in Figure 2.13 are characteristic for any given soil. Sets of measurements carried out on separate samples from the same homogeneous soil would show only the usual statistical variations associated with spatially separated sampling points. The curves are often called the characteristic curves. In the saturated zone we have the two fundamental hydraulic parameters K₀ and n; in the unsaturated zone these become the functional relationships K(ψ) and θ(ψ). More succinctly,

$\theta = \theta(\psi) \hspace{1cm} \psi < \psi_a$
(2.42)

$\theta = n \hspace{1cm} \psi \geq \psi_a$

$K = K(\psi) \hspace{1cm} \psi < \psi_a$
(2.43)

$K = K_0 \hspace{1cm} \psi \geq \psi_a$

Figure 2.14 shows some hypothetical single-valued characteristic curves (i.e., without hysteresis) that are designed to show the effect of soil texture on the shape of the curves. For a more complete description of the physics of moisture retention in unsaturated soils, the reader is directed to White et al. (1971).

**Figure 2.14** Single-valued characteristic curves for three hypothetical soils. (a) Uniform sand; (b) silty sand; (c) silty clay.

Saturated, Unsaturated, and Tension-Saturated Zones

It is worthwhile at this point to summarize the properties of the saturated and unsaturated zones as they have been unveiled thus far. For the saturated zone, we can state that:

It occurs below the water table.
The soil pores are filled with water, and the moisture content θ equals the porosity n.
The fluid pressure p is greater than atmospheric, so the pressure head ψ (measured as gage pressure) is greater than zero.
The hydraulic head h must be measured with a piezometer.
The hydraulic conductivity K is a constant; it is not a function of the pressure head ψ.

For the unsaturated zone (or, as it is sometimes called, the zone of aeration or the vadose zone):

It occurs above the water table and above the capillary fringe.
The soil pores are only partially filled with water; the moisture content θ is less than the porosity n.
The fluid pressure p is less than atmospheric; the pressure head ψ is less than zero.
The hydraulic head h must be measured with a tensiometer.
The hydraulic conductivity K and the moisture content θ are both functions of the pressure head ψ.

In short, for saturated flow, ψ > 0, θ = n, and K = K₀; for unsaturated flow, ψ < 0, θ = θ(ψ), and K = K(ψ).

The capillary fringe fits into neither of the groupings above. The pores there are saturated, but the pressure heads are less than atmospheric. A more descriptive name that is now gaining acceptance is the tension-saturated zone. An explanation of its seemingly anomalous properties can be discovered in Figure 2.13. It is the existence of the air entry pressure head ψ_a < 0 on the characteristic curves that is responsible for the existence of a capillary fringe. ψ_a is the value of ψ that will exist at the top of the tension-saturated zone, as shown by ψ_A for point A in Figure 2.12(d). Since ψ_A has greater negative values in clay soils than it does in sands, these fine-grained soils develop thicker tension-saturated zones than do coarse-grained soils.

Some authors consider the tension-saturated zone as part of the saturated zone, but in that case the water table is no longer the boundary between the two zones. From a physical standpoint it is probably best to retain all three zones—saturated, tension-saturated, and unsaturated—in one’s conception of the complete hydrologic system.

A point that follows directly from the foregoing discussion in this section may warrant a specific statement. In that fluid pressures are less than atmospheric, there can be no natural outflow to the atmosphere from an unsaturated or tension-saturated face. Water can be transferred from the unsaturated zone to the atmosphere by evaporation and transpiration, but natural outflows, such as springs on streambanks or inflows to well bores, must come from the saturated zone. The concept of a saturated seepage face is introduced in Section 5.5 and its importance its relation to hillslope hydrology is emphasized in Section 6.5.

Perched and Inverted Water Tables

The simple hydrologic configuration that we have considered thus far, with a single unsaturated zone overlying the main saturated groundwater body, is a common one. It is the rule where homogeneous geologic deposits extend to some depth. Complex geological environments, on the other hand, can lead to more complex saturated-unsaturated conditions. The existence of a low-permeability clay layer in a high-permeability sand formation, for example, can lead to the formation of a discontinuous saturated lense, with unsaturated conditions existing both above and below. If we consider the line ABCDA in Figure 2.15 to be the ψ = 0 isobar, we can refer to the ABC portion as a perched water table and ADC as an inverted water table. EF is the true water table.

**Figure 2.15** Perched water table *ABC*, inverted water table *ADC*, and true water table EF.

Saturated conditions can be discontinuous in time as well as space. Heavy rainfall can lead to the formation of a temporary saturated zone at the ground surface, its lower boundary being an inverted water table underlain by unsaturated conditions. Saturated zones of this type dissipate with time under the influence of downward percolation and evaporation from the surface. In Chapter 6 we will examine the interactions of rainfall and infiltration in saturated-unsaturated systems in greater detail.

Multiphase Flow

The approach to unsaturated flow outlined in this section is the one used almost universally by soil physicists, but it is, at root, an approximate method. Unsaturated flow is actually a special case of multiphase flow through porous media, with two phases, air and water, coexisting in the pore channels. Let θ_w be the volumetric the moisture content (previously denoted by θ) and θ_a be the volumetric air content, defined analogously to θ_w. There are now two fluid pressures to consider: p_w for the water phase and p_a for the air phase; and two pressure heads, ψ_w and ψ_a. Each soil now possesses two characteristic curves of fluid content versus pressure head, one for the water, θ_w(ψ_w) and one for the air, θ_a(ψ_a).

When it comes to conductivity relationships, it makes sense to work with the permeability k [Eq. (2.28)] rather than the hydraulic conductivity K, since it is independent of the fluid and K is not. The flow parameters k_w and k_a are called the effective permeabilities of the medium to water and air. Each soil has two characteristic curves of effective permeability versus pressure head, one for water, k_w(ψ_w) and one for air, k_a(ψ_a).

The single-phase approach to unsaturated flow leads to techniques of analysis that are accurate enough for almost all practical purposes, but there are some unsaturated flow problems where the multiphase flow of air and water must be considered. These commonly involve cases where a buildup in air pressure in the entrapped air ahead of a wetting front influences the rate of propagation of the front through a soil. Wilson and Luthin (1963) encountered the effects experimentally, Youngs and Peck (1964) provide a theoretical discussion, and McWhorter (1971) presents a complete analysis. As will be shown in Section 6.8, air entrapment also influences water-table fluctuations. Bianchi and Haskell (1966) discuss air entrapment problems in a field context, and Green et al. (1970) describe a field application of the multiphase approach to the analysis of a subsurface flow system.

Much of the original research on multiphase flow through porous media was carried out in the petroleum industry. Petroleum reservoir engineering involves the analysis of three-phase flow of oil, gas, and water. Pirson (1958) and Amyx et al. (1960) are standard references in the field. Stallman (1964) provides an interpretive review of the petroleum multiphase contributions as they pertain to groundwater hydrology.

The two-phase analysis of unsaturated flow is an example of immiscible displacement; that is, the fluids displace each other without mixing, and there is a distinct fluid-fluid interface within each pore. The simultaneous flow of two fluids that are soluble in each other is termed miscible displacement, and in such cases a distinct fluid-fluid interface does not exist. Bear (1972) provides an advanced theoretical treatment of both miscible and immiscible displacement in porous media. In this text, the only examples of immiscible displacement are those that have been discussed in this subsection. In the rest of the text, unsaturated flow will be treated as a single-phase problem using the concepts and approach of the first part of this section. The most common occurrences of miscible displacement in groundwater hydrology involve the mixing of two waters with different chemistry (such as seawater and fresh-water, or pure water and contaminated water). The transport processes associated with miscible displacement and the techniques of analysis of groundwater contamination will be discussed in Chapter 9.

2.7 Aquifers and Aquitards

Of all the words in the hydrologic vocabulary, there are probably none with more shades of meaning than the term aquifer. It means different things to different people, and perhaps different things to the same person at different times. It is used to refer to individual geologic layers, to complete geologic formations, and even to groups of geologic formations. The term must always be viewed in terms of the scale and context of its usage.

Aquifers, Aquitards, and Aquicludes

An aquifer is best defined as a saturated permeable geologic unit that can transmit significant quantities of water under ordinary hydraulic gradients. An aquiclude is defined as a saturated geologic unit that is incapable of transmitting significant quantities of water under ordinary hydraulic gradients.

An alternative pair of definitions that are widely used in the water-well industry states that an aquifer is permeable enough to yield economic quantities of water to wells, whereas aquicludes are not.

In recent years the term aquitard has been coined to describe the less-permeable beds in a stratigraphic sequence. These beds may be permeable enough to transmit water in quantities that are significant in the study of regional groundwater flow, but their permeability is not sufficient to allow the completion of production wells within them. Most geologic strata are classified as either aquifers or aquitards; very few formations fit the classical definition of an aquiclude. As a result, there is a trend toward the use of the first two of these terms at the expense of the third.

The most common aquifers are those geologic formations that have hydraulic conductivity values in the upper half of the observed range (Table 2.2): unconsolidated sands and gravels, permeable sedimentary rocks such as sandstones and limestones, and heavily fractured volcanic and crystalline rocks. The most common aquitards are clays, shales, and dense crystalline rocks. In Chapter 4, the principal aquifer and aquitard types will be examined more fully within the context of a discussion on geological controls on groundwater occurrence.

The definitions of aquifer and aquitard are purposely imprecise with respect to hydraulic conductivity. This leaves open the possibility of using the terms in a relative sense. For example, in an interlayered sand-silt sequence, the silts may be considered aquitards, whereas in a silt-clay system, they are aquifers.

Aquifers are often called by their stratigraphic names. The Dakota Sandstone for example, owes its geological fame largely to Meinzer’s (1923) assessment of it properties as an aquifer. Two other well-known North American aquifers are the St. Peter Sandstone in Illinois and the Ocala Limestone in Florida. A summary of the principal aquifer systems in the United States can be found in McGuinnes (1963) and Maxey (1964), who build on the earlier compilations of Meinzer (1923) Tolman (1937), and Thomas (1951). Brown (1967) provides information on Canada’s major aquifers.

In the ideal world of analysis where many of the expository sections of this book must reside, aquifers tend to appear as homogeneous, isotropic formations of constant thickness and simple geometry. We hope the reader will bear in mind that the real world is somewhat different. The hydrogeologist constantly faces complex aquifer-aquitard systems of heterogeneous and anisotropic formations rather than the idealized cases pictured in texts. It will often seem that the geological processes have maliciously conspired to maximize the interpretive and analytic difficulties.

Confined and Unconfined Aquifers

A confined aquifer is an aquifer that is confined between two aquitards. An unconfined aquifer, or water-table aquifer, is an aquifer in which the water table forms the upper boundary. Confined aquifers occur at depth, unconfined aquifers near the ground surface (Figure 2.16). A saturated lense that is bounded by a perched water table (Figure 2.15) is a special case of an unconfined aquifer. Such lenses are sometimes called perched aquifers.

In a confined aquifer, the water level in a well usually rises above the top of the aquifer. If it does, the well is called an artesian well and the aquifer is said to exist under artesian conditions. In some cases the water level may rise above the ground surface, in which case the well is known as a flowing artesian well and the aquifer is said to exist under flowing artesian conditions. In Section 6.1, we will examine the topographic and geologic circumstances that lead to flowing artesian conditions. The water level in a well in an unconfined aquifer rests at the water table.

Potentiometric Surface

For confined aquifers, which are extensively tapped by wells for water supply, there has grown up a traditional concept that is not particularly sound but which is firmly entrenched in usage. If the water-level elevations in wells tapping a confined aquifer are plotted on a map and contoured, the resulting surface, which is actually a map of the hydraulic head in the aquifer, is called a potentiometric surface. A potentiometric map of an aquifer provides an indication of the directions of groundwater flow in the aquifer.

The concept of a potentiometric surface is only rigorously valid for horizontal flow in horizontal aquifers. The condition of horizontal flow is met only in aquifers with hydraulic conductivities that are much higher than those in the associated confining beds. Some hydrogeological reports contain potentiometric surface maps based on water-level data from sets of wells that bottom near the same elevation but that are not associated with a specific well-defined confined aquifer. This type of potentiometric surface is essentially a map of hydraulic head contours on a two-dimensional horizontal cross section taken through the three-dimensional hydraulic head pattern that exists in the subsurface in that area. If there are vertical components of flow, as there usually are, calculations and interpretations based on this type of potentiometric surface can be grossly misleading.

It is also possible to confuse a potentiometric surface with the water table in areas where both confined and unconfined aquifers exist. Figure 2.16 schematically distinguishes between the two. In general, as we shall see from the flow nets in Chapter 6, the two do not coincide.

2.8 Steady-State Flow and Transient Flow

Stead-state flow occurs when at any point in a flow field the magnitude and direction of the flow velocity are constant with time. Transient flow (or unsteady flow, or nonsteady flow) occurs when at any point in a flow field the magnitude or direction of the flow velocity changes with time.

Figure 2.17(a) shows a steady-state groundwater flow pattern (dashed equipotentials, solid flowline) through a permeable alluvial deposit beneath a concrete dam. Along the line AB, the hydraulic head h_AB = 1000 m. It is equal to the elevation of the surface of the reservoir above AB. Similarly, h_CD = 900 m (the elevation of the tailrace pond above CD). The hydraulic head drop Δh across the system is 100 m. If the water level in the reservoir above AB and the water level in the tailrace pond above CD do not change with time, the flow net beneath the dam will not change with time. The hydraulic head at point E, for example, will be h_E = 950 m and will remain constant. Under such circumstances the velocity $v = -K$ $\partial h/\partial l$ will also remain constant through time. In a steady-state flow system the velocity may vary from point to point, but it will not vary with time at any given point.

**Figure 2.17** Steady-state and transient groundwater flow beneath a dam.

Let us now consider the transient flow problem schematically shown in Figure 2.17(b). At time t₀ the flow net beneath the dam will be identical to that of Figure 2. 17(a) and h_E will be 950 m. If the reservoir level is allowed to drop over the period t₀ to t₁, until the water levels above and below the dam are identical at time t₁, the ultimate conditions under the dam will be static with no flow of water from the upstream to the downstream side. At point E the hydraulic head h_E will undergo a time-dependent decline from h_E = 950 m at time t₀ to its ultimate value of h_E = 900 m. There may well be time lag in such a system so that h_E will not necessarily reach the value h_E = 900 m until sometime after t = t₁.

One important difference between steady and transient systems lies in the relation between their flowlines and pathlines. Flowlines indicate the instantaneous directions of flow throughout a system (at all times in a steady system, or at a given instant in time in a transient system). They must be orthogonal to the equipotential lines throughout the region of flow at all times. Pathlines map the route that an individual particle of water follows through a region of flow during a steady or transient event. In a steady flow system a particle of water that enters the system at an inflow boundary will flow toward the outflow boundary along a pathline that coincides with a flowline such as that shown in Figure 2.17(a). In a transient flow system, on the other hand, pathlines and flowlines do not coincide. Although a flow net can be constructed to describe the flow conditions at any given instant in time in a transient system, the flowlines shown in such a snapshot represent only the directions of movement at that instant in time. In that the configuration of flowlines changes with time, the flowlines cannot describe, in themselves, the complete path of a particle of water as it traverses the system. The delineation of transient pathlines has obvious importance in the study of groundwater contamination.

A groundwater hydrologist must understand the techniques of analysis for both steady-state flow and transient flow. In the final sections of this chapter the equations of flow will be developed for each type of flow, under both saturated and unsaturated conditions. The practical methodology that is presented in later chapters is often based on the theoretical equations, but it is not usually necessary for the practicing hydrogeologist to have the mathematics at his or her fingertips. The primary application of steady-state techniques in groundwater hydrology is in the analysis of regional groundwater flow. An understanding of transient flow is required for the analysis of well hydraulics, groundwater recharge, and many of the geochemical and geotechnical applications.

2.9 Compressibility and Effective Stress

The analysis of transient groundwater flow requires the introduction of the concept of compressibility, a material property that describes the change in volume, or strain, induced in a material under an applied stress. In the classical approach to the strength of elastic materials, the modulus of elasticity is a more familiar material property. It is defined as the ratio of the change in stress d_σ to the resulting change in strain dε. Compressibility is simply the inverse of the modulus of elasticity. It is defined as strain/stress, dε/dσ, rather than stress/strain, dσ/dε. The term is utilized for both elastic materials and nonelastic materials. For the flow of water through porous media, it is necessary to define two compressibility terms, one for the water and one for the porous media.

Compressibility of Water

Stress is imparted to a fluid through the fluid pressure p. An increase in pressure dp leads to a decrease in the volume V_w of a given mass of water. The compressibility of water β is therefore defined as

$\beta = \frac{-dV_w/V_w}{dp}$ (2.44)

The negative sign is necessary if we wish β to be a positive number.

Equation (2.44) implies a linear elastic relationship between the volumetric strain dV_w/V_w and the stress induced in the fluid by the change in fluid pressure dp. The compressibility, β is thus the slope of the line relating strain to stress for water, and this slope does not change over the range of fluid pressures encountered in groundwater hydrology (including those less than atmospheric that are encountered in the unsaturated zone). For the range of groundwater temperatures that are usually encountered, temperature has a small influence on β so that for most practical purposes we can consider β a constant. The dimensions of β are the inverse of those for pressure or stress. Its value can be taken as 4.4 × 10^–10 m²/N (or Pa^–1).

For a given mass of water it is possible to rewrite Eq. (2.44) in the form

$\beta = \frac{dp/p}{dp}$ (2.45)

where ρ is the fluid density. Integration of Eq. (2.45) yields the equation of state for water:

$\rho = \rho_0 \text{exp}[\beta(p - p_0)]$ (2.46)

where ρ₀ is the fluid density at the datum pressure p₀. For p₀ atmospheric, Eq. (2.46) can be written in terms of gage pressures as

$\rho = \rho_0 e^{\beta p}$ (2.47)

An incompressible fluid is one for which β = 0 and ρ = ρ₀ = constant.

Effective Stress

Let us now consider the compressibility of the porous medium. Assume that a stress is applied to a unit mass of saturated sand. There are three mechanisms by which a reduction in volume can be achieved: (1) by compression of the water in the pores, (2) by compression of the individual sand grains, and (3) by a rearrangement of the sand grains into a more closely packed configuration. The first of these mechanisms is controlled by the fluid compressibility β. Let us assume that the second mechanism is negligible, that is, that the individual soil grains are incompressible. Our task is to define a compressibility term that will reflect the third mechanism.

To do so, we will have to invoke the principle of effective stress. This concept was first proposed by Terzaghi (1925), and has been analyzed in detail by Skempton (1961). Most soil mechanics texts, such as those by Terzaghi and Peck (1967) and Scott (1963), provide a full discussion.

For our purposes, consider the stress equilibrium on an arbitrary plane through a saturated geological formation at depth (Figure 2.18). σ_T is the total stress acting downward on the plane. It is due to the weight of overlying rock and water. This stress is borne in part by the granular skeleton of the porous medium and in part by the fluid pressure p of the water in the pores. The portion of the total stress that is not borne by the fluid is called the effective stress σ_e. It is this stress that is actually applied to the grains of the porous medium. Rearrangement of the soil grains and the resulting compression of the granular skeleton is caused by changes in the effective stress, not by changes in the total stress. The two are related by the simple equation

$\sigma_T = \sigma_e + p$ (2.48)

or, in terms of the changes,

$d\sigma_T = d\sigma_e + dp$ (2.49)

**Figure 2.18** Total stress, effective stress, and fluid pressure on an arbitrary plane through a saturated porous medium.

Many of the transient subsurface flow problems that must be analyzed do not involve changes in the total stress. The weight of rock and water overlying each point in the system often remains essentially constant through time. In such cases, dσ_T = 0 and

$d\sigma_e = -dp$ (2.50)

Under these circumstances, if the fluid pressure increases, the effective stress decreases by an equal amount; and if the fluid pressure decreases, the effective stress increases by an equal amount. For cases where the total stress does not change with time, the effective stress at any point in the system, and the resulting volumetric deformations there, are controlled by the fluid pressures at that point. Since p = ρgψ and ψ = h – z (z being constant at the point in question), changes in the effective stress at a point are in effect governed by changes in the hydraulic head at that point:

$d\sigma_e = -\rho g \hspace{1mm} d\psi = -\rho g \hspace{1mm} dh$ (2.51)

Compressibility of a Porous Medium

The compressibility of a porous medium is defined as

$\alpha = \frac{-dV_T/V_T}{d\sigma_e}$ (2.52)

where V_T is the total volume of a soil mass and dσ_e the change in effective stress.

Recall that V_T = V_S + V_v, where V_S is the volume of the solids and V_v is the volume of the water-saturated voids. An increase in effective stress dσ_e produces a reduction dV_T in the total volume of the soil mass. In granular materials this reduction occurs almost entirely as a result of grain rearrangements. It is true that individual grains may themselves be compressible, but the effect is usually considered to be negligible. In general, dV_T = dV_S + dV_v; but for our purposes we will assume that dV_S = 0 and dV_T = dV_v.

Consider a sample of saturated soil that has been placed in a laboratory loading cell such as the one shown in Figure 2.19(a). A total stress σ_T = L/A can be applied to the sample through the pistons. The sample is laterally confined by the cell walls, and entrapped water is allowed to escape through vents in the pistons to an external pool held at a constant known fluid pressure. The volumetric reduction in the size of the soil sample is measured at several values of L as L is increased in a stepwise fashion. At each step, the increased total stress dσ_T is initially borne by the water under increased fluid pressures, but drainage of the water from the sample to the external pool slowly transfers the stress from the water to the granular skeleton. This transient process is known as consolidation, and the time required for the consolidation process to reach hydraulic equilibrium at each L can be considerable. Once attained, however, it is known that dp = 0 within the sample, and from Eq. (2.49), dσ_e = dσ_T = dL/A. If the soil sample has an original void ratio e₀ (where e = V_v/V_S) and an original height b [Figure 2.19(a)], and assuming that dV_T = dV_v, Eq. (2.52) can be written as

$\alpha = \frac{-db/b}{d\sigma_e} = \frac{-de(1+e_0)}{d\sigma_e}$ (2.53)

The compressibility α is usually determined from the slope of a strain-stress plot in the form of e versus σ_e. The curve AB in Figure 2.19(b) is for loading (increasing σ_e), BC is for unloading (decreasing σ_e). In general, the strain-stress relation is neither linear nor elastic. In fact, for repeated loadings and unloadings, many fine-grained soils show hysteretic properties [Figure 2.19(c)]. The soil compressibility α, unlike the fluid compressibility β, is not a constant; it is a function of the applied stress and it is dependent on the previous loading history.

**Figure 2.19** (a) Laboratory loading cell for the determination of soil compressibility; (b), (c), and (d) schematic curves of void ratio versus effective stress.

Figure 2.19(d) provides a schematic comparison of the e – σ_e curves for clay and sand. The lesser slope for the sand curve implies a smaller α, and its linearity implies an α value that stays constant over a wide range of σ_e. In groundwater systems, the time-dependent fluctuations in σ_e are often quite small, so that even for clays, a constant α can have some meaning. Table 2.5 is a table of compressibility values that indicates the ranges of values that have been measured for various types of geologic materials.

Table 2.5 Range of Values of Compressibility*

	Compressibility, α (m²/N or Pa^–1)
Clay	10^–3–10^–8
Sand	10^–7–10^–9
Gravel	10^–8–10^–10
Jointed rock	10^–8–10^–10
Sound rock	10^–9–10^–11
Water (β)	4.4 × 10^–10

*See Table A1.3, Appendix I, for conversion factors.

Original sources of compressibility data include Domenico and Mifflin (1965) and Johnson et al. (1968). The dimensions of α like β, are the inverse of those for stress. Values are expressed in SI units of m²/N or Pa^–1. Note that the compressibility of water is of the same order of magnitude as the compressibility of the less-compressible geologic materials.

As noted in Figure 2.19(b) and (c), the compressibility of some soils in expansion (expansibility) is much less than in compression. For clays, the ratio of the two α’s is usually on the order of 10:1; for uniform sands, it approaches 1:1. For soils that have compressibility values that are significantly less in expansion than compression, volumetric deformations that occur in response to increasing effective stress [perhaps due to decreasing hydraulic heads as suggested by Eq. (2.51)] are largely irreversible. They are not recovered when effective stresses subsequently decrease. In a clay-sand aquifer-aquitard system, the large compactions that can occur in the clay aquitards (due to the large α values) are largely irrecoverable; whereas the small deformations that occur in the sand aquifers (due to the small α values) are largely elastic.

Aquifer Compressibility

The concept of compressibility inherent in Eq. (2.53) and in Figures 2.18 and 2.19 is one-dimensional. In the field, at depth, a one-dimensional concept has meaning if it is assumed that the soils and rocks are stressed only in the vertical direction. The total vertical stress σ_T at any point is due to the weight of overlying rock and water; the neighboring materials provide the horizontal confinement. The effective vertical stress σ_e is equal to σ_T – p. Under these conditions the aquifer compressibility α is defined by the first equality of Eq. (2.53), where b is now the aquifer thickness rather than a sample height. The parameter α is a vertical compressibility. If it is to be determined with a laboratory apparatus like that of Figure 2.19(a), the soil cores must be oriented vertically and loading must be applied at right angles to any horizontal bedding. Within an aquifer, α may vary with horizontal position; that is, on may be heterogeneous with α = (x, y).

In the most general analysis, it must be recognized that the stress field existing at depth is not one-dimensional but three-dimensional. In that case, aquifer compressibility must be considered as an anisotropic parameter. The vertical compressibility α is then invoked by changes in the vertical component of effective stress, and the horizontal compressibilities are invoked by changes in the horizontal components of effective stress. Application of the concepts of three-dimensional stress analysis in the consideration of fluid flow through porous media is an advanced topic that cannot be pursued here. Fortunately, for many practical cases, changes in the horizontal stress field are very small, and in most analyses it can be assumed that they are negligible. It is sufficient for our purposes to think of the aquifer compressibility α as a single isotropic parameter, but it should be kept in mind that it is actually the compressibility in the vertical direction, and that this is the only direction in which large changes in effective stress are anticipated.

To illustrate the nature of the deformations that can occur in compressible aquifers, consider the aquifer of thickness b shown in Figure 2.20. If the weight of overlying material remains constant and the hydraulic head in the aquifer is decreased by an amount -dh, the increase in effective stress dσ_e is given by Eq. (2.51) as ρg dh, and the aquifer compaction, from Eq. (2.53) is

$db = -\alpha b \hspace{1mm} d\sigma_e = -\alpha b \hspace{1mm} \rho g \hspace{1mm} dh$ (2.54)

The minus sign indicates that the decrease in head produces a reduction in the thickness b.

**Figure 2.20** Aquifer compaction caused by groundwater pumping.

One way that the hydraulic head might be lowered in an aquifer is by pumping from a well. Pumping induces horizontal hydraulic gradients toward the well in the aquifer, and as a result the hydraulic head is decreased at each point near the well. In response, effective stresses are increased at these points, and aquifer compaction results. Conversely, pumping water into an aquifer increases hydraulic heads, decreases effective stresses, and causes aquifer expansion. If the compaction of an aquifer-aquitard system due to groundwater pumping is propagated to the ground surface, the result is land subsidence. In Section 8.12 this phenomenon is considered in detail.

Effective Stress in the Unsaturated Zone

The first equality in Eq. (2.51) indicates that the relationship between the effective stress σ_e and the pressure head ψ should be linear. This relation, and the concept of Figure 2.18 on which it is based, holds in the saturated zone, but there is now abundant evidence to suggest that it does not hold in the unsaturated zone (Narasimhan, 1975). For unsaturated flow, Bishop and Blight (1963) suggest that Eq. (2.51) should be modified to read

$d\sigma_e = -\rho g\chi \hspace{1mm} d\psi$ (2.55)

where the parameter χ depends on the degree of saturation, the soil structure, and the wetting-drying history of the soil. Curve ABC in Figure 2.21 shows such a relationship schematically. For ψ > 0, χ = 1; for ψ < 0, χ ≤ 1; and for $\psi \ll 0$ , χ = 0.

**Figure 2.21** Relationship between effective stress and pressure head in the saturated and unsaturated zones (after Narasimhan, 1975).

The χ approach is an empirical one and its use reflects the fact that the capacity of fluid pressures less than atmospheric to support a part of the total stress in an unsaturated flow field is not yet fully understood. As a first approximation, it is not unreasonable to suppose that they have no such capacity, as suggested by curve ABD in Figure 2.21. Under this assumption, for ψ < 0, χ = 0, dσ_e = dσ_T, and changes in pressure head (or moisture content) in the unsaturated zone do not lead to changes in effective stress.

The definition of the compressibility of a porous medium in the unsaturated zone is still given by Eq. (2.52) just as it is in the saturated zone, but the influence of the fluid pressure on the effective stress is now considered to be muted or non-existent.

2.10 Transmissivity and Storativity

There are six basic physical properties of fluid and porous media that must be known in order to describe the hydraulic aspects of saturated groundwater flow. These six have all now been introduced. They are, for the water, density ρ, viscosity μ, and compressibility β; and for the media, porosity n (or void ratio e), permeability k, and compressibility α. All the other parameters that are used to describe the hydrogeologic properties of geologic formations can be derived from these six. For example, we have seen from Eq. (2.28) that the saturated hydraulic conductivity K is a combination of k, ρ, and μ. In this section, we will consider the concepts of specific storage S_s, storativity S, and transmissivity T.

Specific Storage

The specific storage S_s of a saturated aquifer is defined as the volume of water that a unit volume of aquifer releases from storage under a unit decline in hydraulic head. From Section 2.9 we now know that a decrease in hydraulic head h infers a decrease in fluid pressure p and an increase in effective stress σ_e. The water that is released from storage under conditions of decreasing h is produced by two mechanisms: (1) the compactions of the aquifer caused by increasing σ_e, and (2) the expansion of the water caused by decreasing p. The first of these mechanisms is controlled by the aquifer compressibility α and the second by the fluid compressibility β.

Let us first consider the water produced by the compaction of the aquifer. The volume of water expelled from the unit volume of aquifer during compaction will be equal to the reduction in volume of the unit volume of aquifer. The volumetric reduction dV_T will be negative, but the amount of water produced dV_w will be positive, so that, from Eq. (2.52),

$dV_W = -dV_T = \alpha V_Td\sigma_e$ (2.56)

For a unit volume, V_T = 1, and from Eq. (2.51), dσ_e = ρg dh. For a unit decline in hydraulic head, dh = -1, and we have

$dV_W = \alpha \rho g$ (2.57)

Now consider the volume of water produced by the expansion of the water. From Eq. (2.44),

$dV_W = -\beta V_W dp$ (2.58)

The volume of water V_w in the total unit volume V_T is nV_T, where n is the porosity. With V_T and dp = ρg dψ = ρg d(h – z) = ρg dh, Eq. (2.58) becomes, for dh = -1,

$dV_W =\beta n \rho g$ (2.59)

The specific storage S_s is the sum of the two terms given by Eqs. (2.57) and (2.59):

$S_s = \rho g (\alpha + n\beta)$ (2.60)

A dimensional inspection of this equation shows S_s to have the peculiar dimensions of [L]^-1. This also follows from the definition of S_s as a volume per volume per unit decline in head.

Transmissivity and Storativity of a Confined Aquifer

For a confined aquifer of thickness b, the transmissivity (or transmissibility) T is defined as

$T = Kb$ (2.61)

and the storativity (or storage coefficient) S is defined as

$S = S_sb$ (2.62)

If we substitute Eq. (2.60) in Eq. (2.62), the expanded definition of S is seen to be

$S = \rho gb(\alpha + n\beta)$ (2.63)

The storativity of a saturated confined aquifer of thickness b can be defined in words as the volume of water that an aquifer releases from storage per unit surface area of aquifer per unit decline in the component of hydraulic head normal to that surface. The hydraulic head for a confined aquifer is usually displayed in the form of a potentiometric surface, and Figure 2.22(a) illustrates the concept of storativity in this light.

**Figure 2.22** Schematic representation of the storativity in (a) confined and (b) unconfined aquifers (after Ferris et al., 1962).

In that the hydraulic conductivity K has dimensions [L/T], it is clear from Eq.(2.61) that the transmissivity T has dimensions [L²/T]. The SI metric unit is m²/s. T and S are widely used terms in the North American water well industry and are often expressed in FPS engineering units. If K is expressed in gal/day/ft², then T has units of gal/day/ft. The range of values of T can be calculated by multiplying the pertinent K values from Table 2.2 by the range of reasonable aquifer thicknesses, say 5–100 m. Transmissivities greater than 0.015 m²/s (or 0.16 ft²/s or 100,000 gal/day/ft) represent good aquifers for water well exploitation. Storativities are dimensionless. In confined aquifers, they range in value from 0.005 to 0.00005. Reference to the definition of S, coupled with a realization of its range of values, makes it clear that large head changes over extensive areas are required to produce substantial water yields from confined aquifers.

Transmissivities and storativities can be specified for aquitards as well as aquifers. However, in most applications, the vertical hydraulic conductivity of an aquitard has more significance than its transmissivity. It might also be noted that in clay aquitards, $\alpha \gg \beta$ , and the nβ term in the definition of storativity [Eq. (2.63)] and specific storage [Eq. (2.60)] becomes negligible.

It is possible to define a single formation parameter that couples the transmission properties T or K, and the storage properties S or S_s. The hydraulic diffusivity D is defined as

$D = \frac{T}{S}=\frac{K}{S_s}$ (2.64)

The term is not widely used in practice.

The concepts of transmissivity T and storativity S were developed primarily for the analysis of well hydraulics in confined aquifers. For two-dimensional, horizontal flow toward a well in a confined aquifer of thickness b, the terms are well defined; but they lose their meaning in many other groundwater applications. If a groundwater problem has three-dimensional overtones, it is best to revert to the use of hydraulic conductivity K and specific storage S_s; or perhaps even better, to the fundamental parameters permeabilityk, porosity n, and compressibility α.

Transmissivity and Specific Yield in Unconfined Aquifers

In an unconfined aquifer, the transmissivity is not as well defined as in a confined aquifer, but it can be used. It is defined by the same equation [Eq. (2.61)], but b is now the saturated thickness of the aquifer or the height of the water table above the top of the underlying aquitard that bounds the aquifer.

The storage term for unconfined aquifers is known as the specific yield S_y. It is defined as the volume of water that an unconfined aquifer releases from storage per unit surface area of aquifer per unit decline in the water table. It is sometimes called the unconfined storativity. Figure 2.22(b) illustrates the concept schematically.

The idea of specific yield is best visualized with reference to the saturated-unsaturated interaction it represents. Figure 2.23 shows the water-table position and the vertical profile of moisture content vs. depth in the unsaturated zone at two times, t₁ and t₂. The crosshatched area represents the volume of water released from storage in a column of unit cross section. If the water-table drop represents a unit decline, the crosshatched area represents the specific yield.

**Figure 2.23** Concept of specific yield viewed in terms of the unsaturated moisture profiles above the water table.

The specific yields of unconfined aquifers are much higher than the storativities of confined aquifers. The usual range of S_y is 0.01–0.30. The higher values reflect the fact that releases from storage in unconfined aquifers represent an actual dewatering of the soil pores, whereas releases from storage in confined aquifers represent only the secondary effects of water expansion and aquifer compaction caused by changes in the fluid pressure. The favorable storage properties of unconfined aquifers make them more efficient for exploitation by wells. When compared to confined aquifers, the same yield can be realized with smaller head changes over less extensive areas.

Storage in the Unsaturated Zone

In an unsaturated soil, changes in moisture content θ, such as those shown in Figure 2.23, are accompanied by changes in the pressure head ψ, through the θ(ψ) relationship displayed on the characteristic curve of Figure 2.13(a). The slope of this characteristic curve represents the unsaturated storage property of a soil. It is called the specific moisture capacity C and is defined as

$C = \frac{d\theta}{d\psi}$ (2.65a)

An increase of dψ in the pressure head (say, from –200 cm to –100 cm on Figure 2.13) must be accompanied by an increase of dθ in the moisture stored in the unsaturated soil. Since θ(ψ) is nonlinear and hysteretic, so too, is C. It is not a constant; it is a function of the pressure head ψ : C = C(ψ). In the saturated zone, in fact for all ψ > ψ_a, the moisture content θ is equal to the porosity n, a constant, so that C = 0. A parallel formulation to Eq. (2.42) for C is

$C = C(\psi) \hspace{1cm} \psi < \psi_a$
(2.65b)

$C = 0 \hspace{1cm} \psi \geq \psi_a$

The transmission and storage properties of an unsaturated soil are fully specified by the characteristic curve K(ψ) and one of the two curves θ(ψ) or C(ψ).

In an analogous manner to Eq. (2.64), the soil-water diffusivity can be defined as

$D(\psi) = \frac{K(\psi)}{C(\psi)}$ (2.66)

2.11 Equations of Groundwater Flow

In almost every field of science and engineering the techniques of analysis are based on an understanding of the physical processes, and in most cases it is possible to describe these processes mathematically. Groundwater flow is no exception. The basic law of flow is Darcy’s law, and when it is put together with an equation of continuity that describes the conservation of fluid mass during flow through a porous medium, a partial differential equation of flow is the result. In this section, we will present brief developments of the equations of flow for (l) steady-state saturated flow, (2) transient saturated flow, and (3) transient unsaturated flow. All three of the equations of flow are well known to mathematicians, and mathematical techniques for their manipulation are widely available and in common use in science and engineering. Generally, the equation of flow appears as one component of a boundary-value problem, so in the last part of this section we will explore this concept.

In that so many of the standard techniques of analysis in groundwater hydrology are based on boundary-value problems that involve partial differential equations, it is useful to have a basic understanding of these equations as one proceeds to learn the various techniques. Fortunately, it is not an absolute requirement. In most cases, the techniques can be explained and understood without returning at every step to the fundamental mathematics. The research hydrogeologist must work with the equations of flow on a daily basis; the practicing hydrogeologist can usually avoid the advanced mathematics if he or she so desires.

Steady-State Saturated Flow

Consider a unit volume of porous media such as that shown in Figure 2.24. Such an element is usually called an elemental control volume.

**Figure 2.24** Elemental control volume for flow through porous media.

The law of conservation of mass for steady-state flow through a saturated porous medium requires that the rate of fluid mass flow into any elemental control volume be equal to the rate of fluid mass flow out of any elemental control volume. The equation of continuity that translates this law into mathematical form can be written, with reference to Figure 2.24. as

$-\frac{\partial(\rho v_x)}{\partial x} -\frac{\partial(\rho v_y)}{\partial y} -\frac{\partial(\rho v_z)}{\partial z} = 0$ (2.67)

A quick dimensional analysis on the ρv terms will show them to have the dimensions of a mass rate of flow across a unit cross-sectional area of the elemental control volume. If the fluid is incompressible, ρ(x, y, z) = constant and the ρ’s can be removed from Eq. (2.67). Even if the fluid is compressible and ρ(x, y, z) ≠ constant, it can be shown that terms of the form ρ ∂v_x/∂_x are much greater than terms of the form v_x ∂ρ/∂x, both of which arise when the chain rule is used to expand Eq. (2.67). In either case, Eq. (2.67) simplifies to

$-\frac{\partial v_x}{\partial x} -\frac{\partial v_y}{\partial y} -\frac{\partial v_z}{\partial z} = 0$ (2.68)

Substitution of Darcy’s law for v_x, v_y, and v_z in Eq. (2.68) yields the equation of flow for steady-state flow through an anisotropic saturated porous medium:

$\frac{\partial}{\partial x}\left(K_x \frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y}\left(K_y \frac{\partial h}{\partial y}\right) + \frac{\partial}{\partial z}\left(K_z \frac{\partial h}{\partial z}\right) = 0$ (2.69)

For an isotropic medium, K_x = K_y = K_z and if the medium is also homogeneous, then K(x, y, z) = constant. Equation (2.69) then reduces to the equation of flow for steady-state flow through a homogeneous, isotropic medium:

$\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2} + \frac{\partial^2h}{\partial z^2} = 0$ (2.70)

Equation (2.70) is one of the most basic partial differential equations known to mathematicians. It is called Laplace’s equation. The solution of the equation is a function h(x, y, z) that describes the value of the hydraulic head h at any point in a three dimensional flow field. A solution to Eq. (2.70) allows us to produce a contoured equipotential map of h, and with the addition of flowlines, a flow net.

For steady-state, saturated flow in a two-dimensional flow field, say in the xy plane, the central term of Eq. (2.70) would drop out and the solution would be a function h(x, z).

Transient Saturated Flow

The law of conservation of mass for transient flow in a saturated porous medium requires that the net rate of fluid mass flow into any elemental control volume be equal to the time rate of change of fluid mass storage within the element. With reference to Figure 2.24, the equation of continuity takes the form

$-\frac{\partial(\rho v_x)}{\partial x} -\frac{\partial(\rho v_y)}{\partial y} -\frac{\partial(\rho v_z)}{\partial z} = \frac{\partial(\rho n)}{\partial t}$ (2.71)

or, expanding the right-hand side,

$-\frac{\partial(\rho v_x)}{\partial x} -\frac{\partial(\rho v_y)}{\partial y} -\frac{\partial(\rho v_z)}{\partial z} = n\frac{\partial \rho}{\partial t} + \rho\frac{\partial n}{\partial t}$ (2.72)

The first term on the right-hand side of Eq. (2.72) is the mass rate of water produced by an expansion of the water under a change in its density ρ. The second term is the mass rate of water produced by the compaction of the porous medium as reflected by the change in its porosity n. The first term is controlled by the compressibility of the fluid β and the second term by the compressibility of the aquifer α. We have already carried out the analysis (in Section 2.10) that is necessary to simplify the two terms on the right of Eq. (2.72). We know that the change in ρ and the change in n are both produced by a change in hydraulic head h, and that the volume of water produced by the two mechanisms for a unit decline in head is S_s, where S_s is the specific storage given by S_S = ρg(α + nβ). The mass rate of water produced (time rate of change of fluid mass storage) is ρS_S ∂h/∂t, and Eq. (2.72) becomes

$-\frac{\partial(\rho v_x)}{\partial x} -\frac{\partial(\rho v_y)}{\partial y} -\frac{\partial(\rho v_z)}{\partial z} = \rho S_s\frac{\partial h}{\partial t}$ (2.73)

Expanding the terms on the left-hand side by the chain rule and recognizing that terms of the form ρ ∂v_x/∂_x are much greater than terms of the form v_x∂ρ/∂x allows us to eliminate ρ from both sides of Eq. (2.73). Inserting Darcy’s law, we obtain

This is the equation of how for transient flow through a saturated anisotropic porous medium. If the medium is homogeneous and isotropic, Eq. (2.74) reduces to

$\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2} + \frac{\partial^2h}{\partial z^2} = \frac{S_s}{K}\frac{\partial h}{\partial t}$ (2.75)

or expanding S_s,

$\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2} + \frac{\partial^2h}{\partial z^2} = \frac{\rho g(\alpha+n\beta)}{K}\frac{\partial h}{\partial t}$ (2.76)

Equation (2.76) is known as the diffusion equation. The solution h(x, y, z, t) describes the value of the hydraulic head at any point in a flow field at any time. A solution requires knowledge of the three basic hydrogeological parameters, K, α, and n, and the fluid parameters, ρ and β.

For the special case of a horizontal confined aquifer of thickness b, S = S_sb and T = Kb, and the two-dimensional form of Eq. (2.75) becomes

$\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2} = \frac{S}{T}\frac{\partial h}{\partial t}$ (2.77)

The solution h(x, y, t) describes the hydraulic head field at any point on a horizontal plane through the horizontal aquifer at any time. Solution requires knowledge of the aquifer parameters S and T.

The equation of flow for transient, saturated flow [in any of the forms given by Eqs. (2.74) through (2.77)] rests on the law of flow established by Darcy (1856), on the clarification of the hydraulic potential by Hubbert (1940), and on the recognition of the concepts of aquifer elasticity by Meinzer (1923), and effective stress by Terzaghi (1925). The classical development was first put forward by Jacob (1940) and can be found in its most complete form in Jacob (1950). The development presented in this section, together with the storage concepts of the earlier sections, is essentially that of Jacob.

In recent years there has been considerable reassessment of the classical development. Biot (1955) recognized that in compacting aquifers it is necessary to cast Darcy’s law in terms of a relative velocity of fluid to grains, and Cooper (1966) pointed out the inconsistency of taking a fixed elemental control volume in a deforming medium. Cooper showed that Jacob’s classical development is correct if one views the velocity as relative and the coordinate system as deforming. He also showed that De Wiest’s (1966) attempt to address this problem (which also appears in Davis and De Wiest, 1966) is incorrect. Appendix II contains a more rigorous presentation of the Jacob-Cooper development than has been attempted here.

The classical development, through its use of the concept of vertical aquifer compressibility, assumes that stresses and deformations in a compacting aquifer occur only in the vertical direction. The approach couples a three-dimensional flow field and a one-dimensional stress field. The more general approach, which couples a three-dimensional flow field and a three-dimensional stress field, was first considered by Biot (1941, 1955). Verruijt (1969) provides an elegant summary of the approach.

For almost all practical purposes it is not necessary to consider relative velocities, deforming coordinates, or three-dimensional stress fields. The classical equations of flow presented in this section are sufficient.

Transient Unsaturated Flow

Let us define the degree of saturation θ’ as θ’ = θ/n, where θ is the moisture content and n is the porosity. For flow in an elemental control volume that may be only partially saturated, the equation of continuity must reveal the time rate of change of moisture content as well as the time rate of change of storage due to water expansion and aquifer compaction. The ρn term in Eq. (2.71) must become ρnθ’, and Eq. (2.72) becomes

$-\frac{\partial(\rho v_x)}{\partial x} -\frac{\partial(\rho v_y)}{\partial y} -\frac{\partial(\rho v_z)}{\partial z} = n\theta' \frac{\partial \rho}{\partial t} + \rho\theta' \frac{\partial n}{\partial t} + n\rho \frac{\partial \theta '}{\partial t}$ (2.78)

For unsaturated flow, the first two terms on the right-hand side of Eq. (2.78) are much smaller than the third term. Discarding these two terms, canceling the ρ’s from both sides in the usual way, inserting the unsaturated form of Darcy’s law [Eq. (2.41)], and recognizing that n dθ’ = dθ, leads to

$\frac{\partial}{\partial x}\left[K(\psi)\frac{\partial h}{\partial x}\right] + \frac{\partial}{\partial y}\left[K(\psi)\frac{\partial h}{\partial y}\right] + \frac{\partial}{\partial z}\left[K(\psi)\frac{\partial h}{\partial z}\right] = \frac{\partial \theta}{\partial t}$
(2.79)

It is usual to put Eq. (2.79) into a form where the independent variable is either θ or ψ. For the latter case it is necessary to multiply the top and bottom of the right-hand side by ∂ψ. Then, recalling the definition of the specific moisture capacity C [Eq. (2.65)], and noting that h = ψ + z, we obtain

$\frac{\partial}{\partial x}\left[K(\psi)\frac{\partial \psi}{\partial x}\right] + \frac{\partial}{\partial y}\left[K(\psi)\frac{\partial \psi}{\partial y}\right] + \frac{\partial}{\partial z}\left[K(\psi)(\frac{\partial \psi}{\partial z}+1)\right] = C(\psi)\frac{\partial \psi}{\partial t}$
(2.80)

Equation (2.80) is the ψ-based equation of flow for transient flow through an unsaturated porous medium. It is often called Richards equation, in honor of the soil physicist who first developed it (Richards, 1931). The solution ψ(x, y, z, t) describes the pressure head held at any point in a flow field at any time. It can easily be converted into an hydraulic head solution h(x, y, z, t) through the relation h = ψ + z. Solution requires knowledge of the characteristic curves K(ψ) and C(ψ) or θ(ψ).

The coupling of the unsaturated flow equation [Eq. (2.80)] and the saturated flow equation [Eq. (2.74)] has been attempted by Freeze (1971a) and by Narasimhan (1975). Improvements in the theory underlying saturated-unsaturated systems must await a better understanding of the principle of effective stress in the unsaturated zone.

Boundary-Value Problems

A boundary-value problem is a mathematical model. The technique of analysis inferred by this latter term is a four-step process, involving (1) examination of the physical problem, (2) replacement of the physical problem by an equivalent mathematical problem, (3) solution of the mathematical problem with the accepted techniques of mathematics, and (4) interpretation of the mathematical results in terms of the physical problem. Mathematical models based on the physics of flow usually take the form of boundary-value problems of the type pioneered by the developers of potential field theory and as applied in physics to such problems as the conduction of heat in solids (Carslaw and Jaeger, 1959).

To fully define a transient boundary-value problem for subsurface flow, one needs to know (1) the size and shape of the region of flow, (2) the equation of how within the region, (3) the boundary conditions around the boundaries of the region, (4) the initial conditions in the region, (5) the spatial distribution of the hydrogeologic parameters that control the flow, and (6) a mathematical method of solution. If the boundary-value problem is for a steady-state system, requirement (4) is removed.

Consider the simple groundwater flow problem illustrated in Figure 2.25(a). The region ABCD contains a homogeneous, isotropic porous medium of hydraulic conductivity K₁. The boundaries AB and CD are impermeable; the hydraulic heads on AD and BC are h₀ and h₁ respectively. Assuming steady flow and setting h₀ = 100 m and h₁ = 0 m, we can see by inspection that the hydraulic head at point E will be 50 m. Apparently, we made implicit use of properties (1), (3), and (5) from the list above; our method of solution (6) was one of inspection. It is not clear that we needed to know the equation of flow within the region. If we move to a more difficult problem such as that shown in Figure 2.25(b) (an earthfill dam resting on a sloping base), the value of the hydraulic head at point F does not come so easily. Here we would have to invoke a mathematical method of solution, and it would require that we know the equation of flow.

**Figure 2.25** Two steady-state boundary-value problems in the xy plane.

The methods of solution can be categorized roughly into five approaches: (1) solution by inspection, (2) solution by graphical techniques, (3) solution by analog model, (4) solution by analytical mathematical techniques, and (5) solution by numerical mathematical techniques. We have just seen an example of solution by inspection. The methods of flow-net construction presented in Chapter 5 can be viewed as graphical solutions to boundary-value problems. Electrical analog models are discussed in Sections 5.2 and 8.9. Numerical solutions are the basis of modern computer simulation techniques as described in Sections 5.3 and 8.8.

The most straightforward approach to the solution of boundary-value problems is that of analytical solutions. Many of the standard groundwater techniques presented later in the text are based on analytical solutions, so it is pertinent to examine a simple example. Consider, once again, the boundary-value problem of Figure 2.25(a). The analytical solution is

$h(x, y) = h_0 - (h_0 - h_1) \frac{x}{x_L}$ (2.81)

This is the equation of a set of equipotential lines traversing the field ABCD parallel to the boundaries AD and BC. Since the equipotentials are parallel to the y axis, his not a function of y and y does not appear on the right-hand side of Eq. (2.81). At point E, x/x_L = 0.5, and if h₀ = 100 m and h₁ = 0 m as before, then h_E from Eq. (2.81) is 50 m, as expected. In Appendix III, the separation-of-variables technique is used to obtain the analytical solution Eq. (2.81), and it is shown that the solution satisfies the equation of flow and the boundary conditions.

2.12 Limitations of the Darcian Approach

Darcy’s law provides an accurate description of the flow of groundwater in almost all hydrogeological environments. In general, Darcy’s law holds (1) for saturated flow and for unsaturated how, (2) for steady-state flow and for transient flow, (3) for flow in aquifers and for how in aquitards, (4) for flow in homogeneous systems and for how in heterogeneous systems, (5) for flow in isotropic media and for flow in anisotropic media, and (6) for how in both rocks and granular media. In this text, we will assume that Darcy’s law is a valid basis for our quantitative analyses.

Despite this soothing statement, or perhaps because of it, it is important that we examine the theoretical and practical limitations of the Darcian approach. It is necessary to look at the assumptions that underlie our definition of a continuum; examine the concepts of microscopic and macroscopic flow; investigate the upper and lower limits of Darcy’s law; and consider the particular problems associated with how in fractured rock.

Darcian Continuum and Representative Elementary Volume

In Section 2.1, it was noted that the definition of Darcy’s law requires the replacement of the actual ensemble of grains that make up a porous medium by a representative continuum. It was further stated that this continuum approach is carried out at a macroscopic rather than a microscopic scale. If Darcy’s law is a macroscopic law, there must be a lower limit to the size of an element of porous media for which the law is valid. Hubbert (1940) has addressed this problem. He defined the term macroscopic with the aid of Figure 2.26. This diagram is a hypothetical plot of the porosity of a porous medium as it might be measured on samples of increasing volume V₁, V₂, . . . , taken at a point P within a porous medium.

**Figure 2.26** Microscopic and macroscopic domains and the representative elementary volume V₃ (after Hubbert, 1956; Bear, 1972).

Bear (1972) defines the volume V₃ in Figure 2.26 as the representative elementary volume. He notes that it is a volume that must be larger than a single pore. In fact, it must include a sufficient number of pores to permit the meaningful statistical average required in the continuum approach. Below this volume there is no single value that can represent the porosity at P. Throughout this text the values of porosity, hydraulic conductivity, and compressibility refer to measurements that could be carried out on a sample larger than the representative elementary volume. In a more practical sense, they refer to values that can be measured on the usual sizes of cored soil samples. Where the scale of analysis involves volumes, such as V₅ in Figure 2.26, that may encompass more than one stratum in heterogeneous media, the scale is sometimes called megascopic.

The development of each of the equations of flow presented in Section 2.11 included the invocation of Darcy’s law. It must be recognized, then, that the methods of analysis that are based on boundary-value problems involving these equations apply on a macroscopic scale, at the level of the Darcian continuum. There are some groundwater phenomena, such as the movement of a tracer through a porous medium, that cannot be analyzed on this scale. It is therefore necessary to examine the interrelationship that exists between the Darcy velocity (or specific discharge) defined for the macroscopic Darcian continuum and {the microscopic velocities that exist in the liquid phase of the porous medium.

Specific Discharge, Macroscopic Velocity, and Microscopic Velocity

Our development will be more rigorous if we first differentiate, as Bear (1972) has done, between the volumetric porosity, n, which was defined in Section 2.5, and the areal porosity, n_A, which can be defined for any areal cross section through a unit volume, as n_A = A_v/A_T, where A_v is the area occupied by the voids and A_T is the total area. As suggested by Figure 2.27(a), various cross sections within a given unit volume may exhibit differing areal porosities n_A₁, n_A₂, . . . The volumetric porosity, n, is an average of the various possible areal porosities, n_Ai.

**Figure 2.27** Concepts of (a) areal porosity and (b) average linear flow velocity.

For any cross section A, the specific discharge, v, is defined from Eq. (2.1) as

$v = \frac{Q}{A}$

In that the volumetric flux Q is divided by the full cross-sectional area (voids and solids alike), this velocity is identified as being pertinent to the macroscopic continuum approach. In actual fact, the flow passes through only that portion of the cross-sectional area occupied by voids. For cross section A₁ we can define a velocity $\bar{v}_1 = Q/n_{A1}A$ that represents the volumetric flux divided by the actual cross-sectional area through which flow occurs. For the various sections A₁, A₂, . . . we can define $\bar{v}_1$ , $\bar{v}_2$ , . . . . If we denote their average by $\bar{v}$ then

$\bar{v} = \frac{Q}{nA} = \frac{v}{n} = -\frac{-K}{n}\frac{\partial h}{\partial l}$
(2.82)

The velocity $\bar{v}$ is known under a variety of names. We will refer to it as the average linear velocity. In that Q, n, and A are measurable macroscopic terms, so is $\bar{v}$ . It should be emphasized that $\bar{v}$ does not represent the average velocity of the water particles traveling through the pore spaces. These true, microscopic velocities are generally larger than $\bar{v}$ , because the water particles must travel along irregular paths that are longer than the linearized path represented by $\bar{v}$ . This is shown schematically in Figure 2.27(b). The true, microscopic velocities that exist in the pore channels are seldom of interest, which is indeed fortunate, for they are largely indeterminate. For all the situations that will be considered in this text, the Darcy velocity v and the average linear velocity $\bar{v}$ will suffice.

As a basis for further explanation of $\bar{v}$ , consider an experiment where a tracer is used to determine how much time is required for the bulk mass of groundwater to move a short but significant distance AB along a flow path. $\bar{v}$ is then defined as the ratio of travel distance to travel time, where the travel distance is defined as the linear distance from A to B and the travel time is the time required for the tracer to travel from A to B. In light of this conceptualization of $\bar{v}$ , Nelson (1968) has suggested a slightly different form of Eq. (2.82):

$\bar{v} = \frac{Q}{\epsilon nA} = \frac{v}{\epsilon n}$ (2.83)

where ε is an empirical constant dependent on the characteristics of the porous medium. Data obtained in laboratory experiments by Ellis et al. (1968) using relatively uniform sands indicate values of ε in the range 0.98–1.18. Values of ε for nonuniform sands and for other materials do not exist at present. In studies of groundwater tracers and groundwater contamination the almost universal unstated assumption is that ε = 1. For granular media this probably introduces little error. In fractured media the assumption may have less validity.

Upper and Lower Limits of Darcy’s Law

Even if we limit ourselves to the consideration of specific discharge on a macroscopic scale through the Darcian continuum, there may be limitations on the applicability of Darcy’s law. Darcy’s law is a linear law. If it were universally valid, a plot of the specific discharge v versus the hydraulic gradient dh/dl would reveal a straight-line relationship for all gradients between 0 and $\infty$ . For flow through granular materials there are at least two situations where the validity of this linear relationship is in question. The first concerns flow through low-permeability sediments under very low gradients and the second concerns large flows through very high permeability sediments. In other words, there may be both a lower limit and an upper limit to the range of validity of Darcy’s law. It has been suggested that a more general form of the porous media flow law might be

$v = -K\left(\frac{dh}{dl}\right)^m$ (2.84)

If m = 1, as it does in all the common situations, the flow law is linear and is called Darcy’s law; if m ≠ 1 the flow law is not linear and should not be called Darcy’s law.

For fine-grained materials of low permeability, it has been suggested on the basis of laboratory evidence that there may be a threshold hydraulic gradient below which flow does not take place. Swartzendruber (1962) and Bolt and Groenevelt (1969) review the evidence and summarize the various hypotheses that have been put forward to explain the phenomenon. As yet, there is no agreement on the mechanism, and the experimental evidence is still open to some doubt. In any event, the phenomenon is of very little practical importance; at the gradients being considered as possible threshold gradients, flow rates will be exceedingly small in any case.

Of greater practical importance is the upper limit on the range of validity of Darcy’s law. It has been recognized and accepted for many years (Rose, 1945; Hubbert, 1956) that at very high rates of flow, Darcy’s law breaks down. The evidence is reviewed in detail by both Todd (1959) and Bear (1972). The upper limit is usually identified with the aid of the Reynolds number R_e, a dimensionless number that expresses the ratio of inertial to viscous forces during flow. It is widely used in fluid mechanics to distinguish between laminar flow at low velocities and turbulent flow at high velocities. The Reynolds number for flow through porous media is defined as

$R_e = \frac{\rho vd}{\mu}$ (2.85)

where ρ and μ are the fluid density and viscosity, v the specific discharge, and d a representative length dimension for the porous medium, variously taken as a mean pore dimension, a mean grain diameter, or some function of the square root of the permeability k. Bear (1972) summarizes the experimental evidence with the statement that “Darcy’s law is valid as long as the Reynolds number based on average grain diameter does not exceed some value between 1 and 10” (p. 126). For this range of Reynolds numbers, all flow through granular media is laminar.

Flow rates that exceed the upper limit of Darcy’s law are common in such important rock formations as karstic limestones and dolomites, and cavernous volcanics. Darcian flow rates are almost never exceeded in nonindurated rocks and granular materials. Fractured rocks (and we will use this term to refer to rocks rendered more permeable by joints, fissures, cracks, or partings of any genetic origin) constitute a special case that deserves separate attention.

Flow in Fractured Rocks

The analysis of flow in fractured rocks can be carried out either with the continuum approach that has been emphasized thus far in this text or with a noncontinuum approach based on the hydraulics of flow in individual fractures. As with granular porous media, the continuum approach involves the replacement of the fractured media by a representative continuum in which spatially defined values of hydraulic conductivity, porosity, and compressibility can be assigned. This approach is valid as long as the fracture spacing is sufficiently dense that the fractured media acts in a hydraulically similar fashion to granular porous media. The conceptualization is the same, although the representative elementary volume is considerably larger for fractured media than for granular media. If the fracture spacings are irregular in a given direction, the media will exhibit trending heterogeneity. If the fracture spacings are different in one direction than they are in another, the media will exhibit anisotropy. Snow (1968, 1969) has shown that many fracture-flow problems can be solved using standard porous-media techniques utilizing Darcy’s law and an anisotropic conductivity tensor.

If the fracture density is extremely low, it may be necessary to analyze flow in individual fissures. This approach has been used in geotechnical applications where rock-mechanics analyses indicate that slopes or openings in rock may fail on the basis of fluid pressures that build up on individual critical fractures. The methods of analysis are based on the usual fluid mechanics principles embodied in the Navier-Stokes equations. These methods will not be discussed here. Wittke (1973) provides an introductory review.

Even if we limit ourselves to the continuum approach there are two further problems that must be addressed in the analysis of flow through fractured rock. The first is the question of non-Darcy flow in rock fractures of wide aperture. Sharp and Maini (1972) present laboratory data that support a nonlinear flow law for fractured rock. Wittke (1973) suggests that separate flow laws be specified for the linear-laminar range (Darcy range), a nonlinear laminar range, and a turbulent range. Figure 2.28 puts these concepts into the context of a schematic curve of specific discharge vs. hydraulic gradient. In wide rock fractures, the specific discharges and Reynolds numbers are high, the hydraulic gradients are usually less than 1, and the exponent m in Eq. (2.84) is greater than 1. These conditions lead to a downward deflection in the curve in Figure 2.28.

**Figure 2.28** Range of validity of Darcy’s law.

The second problem concerns the interaction of the three-dimensional stress field and the three-dimensional fluid flow field in rock. The general theoretical requirement for the coupling of these two fields was briefly discussed in Section 2.11, and reference was made there to the classic work of Biot (1941, 1955) for flow through porous media. For fractured rock, however, there is a further complication. Because the porosity of fractured rock is so low, the expansions and contractions of the fracture apertures that occur under the influence of changes in stress affect the values of hydraulic conductivity, K. The interaction between the fluid pressure p(x, y, z, t), or the hydraulic head h(x, y, z, t), and the effective stress σ_e(x, y, z, t) is thus complicated by the fact that K must be represented by a function, K(σ_e). The analysis of such systems, and the experimental determination of the nature of the K(σ_e) function, is a continuing subject of research in the fields of rock mechanics and groundwater hydrology.

Many researchers involved in the application of groundwater theory in rock mechanics have proposed formulas that relate the fracture porosity n_f and the hydraulic conductivity K of jointed rocks to the joint geometry. Snow (1968) notes that for a parallel array of planar joints of aperture, b, with N joints per unit distance across the rock face, n_f = Nb, and

$K = \left(\frac{\rho g}{\mu}\right) \left(\frac{Nb^3}{12}\right)$ (2.86)

$k = \frac{Nb^3}{12}$ (2.87)

where k is the permeability of the rock. N and b have dimensions 1/L and L, respectively, so that k comes out in units of L², as it should. Equation (2.86) is based on the hydrodynamics of flow in a set of planar joints. It holds in the linear-laminar range where Darcy’s law is valid. It must be applied to a block of rock of sufficient size that the block acts as a Darcian continuum. A permeability k, calculated with Eq. (2.87), can be considered as the permeability of an equivalent porous medium; one that acts hydraulically like the fractured rock.

Snow (1968) states that a cubic system of like fractures creates an isotropic system with a porosity n_f = 3Nb and a permeability twice the permeability that any one of its sets would contribute; that is, k = Nb³/6. Snow (1969) also provides predictive interrelationships between porosity and the anisotropic permeability tensor for three-dimensional joint geometries in which fracture spacings or apertures differ with direction. Sharp and Maini (1972) provide further discussion of the hydraulic properties of anisotropic jointed rock.

2.13 Hydrodynamic Dispersion

It is becoming increasingly common in the investigation of groundwater flow systems to view the flow regime in terms of its ability to transport dissolved substances known as solutes. These solutes may be natural constituents, artificial tracers, or contaminants. The process by which solutes are transported by the bulk motion of the flowing groundwater is known as advection. Owing to advection, nonreactive solutes are carried at an average rate equal to the average linear velocity, $\bar{v}$ , of the water. There is a tendency, however, for the solute to spread out from the path that it would be expected to follow according to the advective hydraulics of the flow system. This spreading phenomenon is called hydrodynamic dispersion. It causes dilution of the solute. It occurs because of mechanical mixing during fluid advection and because of molecular diffusion due to the thermal-kinetic energy of the solute particles. Diffusion, which is a dispersion process of importance only at low velocities, is described in Section 3.4. Emphasis in the present discussion is on dispersion that is caused entirely by the motion of the fluid. This is known as mechanical dispersion (or hydraulic dispersion). Figure 2.29 shows a schematic example of the results of this dispersive process in a homogeneous, granular medium.

Mechanical dispersion is most easily viewed as a microscopic process. On the microscopic scale, dispersion is caused by three mechanisms (Figure 2.30). The first occurs in individual pore channels because molecules travel at different velocities at different points across the channel due to the drag exerted on the fluid by the roughness of the pore surfaces. The second process is caused by the difference in pore sizes along the flow paths followed by the water molecules. Because of differences in surface area and roughness relative to the volume of water in individual pore channels, different pore channels have different bulk fluid velocities.

**Figure 2.29** Schematic representation of the dilution process caused by mechanical dispersion in granular porous media.

**Figure 2.30** Processes of dispersion on a microscopic scale.

The third dispersive process is related to the tortuosity, branching, and interfingering of pore channels. The spreading of the solute in the direction of bulk flow is known as longitudinal dispersion. Spreading in directions perpendicular to the flow is called transverse dispersion. Longitudinal dispersion is normally much stronger than lateral dispersion.

Dispersion is a mixing process. Qualitatively, it has a similar effect to turbulence in surface-water regimes. For porous media, the concepts of average linear velocity and longitudinal dispersion are closely related. Longitudinal dispersion is the process whereby some of the water molecules and solute molecules travel more rapidly than the average linear velocity and some travel more slowly. The solute therefore spreads out in the direction of flow and declines in concentration.

When a tracer experiment is set up in the laboratory, the only dispersion that can be measured is that which is observable at the macroscopic scale. It is assumed that this macroscopic result has been produced by the microscopic processes described above. Some investigators believe that heterogeneities on the macroscopic scale can cause additional dispersion to that caused by the microscopic processes. The concept of macroscopic dispersion is still not well understood. Dispersive processes are pursued further in Chapter 9.