Appendices | Freeze and Cherry Groundwater Book

Appendices

Appendix I: Elements of Fluid Mechanics

The analysis of groundwater flow requires an understanding of the elements of fluid mechanics. Albertson and Simons (1964) provide a useful short review; Streeter (1962) and Vennard (1961) are standard texts. Our purpose here is simply to a standtroduce the basic fluid properties: mass density, weight density, compressibility, and viscosity; and to examine the concepts of fluid pressure and pressure head.

An examination of the principles of fluid mechanics must begin with a review of the mechanics of materials in general. Table A1.1 provides a list of the basic mechanical properties of matter, together with their dimensions and units in the SI metric system. The SI system has as its basic dimensions: mass, length, and time; with basic SI units: kilogram (kg), meter (m), and second (s). All other properties are measured in units that are derived from this basic set. Some of these properties are so widely encountered, and their basic dimensions so complex, that special SI names have been coined for their derived units. As noted on Table A1.1, force and weight are measured in newtons (N), pressure and stress in N/m² or pascals (Pa), and work and energy in joules (J).

Table A1.1 Definitions, Dimensions, and SI Units for Basic Mechanical Properties

					Dimension of unit
Property	Symbol	Definition	SI unit	SI symbol	Derived	Basic
Mass	M		kilogram	kg		kg
Length	l		meter	m		m
Time	t		second	s		s
Area	A	A = l²				m²
Volume	V	V = l³				m³
Velocity	v	v = l/t				m/s
Acceleration	a	a = l/t²				m/s²
Force	F	F = Ma	newton	N	N	kg·m/s²
Weight	w	W = Mg	newton	N	N	kg·m/s²
Pressure	p	P = F/A	pascal	Pa	N/m²	kg/m·s²
Work	W	W = Fl	joule	J	N·m	kg·m²/s²
Energy		Work done	joule	J	N·m	kg·m²/s²
Mass density	ρ	ρ = M/V				kg/m³
Weight density	γ	γ = w/V			N/m³	kg/m²·s²
Stress	σ, τ	Internal response to external p	pascal	Pa	N/m²	kg/m·s²
Strain	ϵ	ϵ = ΔV/V				Dimensionless
Young’s modulus	E	Hooke’s law			N/m²	kg/m·s²

Table A1.2 provides an SI analysis of certain fluid properties and groundwater terms that occur in this text. Each is more fully described in Chapter 2.

Much of the technology associated with groundwater resource development in North America is still married to the FPS (foot-pound-second) system of units.

Table A1.2 Definitions, Dimensions, and SI Units for fluid Properties and Groundwater Terms

					Dimensions of unit
Property	Symbol	Definition	SI unit	SI symbol	Derived	Basic
Volume	V	V = l³	liter (= m³ × 10^-3)	ℓ	ℓ	m³
Discharge	Q	Q = l³/t			ℓ/s	m³/s
Fluid pressure	p	p = F/A	pascal	Pa	N/m²	kg/m·s²
Head	h					m
Mass density	ρ	ρ = M/V				kg/m³
Dynamic viscosity	μ	Newton’s law	centipoise (= N·s/m² × 10^-6)	cP	cP, N·s/m²	kg/m·s
Kinematic viscosity	v	v = μ/ρ	centistoke (=m²/s × 10^-6)	cSt	cSt	m²/s
Compressibility	α, β	α = 1/E			m²/N
Hydraulic conductivity	K	Darcy’s law			cm/s
Permeability	k	k = Kμ/pg			cm²	m²
Porosity	n					Dimensionless
Specific storage	S_S	S_S = pg(α + nβ)				1/m
Storativity	S	S = S_Sb*				Dimensionless
Transmissivity	T	T = Kb*				m²/s

*b, thickness of confined aquifer (see Section 2.10).

Table A1.3 Conversion Factors FPS (foot-pound-second)
System of Units to SI Units

	Multiply	By	To obtain
Length	ft	3.048 × 10^–1	m
	ft	3.048 × 10	cm
	ft	3.048 × 10^–4	km
	mile	1.069 × 10³	m
	mile	1.069	km
Area	ft²	9.290 × 10^–2	m²
	mi²	2.590	km²
	acre	4.407 × 10³	m²
	acre	4.407 × 10^–3	km²
Volume	ft³	2.832 × 10^–2	m³
	U.S. gal	3.785 × 10^–3	m³
	U.K. gal	4.546 × 10^–3	m³
	ft³	2.832 × 10
	U.S. gal	3.785
	U.K. gal	4.546
Velocity	ft/s	3.048 × 10^–1	m/s
	ft/s	3.048 × 10	cm/s
	mi/h	4.470 × 10^–1	m/s
	mi/h	1.609	km/h
Acceleration	ft/s²	3.048 × 10^–1	m/s²
Mass	lb_m*	4.536 × 10^–1	kg
	slug*	1.459 × 10	kg
	ton	1.016 × 10³	kg
Force and weight	lb_f*	4.448	N
	poundal	1.383 × 10^–1	N
Pressure and stress	psi	6.895 × 10³	Pa or N/m²
	lb_f/ft²	4.788 × 10^–1	Pa
	poundal/ft²	1.488	Pa
	atm	1.013 × 10⁵	Pa
	in Hg	3.386 × 10³	Pa
	mb	1.000 × 10²	Pa
Work and energy	ft-lb_f	1.356	J
	ft-poundal	4.214 × 10^–2	J
	Btu	1.055 × 10^–3	J
	calorie	4.187	J
Mass density	lb/ft³	1.602 × 10	kg/m³
	slug/ft³	5.154 × 10²	kg/m³
Weight density	lb_f/ft³	1.571 × 10²	N/m³
Discharge	ft³/s	2.832 × 10^–2	m³/s
	ft³/s	2.832 × 10	ℓ/s
	U.S. gal/min	6.309 × 10^–5	m³/s
	U.K. gal/min	7.576 × 10^–5	m³/s
	U.S. gal/min	6.309 × 10^–2	ℓ/s
	U.K. gal/min	7.576 × 10^–2	ℓ/s
Hydraulic conductivity (see also Table 2.3)	ft/s	3.048 × 10^–1	m/s
	U.S. gal/day/ft²	4.720 × 10^–7	m/s
Transmissivity	ft²/s	9.290 × 10^–2	m²/s
	U.S. gal/day/ft	1.438 × 10^–7	m²/s

*A body whose mass is 1 lb mass (lb_m) has a weight of 1 lb force (lb_f). 1 lb_f is the force required to accelerate a body of 1 lb_m to an acceleration of g = 32.2 ft/s². A slug is the unit of mass which, when acted upon by a force of 1 lb_f, acquires an acceleration of 1 ft/s².

Table A1.3 provides a set of conversion factors for converting FPS units to SI units.

The mass density (or simply, density) ρ of a fluid is defined as its mass per unit volume (Table A1.1). The weight density (or specific weight, or unit weight) γ of a fluid is defined as its weight per unit volume. The two parameters are related by

$\gamma = \rho g$ (A1.1)

For water, ρ = 1.0 g/cm³ = 1000 kg/m³; γ = 9.8 × 10³ N/m³. In the FPS system, γ = 62.4 lb_f/ft³.

The specific gravity G of any material is the ratio of its density (or specific weight) to that of water. For water, G = 1.0; for most soils and rocks, G = 2.65.

The viscosity of a fluid is the property that allows fluids to resist relative motion and shear deformation during flow. The more viscous the fluid, the greater the shear stress at any given velocity gradient. According to Newton’s law of viscosity,

$\tau = \mu \frac{dv}{dy}$ (A1.2)

where τ is the shear stress, dv/dy the velocity gradient, and μ the viscosity, or dynamic viscosity. The kinematic viscosity v is given by

$v = \frac{\mu}{\rho}$ (A1.3)

where ρ is the fluid density.

The compressibility of a fluid reflects its stress-strain properties. Stress is the internal response of a material to an external pressure. For fluids, stress is imparted through the fluid pressure. Strain is a measure of the linear or volumetric deformation of a stressed material. For fluids, it takes the form of reduced volume (and increased density) under increasing fluid pressures. The compressibility of water β is fully discussed in Section 2.9. It is defined by Eq. (2.44).

The density, viscosity, and compressibility of water are functions of temperature and pressure (Dorsey, 1940; Weast, 1972). In general, their variation is not great, and for the range of pressures and temperatures that occur in most ground-water applications, it is common to consider them as constants. At 15.5°C, ρ = 1.0 g/cm³, μ = 1.124 cP, and β = 4.4 × 10^–2 m²/N.

The fluid pressure p at any point in a standing body of water is the force per unit area which acts at that point. Under hydrostatic conditions, the fluid pressure at a point reflects the weight of the column of water overlying a unit cross-sectional area around the point. It is possible to express pressure relative to absolute zero pressure, but more commonly it is expressed relative to atmospheric pressure. In the latter case it is called gage pressure, as this is the pressure reading that is obtained on gages zeroed to the atmosphere.

The pressure head ψ at a point in a fluid is the height that a column of water would attain in a manometer placed at that point. In a standing body of water, ψ is equal to the depth of the point of measurement below the surface. If p is expressed as a gage pressure, ψ is defined by the relationship

$p = \rho g\psi = \gamma\psi$ (A1.4)

In effect, the pressure head ψ is a measurement of the fluid pressure p.

Fluid pressures are also developed in groundwater as it flows through porous geological formations and soils. In Section 2.2, the elements of fluid mechanics presented in this appendix are applied in the development of groundwater flow theory.

Appendix II: Equation of Flow for Transient Flow Through Deforming Saturated Media

A rigorous development of the equation of flow for transient flow in saturated media must recognize the fact that transient changes in fluid pressure lead to deformations in the granular skeleton of a porous medium, and these deformations imply that the medium, as well as the water, is in motion. This realization creates the need for two refinements to the classical derivation presented by Jacob (1940) and followed in Section 2.11 of this text. First, as recognized by Biot (1955), it is necessary to cast Darcy’s law in terms of the relative velocity of fluid to grains. And second, as recognized by Cooper (1966), it is necessary to consider the conservation of mass for the medium as well as for the fluid in the elemental control volume. One can develop the continuity relationships in one of three ways: (1) by considering a deforming elemental volume in deforming coordinates, (2) by considering a deforming elemental volume in fixed coordinates, or (3) by considering a fixed elemental volume in fixed coordinates. Following Gambolati and Freeze (1973), we will use a fixed elemental volume in fixed coordinates. The approach requires the use of vector notation and the material derivative (total derivative, substantial derivative). If these concepts are not familiar, Aris (1962) and Wills (1958) provide introductory treatments. The development will be presented here for a homogeneous, isotropic medium with hydraulic conductivity K, porosity n, and vertical compressibility α. The same approach is easily adapted to heterogeneous and anisotropic media.

In vector notation, the three-dimensional form of Darcy’s law [Eq. (2.34)] is

$\vec{v} = -K\nabla h$ (A2.1)

where $\vec{v} = (v_x, v_y, v_z)$ is the relative velocity of fluid to grains, and ∇h = (∂h/∂x, ∂h/∂y, ∂h/∂z) is the hydraulic gradient. We can expand the vector $\vec{v}$ as

$\vec{v} = n(\vec{v}_w - \vec{v}_s)$ (A2.2)

where $\vec{v}_w$ is the fluid velocity and $\vec{v}_s$ the velocity of the deforming medium.

The equation of state for the water [Eq. (2.47)] is

$\rho = \rho_0 e^\beta^p$ (A2.3)

and that for the soil grains, which are incompressible, is

ρ_s = constant (A2.4)

The equation of continuity for the water is

$-\nabla \cdot [np\vec{v}_w] = \frac{\partial}{\partial t}[np]$ (A2.5)

and that for the soil is

$-\nabla \cdot [(1 - n)\rho_s\vec{v}_s] = \frac{\partial}{\partial t}[(1 - n)\rho_s]$ (A2.6)

In these equations $\nabla \cdot$ is the divergence operator:

$\nabla \cdot = \frac{\partial}{\partial x} + \frac{\partial}{\partial y} + \frac{\partial}{\partial z}$

Expanding Eq. (A2.5), We arrive at

$-\rho\nabla \cdot (n\vec{v}_w) - n\vec{v}_w \cdot \nabla\rho = n\frac{\partial\rho}{\partial t} + \rho\frac{\partial n}{\partial t}$ (A2.7)

If we cancel ρ_s from Eq. (A2.6) and rearrange that equation, we obtain an expression for ∂n/∂t. This can be substituted in Eq. (A2.7) together with an expression for $n\vec{v}_w$ obtained from Eq. (A2.2). After dividing through by ρ and rearranging, Eq. (A2.7) becomes

$-\nabla \cdot \vec{v} - (\frac{v}{\rho}) \cdot \nabla\rho = (\frac{n}{\rho})(\frac{\partial\rho}{\partial t} + \vec{v}_s \cdot \nabla\rho) + \nabla \cdot \vec{v}_s$ (A2.8)

If we use the material derivative, $D/Dt = \partial/\partial t + \vec{v}_w \cdot \nabla$ , Eq. (A2.8) can be written as

$-\nabla \cdot \vec{v} - (\frac{v}{\rho}) \cdot \partial\rho = \frac{n}{\rho}\frac{Dp}{Dt} + \nabla \cdot \vec{v}_s$ (A2.9)

The first term on the right-hand side of Eq. (A2.9) can be related to the compressibility of water β by the relation

$\frac{D\rho}{Dt} = \rho\beta \frac{Dp}{Dt}$ (A2.10)

The material derivative on the right-hand side of Eq. (A2.10) can be replaced by a partial derivative only if the following inequality is satisfied:

$\vec{v}_s \cdot \nabla p \ll \frac{\partial p}{\partial t}$ (A2.11)

Let us further assume, on the left-hand side of Eq. (A2.9), that

$(\frac{\vec{v}}{\rho}) \cdot \nabla\rho \ll \nabla \cdot \vec{v}$ (A2.12)

Then, substituting Eqs. (A2.1) and (A2.10) in Eq. (A2.9) yields

$\nabla \cdot (K\nabla h) = n\beta\frac{\partial p}{\partial t} + \nabla \cdot \vec{v}_s$ (A2.13)

In a three-dimensional stress field, the grain velocity vector $\vec{v}_s = (v_s_x, v_s_y, v_s_z)$ is related to the deformation (or soil displacement) vector $\vec{u} = (u_x, u_y, u_z)$ by

$\vec{v}_s = \frac{D\vec{u}}{Dt}$ (A2.14)

In a one-dimensional vertical stress field,

$v_s_x = v_s_y = u_x = u_y = 0$ (A2.15)

If the conditions of Eq. (A2.15) are satisfied, the final term of Eq. (A2.13) can be expanded (Cooper, 1966; Gambolati and Freeze, 1973; Gambolati, 1973a) as

$\nabla \cdot \vec{v}_s = \frac{\partial v_s_z}{\partial z} = \frac{\partial}{\partial z}(\frac{Du_z}{Dt}) = \frac{D}{Dt}(\frac{\partial u_z}{\partial z}) = \alpha\frac{Dp}{Dt}$ (A2.16)

where α is the vertical compressibility of the porous medium. The change of derivative around the central equality is valid for a position vector but not in general. The material derivative in the right-hand expression of Eq. (A2.16) can be replaced by the partial derivative if Eq. (A2.11) is satisfied. In that case, Eq. (A2.13) becomes

$\nabla \cdot (K\nabla h) = n\beta\frac{\partial p}{\partial t} + \alpha\frac{\partial p}{\partial t}$ (A2.17)

Since ρ = ρg(h – z) and K is a constant, Eq. (A2.17) simplifies to

$\nabla^2h = \frac{\rho g(\alpha +n\beta}{K}\frac{\partial h}{\partial t}$ (A2.18)

Or, recalling that S_s = ρ g(α + nβ) and expanding the vector notation,

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

Equation (A2.19) is identical to Eq. (2.75) developed by Jacob (1940). The more rigorous development makes it clear that the validity of the classical equation of flow rests on the satisfaction of the inequalities of Eqs. (A2.11) and (A2.12) and the stress condition of Eq. (A2.15). It is unlikely that these conditions are always satisfied. Gambolati (1973b) shows that where the rate of consolidation $\vec{v}_s$ exceeds the rate of percolating fluid $\vec{v}_w$ , as it can in thick clay layers of low permeability and high compressibility, the inequalities may not be satisfied. With regard to the stress condition, the $\nabla \cdot \vec{v}_s$ , term at the end of Eq. (A2.13) is really the protruding tip of an iceberg that relates the three-dimensional flow field to the three-dimensional stress field. Biot (1941, 1955) first exposed the interrelationships and Verruijt (1969) provides a clear derivation. Schifiman et al. (1969) provides a comparison of the classical and Biot approaches, and Gambolati (1974) analyses the range of validity of the classical equation of flow.

Appendix III: Example of an Analytical Solution to a Boundary-Value Problem

Consider the simple groundwater flow problem shown in Figure 2.25(a). The equation of flow for steady-state saturated flow in the xy plane is

$\frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial y^2} = 0$ (A3.1)

The mathematical statement of the boundary conditions is

$\frac{\partial h}{\partial y} = 0$ on $y = 0$ and $y = y_L$ (A3.2)

$h = h_0$ on $x = 0$ (A3.3)

$h = h_1$ on $x = x_L$ (A3.4)

We will solve for h(x, y) using the separation-of-variables technique.

Assume that the solution is a product solution of the form

$h(x, y) = X(x) \cdot Y(y)$ (A3.5)

Equation (A3.1) then becomes

$Y\frac{\partial^2x}{\partial x^2} + X\frac{\partial^2y}{\partial y^2} = 0$ (A3.6)

Dividing through by XY yields

$\frac{1}{X}\frac{\partial^2X}{\partial x^2} = \frac{1}{Y}\frac{\partial^2Y}{\partial y^2}$ (A3.7)

The left-hand side is independent of y. Therefore the right-hand side, despite its appearance, must also be independent of y, since it is identically equal to the left-hand side. Similarly, the right-hand side is independent of x, and so too is the left-hand side. If both sides are independent of x and y, each side must equal a constant. Therefore,

$\frac{1}{X}\frac{\partial^2X}{\partial x^2} = G$ and $\frac{1}{Y}\frac{\partial^2Y}{\partial y^2} = G$ (A3.8)

The constant G may be positive, negative, or zero. All three cases lead to product solutions, but only the case G = 0 leads to a solution that is physically meaningful for this problem. Equations (A3.8) become

$\frac{1}{X}\frac{\partial^2X}{\partial x^2} = 0$ and $\frac{1}{Y}\frac{\partial^2Y}{\partial y^2} = 0$ (A3.9)

These are ordinary differential equations whose solutions are well known:

$X = Ax + B$ and $Y = Cy + D$ (A3.10)

The product solution Eq. (A3.5) becomes

$h(x, y) = (Ax + B)(Cy +D)$ (A3.11)

We can evaluate the coefficients A, B, C and D by invoking the boundary conditions. Differentiating Eq. (A3.11) with respect to y yields

$\frac{\partial h}{\partial y} = (Ax + B)C$ (A3.12)

and invocation of Eq. (A3.2) implies that C = 0. From Eq. (A3.11), we are left with

$h(x, y) = (Ax + B)D = Ex + F$ (A3.13)

Invoking the boundary conditions of Eqs. (A3.3) and (A3.4) yields F = h₀ and E = –(h₀ – h₁)/x_L. The solution is, therefore,

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

This equation is identical to Eq. (2.81), presented in Section 2.11 without proof.

It is immediately clear that Eq. (A3.14) satisfies the boundary conditions of Eqs. (A3.3) and (A3.4). Differentiation with respect to y yields zero in satisfaction of Eq. (A3.2). Differentiation with respect to x twice also yields zero, so the solution Eq. (A3.14) satisfies the equation of flow [Eq. (A3.1)].

Appendix IV: Debye-Hückel Equation and Kielland Table For Ion-Activity Coefficients

Debye-Hückel expression for individual ion-activities:

$log$ $\gamma = \frac{-Az^2\sqrt{I}}{1 + \r{a}B\sqrt{I}}$

Values of the Ion-Size Parameter å for Common Ions Encountered in Natural Water:

å × 10⁸	Ion
2.5	$\ce{NH^+_4}$
3.0	K⁺, Cl^–, $\ce{NO^-_3}$
3.5	OH^–, HS^–, $\ce{MnO^-_4}$ , F^–
4.0	SO₄^2–, PO₄^3–, HPO₄^2–
4.0–4.5	Na⁺, $\ce{HCO^-_3}$ , $\ce{H2PO^-_4}$ , $\ce{HSO^-_3}$
4.5	CO₃^2–, SO₃^2–
5	Sr²⁺, Be²⁺, S^2–
6	Ca²⁺, Fe²⁺, Mn²⁺
8	Mg²⁺
9	H⁺, Al³⁺, Fe³⁺

Parameters A and B at 1 Bar

Temperature (°C)	A	B (× 10^–8)
0	0.4883	0.3241
5	0.4921	0.3249
10	0.4960	0.3258
15	0.5000	0.3262
20	0.5042	0.3273
25	0.5085	0.3281
30	0.5130	0.3290
35	0.5175	0.3297
40	0.5221	0.3305
50	0.5319	0.3321
60	0.5425	0.3338

Kielland Table of Ion-Activity Coefficients at 25º Celcius Arranged by the Size of Ions (Based on Debye-Hückel Equation)

Charge	*Size, å**	Ions	l = 0.0005	0.001	0.0025	0.005	0.01	0.025	0.05	0.1
1	2.5	Rb⁺, Cs⁺, Ag⁺, $\ce{NH^+_4}$ , Tl⁺	0.975	0.964	0.945	0.924	0.898	0.85	0.80	0.75
	3	K⁺, Cl^–, Br^–, I^–, CN^–, $\ce{NO^-_3}$ , $\ce{NO^-_2}$ , OH^–, F^–, $\ce{ClO^-_4}$	0.975	0.964	0.945	0.925	0.899	0.85	0.805	0.755
	4	Na⁺, $\ce{IO^-_3}$ , $\ce{HCO^-_3}$ , $\ce{HSO^-_3}$ , $\ce{H2PO^-_4}$ , $\ce{CLO^-_2}$ , $\ce{C2H3O^-_2}$	0.975	0.964	0.947	0.928	0.902	0.86	0.82	0.775
	6	Li⁺, $\ce{C6H3O^-_2}$	0.975	0.965	0.948	0.929	0.907	0.87	0.835	0.80
	9	H⁺	0.975	0.967	0.950	0.933	0.914	0.88	0.86	0.83
2	4.5	Pb²⁺, Hg₂²⁺, SO₄^2–, CrO₄^2–, CO₃^2–, SO₃^2–, C₂O₄^2–, S₂O₃^2–, H citrate^2–	0.903	0.867	0.805	0.742	0.665	0.55	0.455	0.37
	5	Sr²⁺, Ba²⁺, Cd²⁺, Hg²⁺, S^2–, WO₄^2–	0.903	0.868	0.805	0.744	0.67	0.555	0.465	0.38
	6	Ca²⁺, Cu²⁺, Zn²⁺, Sn²⁺, Mn²⁺, Fe²⁺, Ni²⁺, Co²⁺, phthalate^2–	0.905	0.870	0.809	0.749	0.675	0.57	0.485	0.405
	8	Mg²⁺, Be²⁺	0.906	0.872	0.813	0.755	0.69	0.595	0.52	0.45
3	4	PO₄^3–, Fe(CN)₆^3–, Cr(NH₃)₆³⁺	0.796	0.725	0.612	0.505	0.395	0.25	0.16	0.095
3	9	Al³⁺, Fe³⁺, Cr³⁺, Sc³⁺, In³⁺, and rare earths	0.802	0.738	0.632	0.54	0.445	0.325	0.245	0.18

*Note that these sizes are rounded values for the effective size in water solution and are not the size of the simple ions, unhydrated. For a more detailed discussion, see the original paper from which these values were taken [J. Kielland, J. Amer. Chem. Soc., 59, 1675 (1937)].

References

BERNER, R. A. 1971. Principles of Chemical Sedimentology. McGraw-Hill, New York, 240 pp.

KLOTZ, I. M. 1950. Chemical Thermodynamics. Prentice-Hall, Eaglewood Cliffs, N.J., 369 pp.

MANOV, G. G., R. G. BATES, W. J HAMER, and S. F. ACREE. 1943. Values of the constants in the Debye-Hückel equation for activity coefficients. J. Amer. Chem. Soc., 65, pp. 1765–1767.

Appendix IV: Debye-Hückel Equation and Kielland Table For Ion-Activity Coefficients

Debye-Hückel expression for individual ion-activities:

$log$ $\gamma = \frac{-Az^2\sqrt{I}}{1 + \r{a}B\sqrt{I}}$

Values of the Ion-Size Parameter å for Common Ions Encountered in Natural Water:

å × 10⁸	Ion
2.5	$\ce{NH^+_4}$
3.0	K⁺, Cl^–, $\ce{NO^-_3}$
3.5	OH^–, HS^–, $\ce{MnO^-_4}$ , F^–
4.0	SO₄^2–, PO₄^3–, HPO₄^2–
4.0–4.5	Na⁺, $\ce{HCO^-_3}$ , $\ce{H2PO^-_4}$ , $\ce{HSO^-_3}$
4.5	CO₃^2–, SO₃^2–
5	Sr²⁺, Be²⁺, S^2–
6	Ca²⁺, Fe²⁺, Mn²⁺
8	Mg²⁺
9	H⁺, Al³⁺, Fe³⁺

Parameters A and B at 1 Bar

Temperature (°C)	A	B (× 10^–8)
0	0.4883	0.3241
5	0.4921	0.3249
10	0.4960	0.3258
15	0.5000	0.3262
20	0.5042	0.3273
25	0.5085	0.3281
30	0.5130	0.3290
35	0.5175	0.3297
40	0.5221	0.3305
50	0.5319	0.3321
60	0.5425	0.3338

Kielland Table of Ion-Activity Coefficients at 25º Celcius Arranged by the Size of Ions (Based on Debye-Hückel Equation)

Charge	*Size, å**	Ions	l = 0.0005	0.001	0.0025	0.005	0.01	0.025	0.05	0.1
1	2.5	Rb⁺, Cs⁺, Ag⁺, $\ce{NH^+_4}$ , Tl⁺	0.975	0.964	0.945	0.924	0.898	0.85	0.80	0.75
	3	K⁺, Cl^–, Br^–, I^–, CN^–, $\ce{NO^-_3}$ , $\ce{NO^-_2}$ , OH^–, F^–, $\ce{ClO^-_4}$	0.975	0.964	0.945	0.925	0.899	0.85	0.805	0.755
	4	Na⁺, $\ce{IO^-_3}$ , $\ce{HCO^-_3}$ , $\ce{HSO^-_3}$ , $\ce{H2PO^-_4}$ , $\ce{CLO^-_2}$ , $\ce{C2H3O^-_2}$	0.975	0.964	0.947	0.928	0.902	0.86	0.82	0.775
	6	Li⁺, $\ce{C6H3O^-_2}$	0.975	0.965	0.948	0.929	0.907	0.87	0.835	0.80
	9	H⁺	0.975	0.967	0.950	0.933	0.914	0.88	0.86	0.83
2	4.5	Pb²⁺, Hg₂²⁺, SO₄^2–, CrO₄^2–, CO₃^2–, SO₃^2–, C₂O₄^2–, S₂O₃^2–, H citrate^2–	0.903	0.867	0.805	0.742	0.665	0.55	0.455	0.37
	5	Sr²⁺, Ba²⁺, Cd²⁺, Hg²⁺, S^2–, WO₄^2–	0.903	0.868	0.805	0.744	0.67	0.555	0.465	0.38
	6	Ca²⁺, Cu²⁺, Zn²⁺, Sn²⁺, Mn²⁺, Fe²⁺, Ni²⁺, Co²⁺, phthalate^2–	0.905	0.870	0.809	0.749	0.675	0.57	0.485	0.405
	8	Mg²⁺, Be²⁺	0.906	0.872	0.813	0.755	0.69	0.595	0.52	0.45
3	4	PO₄^3–, PO₄^3–	0.796	0.725	0.612	0.505	0.395	0.25	0.16	0.095
3	9	Al³⁺, Fe³⁺, Cr³⁺, Sc³⁺, In³⁺, and rare earths	0.802	0.738	0.632	0.54	0.445	0.325	0.245	0.18

References

BERNER, R. A. 1971. Principles of Chemical Sedimentology. McGraw-Hill, New York, 240 pp.

KLOTZ, I. M. 1950. Chemical Thermodynamics. Prentice-Hall, Eaglewood Cliffs, N.J., 369 pp.

MANOV, G. G., R. G. BATES, W. J HAMER, and S. F. ACREE. 1943. Values of the constants in the Debye-Hückel equation for activity coefficients. J. Amer. Chem. Soc., 65, pp. 1765–1767.

Appendix V: Complementary Error Function (erfc)

$\text{erf}$ $(\beta) = \frac{2}{\sqrt{\mu}}\int_{0}^{\beta}e^{-e^2}d \epsilon$

$\text{erf}$ $(-\beta)$ = $\text{-erf}\beta$

$\text{erfc}$ $(-\beta) = 1 - \text{erf}$ $(\beta)$

$\pmb{\beta}$	erf ( $\pmb{\beta}$ )	erfc ( $\pmb{\beta}$ )
0	0	1.0
0.05	0.056372	0.943628
0.1	0.112463	0.887537
0.15	0.167996	0.832004
0.2	0.222703	0.777297
0.25	0.276326	0.723674
0.3	0.328627	0.671373
0.35	0.379382	0.620618
0.4	0.428392	0.571608
0.45	0.475482	0.524518
0.5	0.520500	0.479500
0.55	0.563323	0.436677
0.6	0.603856	0.396144
0.65	0.642029	0.357971
0.7	0.677801	0.322199
0.75	0.711156	0.288844
0.8	0.742101	0.257899
0.85	0.770668	0.229332
0.9	0.796908	0.203092
0.95	0.820891	0.179109
1.0	0.842701	0.157299
1.1	0.880205	0.119795
1.2	0.910314	0.089686
1.3	0.934008	0.065992
1.4	0.952285	0.047715
1.5	0.966105	0.033895
1.6	0.976348	0.023652
1.7	0.983790	0.016210
1.8	0.989091	0.010909
1.9	0.992790	0.007210
2	0.995322	0.004678
2.1	0.997021	0.002979
2.2	0.998137	0.001863
2.3	0.998857	0.001143
2.4	0.999311	0.000689
2.5	0.999593	0.000407
2.6	0.999764	0.000236
2.7	0.999866	0.000134
2.8	0.999925	0.000075
2.9	0.999959	0.000041
3	0.999978	0.000022

Appendix VI: Development of Finite-Difference Equation for Steady-State Flow in a Homogeneous, Isotropic Medium

The partial differential equation that describes steady-state flow in a two-dimensional, homogeneous, isotropic region of flow (Section 2.11) is Laplace’s equation:

$\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2}$ (A6.1)

To find the finite-difference equation for an interior node in the nodal grid used to discretize the region of flow, we must replace the second-order partial derivatives in Eq. (A6.1) by differences. Let us consider the first term of the equation first. Recall that the definition of the partial derivative with respect to x of a function of two variables h(x, z) is

$\frac{\partial h}{\partial x} = \substack {lim \\ \Delta x \rightarrow 0} \frac{h(x + \Delta x, z) - h(x, z)}{\Delta x}$ (A6.2)

On a digital computer it is impossible to take the limit as $\Delta x \rightarrow 0$ , but it is possible to approximate the limit by assigning to Δx some arbitrarily small value; in fact, we can do so by designing a nodal network with a mesh spacing of Δx.

For any value of z, say z₀, we can expand h(x, z₀) in a Taylor’s expansion about the point (x₀, z₀) as follows:

$h(x, z_0) = h(x_0, z_0) + (x - x_0)\frac{\partial h}{\partial x}(x_0, z_0) + \frac{(x - x_0)^2}{2}\frac{\partial^2h}{\partial x^2}(x_0, z_0) + ...$ (A6.3)

If we let x = x₀ + Δx (this is known as a forward difference) and abandon all the terms of order greater than unity, we can approximate ∂h/∂x by

$\frac{\partial h}{\partial x}(x_0, z_0) = \frac{h(x_0 + \Delta x, z_0) - h(x_0, z_0)}{\Delta x}$ (A6.4)

The abandoned terms of the Taylor expansion represent the truncation error in the finite-difference approximation.

We can obtain a similar expression to Eq. (A6.4) by substituting the backward difference, x = x₀ – Δx, into Eq. (A6.3). This yields

$\frac{\partial h}{\partial x}(x_0, z_0) = \frac{h(x_0, z_0) - h(x_0 - \Delta x, z_0)}{\Delta x}$ (A6.5)

To obtain the approximation for ∂₂h/∂x₂, we write the difference equation in terms of ∂h/∂x using a forward-difference expression:

$\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{\frac{\partial h}{\partial x}(x_0 + \Delta x, z_0) - \frac{\partial h}{\partial x}(x_0, z_0)}{\Delta x}$ (A6.6)

and substitute the backward-difference expression Eq. (A6.5) in Eq (A6.6) to get

$\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{h(x_0 + \Delta x, z_0) - 2h(x_0, z_0) + h(x_0 + \Delta x, z_0)}{(\Delta x)^2}$ (A6.7)

In a similar manner, we can develop a difference expression for ∂₂h/∂x₂ as

$\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{h(x_0, z_0 + \Delta z) - 2h(x_0, z_0) + h(x_0, z_0 - \Delta z)}{(\Delta z)^2}$ (A6.8)

For a square grid, Δx = Δz; adding Eqs. (A6.7) and (A6.8) to for Laplace’s equation yields

$\frac{1}{(\Delta z)^2}[h(x_0 + \Delta x, z_0) + h(x_0 - \Delta x, z_0) + h(x_0, z_0 + \Delta z) + h(x_0, z_0 - \Delta z) + 4h(x_0, z_0) = 0]$ (A6.9)

If we let (x₀, z₀) be the nodal point (i, j), Eq. (A6.9) can be rearranged to yield

$h_{ij} = \frac{1}{4}(h_{i+1,j} + h_{i-1,j} + h_{i,j+1} + h_{i,j-1})$ (A6.10)

which is identical to Eq. (5.24).

Appendix VII: Tóth’s Analytical Solution for Regional Groundwater Flow

Tóth (1962, 1963) has presented two analytical solutions for the boundary-value problem representing steady-state flow in a vertical, two-dimensional, saturated, homogeneous, isotropic flow field bounded on top by a water table and on the other three sides by impermeable boundaries (the shaded cell of Figure 6.1, reproduced in an xz coordinate system in Figure A7.1).

**Figure A7.1** Region of flow for Tóth’s analytical solution.

He first considered the case where the water-table configuration is a straight line of constant slope. For this case, the region of flow in Figure A7.1 is the region ABCEA. Since it is not possible to obtain an analytical solution in a trapezoidal region, Tóth approximated the actual region of flow by the shaded region ABCDA. He projected the hydraulic-head values that exist along the actual water table AE onto the upper boundary AD of the region of solution. The approximation is satisfactory for small α.

The equation of flow is Laplace’s equation:

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

For a region bounded by x = s and z = z₀, the boundary conditions are

$\frac{\partial h}{\partial x}(0, z) = \frac{\partial h}{\partial x}(s, z) = 0$ on AB and CD

$\frac{\partial h}{\partial x}(x, 0) = 0$ on BC (7.2)

$h(x, z_0) = z_0 + cx$ on AD

where c = tan α.

The analytical solution, obtained by separation of variables, is

$h(x, z) = z_0 + \frac{cs}{2} - \frac{4cs}{\pi^2}\sum_{m=0}^{\infty} \frac{\cos[(2m + 1)\pi x/s] \cos h[(2m + 1)\pi z/s]}{(2m +1)^2 \cos h[(2m +1)\pi z_0/s]}$ (7.3)

This equation satisfies both the equation of flow (A7.1) and the boundary conditions (A7.2). When plotted and contoured, it leads to the equipotential net shown in Figure A7.1. Flow is from the recharge boundary DF through the field to the discharge boundary AF.

Tóth also considered the case where the water-table configuration is specified as a sine curve superimposed on the line AE (to represent hummocky topography). The final boundary condition then becomes

$h(x, z_0) = z_0 + cx + \alpha \sin bx$ on AD (A7.4)

where $c = \tan \alpha$ , $a = a'/\cos \alpha$ , and $b = b'/\cos \alpha$ , a’ being the amplitude of the sine curve and b’ being the frequency.

The analytical solution for this case takes the form

$h(x, z) = z_0 + \frac{cs}{2} + \frac{\alpha}{sb}(1 - \cos bs) + 2 \displaystyle\sum_{m=1}^{\infty} \left[\frac{ab(1 - \cos bs \cos m\pi)}{b^2 - \frac{m^2\pi^2}{s^2}} + \frac{cs^2}{m^2\pi^2}(\cos m\pi - 1)\right] \times \frac{\cos (m\pi x/s) \cos h (m\pi z/s)}{s \cdot \cos h (m\pi z_0/s)}$ (A7.5)

The equipotential net described by this function is similar to that of Figure 6.2(b).

Appendix VIII: Numerical Solution of the Boundary-Value Problem Representing One-Dimensional Infiltration Above a Recharging Groundwater Flow System

Consider a one-dimensional, vertical, homogeneous flow field with its upper boundary at the ground surface and its lower boundary below the water table. As outlined in Section 6.4, the equation of flow for such a system is

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

and the boundary conditions are

$\frac{\partial\psi}{\partial t} = \frac{R}{K(\psi)} - 1$ (A8.2)

at the top, and

$\frac{\partial\psi}{\partial z} = \frac{Q}{K_0} - 1$ (A8.3)

at the base. The solution is in terms of the pressure head, $\psi(z, t)$ . The position of the water table at any time is given by the value of z at which $\psi = 0$ . The required input includes the initial conditions, $\psi(z, 0)$ , the rainfall rate R, the groundwater recharge rate Q, and the unsaturated characteristic curves K(ψ) and C(ψ). For $\psi \geq \psi_a$ , K = K₀ and C = 0, $\psi_a$ being the air entry pressure head.

The numerical finite-difference scheme used by Rubin and Steinhardt (1963), Liakopoulos (1965b), and Freeze (1969b) is the implicit scheme of Richtmyer (1957). With this method the (z, t) plane is represented by a rectangular grid of points j = 1, 2, …, L along the z axis, and n = 1, 2, … along the t axis. The distance between the vertical nodes is Δz, and between the time steps, Δt. The solution at any given grid point (j, n) is $\psi_j^n$ . Under this notation the finite-difference form for Eq. (A8.1) for an interior node (j = 2 to j = L – 1), after suitable rearrangement, is

$C(\psi_{\text{III}})\left(\frac{\psi_j^n - \psi_j^{n-1}}{\Delta t}\right) = \frac{1}{\Delta z}\left\{\left[K(\psi_{\text{I}})(1 + \frac{1}{2\Delta z})[\psi_{j+1}^{n-1} + \psi_{j+1}^{n} - \psi_{j}^{n-1} - \psi_j^n]\right]\right\}$ (A8.4)

where the values of $\psi_{\text{I}}$ , $\psi_{\text{II}}$ , and $\psi_{\text{III}}$ , which determine the K and C values to be applied for a particular node at a particular time step, are determined by extrapolation from previous time steps as described by Rubin and Steinhardt (1963).

Similar finite-difference equations, which incorporate the boundary conditions (A8.2) and (A8.3), can be written for the nodes at the upper and lower boundary (j = 1, j = L).

For each time step, the finite-difference approximations constitute a system of L linear algebraic equations in L unknowns. The general form is

$-A_j\psi_{j+1}^n + B_j\psi_j^n - C_j\psi_{j-1}^n = D_j$ (A8.5)

where the coefficients A, B, C, and D vary with the node (j = 1, j = 2 à L – 1, j = L), with the time step (n = 1, n = 2, n > 2), and with the saturation $(\psi \geq \psi_a, \psi < \psi_a)$ . The variables on which A, B, C, and D depend are the ψ values from the previous time step, the boundary values R and Q, and the functional relationships K(ψ) and C(ψ).

The values of $\psi_j^n$ are calculated from the recurrence relation

$\psi_j^n = E_j\psi_{j+1}^n + F_j$ (A8.6)

where

$E_j = \frac{A_j}{B_j - C_j E_{j-1}}$

$E_j = \frac{D_j + C_j F_{j-1}}{B_j - C_j E_{j-1}}$

The E and F coefficients are calculated from j = 1 to j = L, and the ψ’s are back-calculated from j = L to j = 1 using Eq. (A8.6). For j = 1, E₁ = A₁/B₁ and F₁ = D₁/B₁. Complications arise if the upper node j = L is saturated.

If K(ψ) and C(ψ) are hysteretic, they are usually built into the computer program in the form of a table of values representing the main wetting, main drying, and principal scanning curves. The program must locate the correct curve on the basis of whether the node in question is wetting or drying, and by scanning the past history of the node to determine the ψ value at which a change from wetting to drying or vice versa occurred.

Appendix IX: Development of Finite-Difference Equation for Transient Flow in a Heterogeneous, Anisotropic, Horizontal, Confined Aquifer

The partial differential equation that describes transient flow through a saturated anisotropic medium (Section 2.11) is

$\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) = S_s\frac{\partial h}{\partial t}$ (A9.1)

For a horizontal, confined aquifer of thickness, b, the two-dimensional form of Eq. (A9.1) reduces to

$\frac{\partial}{\partial x} \left(T_x\frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y} \left(T_y\frac{\partial h}{\partial y}\right) = S_s\frac{\partial h}{\partial t}$ (A9.2)

where T_x and T_y are the principal components of the transmissivity tensor defined by T = Kb, and S is the storativity defined by S = S_sb. To find the finite-difference equation for an interior node in the nodal grid used to discretize the region of flow, we must replace the partial derivatives in Eq. (A9.2) by differences. Using the definitions developed in Appendix VI and the ijk notation of Section 8.8 and Figure 8.26(c), we can write the finite-difference expression:

$\frac{\partial}{\partial x} \left(T_x\frac{\partial h}{\partial x}\right) \simeq \frac{1}{\Delta x}\left[\left(T_x\frac{\partial h}{\partial x}\right)_{i+1/2,j}^k - \left(T_x\frac{\partial h}{\partial x}\right)_{i-1/2,j}^k\right]$ (A9.3)

where the subscript $(i + \frac{1}{2}, j)$ means that the bracketed quantity is evaluated at the midpoint between nodes (i, j) and (i + 1, j), and the superscript k means that the bracketed quantity is evaluated at time step k. We can further approximate terms on the right-hand side of Eq. (A9.3) by

$\left(T_x\frac{\partial h}{\partial x}\right)_{i+1/2,j}^k = (T_x)_{i+1/2,j}(h_{i+1,j}^k - h_{i,j}^k)\frac{1}{\Delta x}$ (A9.4A)

$\left(T_x\frac{\partial h}{\partial x}\right)_{i-1/2,j}^k = (T_x)_{i-1/2,j}(h_{i,j}^k - h_{i-1,j}^k)\frac{1}{\Delta x}$ (A9.4B)

If we evaluate $(T_x)_{i+1/2,j}$ and $(T_x)_{i-1/2,j}$ by simple averages of the form

$(T_x)_{i+1/2,j} \simeq \frac{(T_x)_{i+1,j} + (T_x)_{i,j}}{2}$ (A9.5)

then these expressions can be substituted in Eqs. (A9.4), and Eqs. (A9.4) can in turn be substituted in Eq. (A9.3) to give

$\frac{\partial}{\partial x} \left(T_x\frac{\partial h}{\partial x}\right) \simeq \frac{1}{2\Delta x^2} \left\{\left[(T_x)_{i+1/2,j} + (T_x)_{i,j}\right]h_{i+1,j}^k - \left[(T_x)_{i+1/2,j} + 2(T_x)_{i,j} + (T_x)_{i-1,j}\right]h_{i,j}^k + \left[(T_x)_{i,j} + (T_x)_{i-1,j}\right]h_{i-1,j}^k \right\}$ (A9.6)

Similarly,

$\frac{\partial}{\partial y} \left(T_y\frac{\partial h}{\partial y}\right) \simeq \frac{1}{2\Delta y^2} \left\{\left[(T_y)_{i+1/2,j} + (T_y)_{i,j}\right]h_{i+1,j}^k - \left[(T_y)_{i+1/2,j} + 2(T_y)_{i,j} + (T_y)_{i-1,j}\right]h_{i,j}^k + \left[(T_y)_{i,j} + (T_y)_{i-1,j}\right]h_{i-1,j}^k \right\}$ (A9.7)

Finally, we can approximate the right-hand side of Eq. (A9.2) by

$S\frac{\partial h}{\partial t} \simeq S_{i,j}\left(\frac{h_{i,j}^k - h_{i,j}^{k-1}}{\Delta t}\right)$ (A9.8)

Substituting Eqs. (A9.6), (A9.7), and (A9.8) for the three terms in Eqs. (A9.6) and gathering terms leads to the general finite-difference equation for an internal node in a heterogeneous, anisotropic aquifer:

$Ah_{i,j}^k = Bh_{i,j-1}^k + Ch_{i+1,j}^k + Dh_{i,j+1}^k + Eh_{i-1,j}^k + F$ (A9.9)

where

$A = \frac{1}{2\Delta x}^2[(T_x)_{i+1,j} + 2(T_x)_{i,j} + (T_x)_{i-1,j}] + \frac{1}{2\Delta y^2}[(T_y)_{i,j+1} + 2(T_y)_{i,j} + (T_y)_{i,j-1}] + \frac{S_{i,j}}{\Delta t}$

$B = \frac{1}{2\Delta y^2}[(T_y)_{i,j} + (T_y)_{i,j-1}]$

$C = \frac{1}{2\Delta x^2}[(T_x)_{i+1,j} + (T_x)_{i,j}]$

$D = \frac{1}{2\Delta y^2}[(T_y)_{i,j+1} + (T_y)_{i,j}]$

$E = \frac{1}{2\Delta x^2}[(T_x)_{i,j} + (T_x)_{i-1,j}]$

$F = \frac{S_{i,j}}{\Delta t} \cdot h_{i,j}^{k-1}$

If the aquifer is homogeneous and isotropic, then T_x = T_y = T for all (i, j) and S_{i, j} = S for all (i, j). Under these conditions, and for a square nodal grid with Δx = Δy, the coefficients of Eq. (A9.9) become

$A = \frac{4T}{\Delta x^2} + \frac {S_{i,j}}{\Delta t}$

$B = C = D = E = \frac{T}{\Delta x^2}$

$F = \frac{S}{\Delta t} \cdot h_{i,j}^{k-1}$

If we divide through by T/Δx², these coefficients are seen to be the same as those developed in a less rigorous way in Section 8.8 and presented in connection with Eq. (8.57).

Trescott et al. (1976) have suggested that there is some advantage to utilizing the harmonic mean rather than the arithmetic mean in Eq. (A9.5). This approach changes the coefficients in the finite-difference equation but does not alter the concepts underlying the development.

Equation (A9.9) is written in terms of the hydraulic head values at five nodes at time step k and one node at time step k – 1. It is known as a backward-difference approximation. Remson et al. (1971) note that there are some computational advantages to using a central-difference approximation, known as the Crank-Nicholson scheme, that utilizes head values at five nodes at time step k and five nodes at time step k – 1. The alternating-direction implicit procedure (ADIP) utilized by Pinder and Bredehoeft (1968) involves two finite-difference equations, one in the xt plane and one in the yt plane. Each uses head values at three nodes at time step k and three nodes at time step k – 1.

Appendix X: Derivation of the Advection-Dispersion Equation for Solute Transport in Saturated Porous Media

The equation derived here is a statement of the law of conservation of mass. The derivation is based on those of Ogata (1970) and Bear (1972). It will be assumed that the porous medium is homogeneous and isotropic, that the medium is saturated, that the flow is steady-state, and that Darcy’s law applies. Under the Darcy assumption, the flow is described by the average linear velocity, which carries the dissolved substance by advection. If this were the only transport mechanism operative, nonreactive solutes being transported by the fluid would move as a plug. In reality, there is an additional mixing process, hydrodynamic dispersion (Section 2.13), which is caused by variations in the microscopic velocity within each pore channel and from one channel to another. If we wish to describe the transport process on a macroscopic scale using macroscopic parameters, yet take into account the effect of microscopic mixing, it is necessary to introduce a second mechanism of transport, in addition to advection, to account for hydrodynamic dispersion.

To establish the mathematical statement of the conservations of mass, the solute flux into and out of a small elemental volume in the porous medium will be considered (Figure A10.1). In Cartesian coordinates the specific discharge v has components (v_x, v_y, v_z) and the average linear velocity $\vec{v} = v/n$ has components $(\vec{v}_x, \vec{v}_y, \vec{v}_z)$ . The rate of advective transport is equal to $\vec{v}$ . The concentration of the solute C is defined as the mass of solute per unit volume of solution. The mass of solute per unit volume of porous media is therefore nC. For a homogeneous medium, the porosity n is a constant, and $\partial(nC)/\partial x = n\partial C/\partial x$ .

**Figure A10.1** Mass balance in a cubic element in space.

The mass of solute transported in the x direction by the two mechanisms of solute transport can be represented as

transport by advection = $\vec{v}_xnC dA$ (A10.1)

transport by dispersion = $nD_x\frac{\partial C}{\partial x} dA$ (A10.2)

where D_x is the dispersion coefficient in the x direction and dA is the elemental cross-sectional area of the cubic element. The dispersion coefficient D_x is related to the dispersivity α_x, and the diffusion coefficient D* by Eq. (9.4):

$D_x = \alpha_x\vec{v}_x + D$ (A10.3)

The form of the dispersive component embodied in Eq. (A10.2) is analogous to Fick’s first law.

If F_x represents the total mass of solute per unit cross-sectional area transported in the x direction per unit time, then

$F_x = \vec{v}_xnC - nD_x\frac{\partial C}{\partial x}$ (A10.4)

The negative sign before the dispersive term indicates that the contaminant moves toward the zone of lower concentration. Similarly, expressions in the other two directions are written

$F_y = \vec{v}_ynC - nD_y\frac{\partial C}{\partial y}$ (A10.5)

$F_z = \vec{v}_znC - nD_z\frac{\partial C}{\partial z}$ (A10.6)

The total amount of solute entering the cubic element (Figure A10.1) is

$F_x dz \hspace{1mm} dy + F_y dz \hspace{1mm} dx + F_z dx \hspace{1mm} dy$

The total amount leaving the cubic element is

$\left(F_x + \frac{\partial F_x}{\partial x}dx\right)dz \hspace{1mm} dy \left(F_y + \frac{\partial F_y}{\partial y}dy\right)dz \hspace{1mm} dx + \left(F_z + \frac{\partial F_z}{\partial z}dx \hspace{1mm} dy\right)$

where the partial terms indicate the spatial change of the solute mass in the specified direction. The difference in the amount entering and leaving the element is, therefore,

$\left(\frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z}\right) dx \hspace{1mm} dy \hspace{1mm} dz$

Because the dissolved substance is assumed to be nonreactive, the difference between the flux into the element and the flux out of the element equals the amount of dissolved substance accumulated in the element. The rate of mass change in the element is

$-n\frac{\partial C}{\partial t} dx \hspace{1mm} dy \hspace{1mm} dz$

The complete conservation of mass expression, therefore, becomes

$\frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z} = n\frac{\partial C}{\partial t}$ (A10.7)

Substitution of expressions (A10.4), (A10.5), and (A10.6) in (A10.7) and cancellation of n from both sides yields:

$\bigg[\frac{\partial}{\partial x}\left(D_x\frac{\partial C}{\partial x}\right) + \frac{\partial}{\partial y}\left(D_y\frac{\partial C}{\partial y}\right) + \frac{\partial}{\partial z}\left(D_z\frac{\partial C}{\partial z}\right)\bigg] -$ $\bigg[\frac{\partial}{\partial x}(\vec{v}_xC) + \frac{\partial}{\partial y}(\vec{v}_yC) + \frac{\partial}{\partial z}(\vec{v}_zC)\bigg] = \frac{\partial C}{\partial t}$ (A10.8)

In a homogeneous medium in which the velocity $\vec{v}$ is steady and uniform (i.e., if it does not vary through time or space), dispersion coefficients D_x, D_y, and D_z, do not vary through space, and Eq. (A10.8) becomes

$\bigg[D_x\frac{\partial^2C}{\partial x^2} + D_y\frac{\partial^2C}{\partial y^2} + D_z\frac{\partial^2C}{\partial z^2} \bigg] - \bigg[ \vec{v}_x\frac{\partial C}{\partial x} + \vec{y}_x\frac{\partial C}{\partial y} + \vec{v}_z\frac{\partial C}{\partial z} \bigg] = \frac{\partial C}{\partial t}$ (A10.9)

In one dimension,

$D_x\frac{\partial^2C}{\partial x^2} - \vec{v}_x\frac{\partial C}{\partial x} = \frac{\partial C}{\partial t}$ (A10.10)

In some applications, the one-dimensional direction is taken as a curvilinear coordinate in the direction of flow along a flowline. The transport equation then becomes

$D_l\frac{\partial^2C}{\partial l^2} - \vec{v}_l\frac{\partial C}{\partial l} = \frac{\partial C}{\partial t}$ (A10.11)

where l is the coordinate direction along the flowline, D_l the longitudinal coefficient of dispersion, and $\vec{v}_l$ , the average linear velocity along the flowline.

For a two-dimensional problem, it is possible to define two curvilinear coordinate directions, S_l and S_t, where S_l is directed along the flowline and S_t is orthogonal to it. The transport equation then becomes

$D_l\frac{\partial^2C}{\partial S_l^2} + D_t\frac{\partial^2C}{\partial S_t^2} - \vec{v}_l\frac{\partial C}{\partial S_l} = \frac{\partial C}{\partial t}$ (A10.12)

where D_l and D_t are the coefficients of dispersion in the longitudinal and transverse directions, respectively.

If v_l varies along the flowline and D_l and D_t vary through space, Eq. (A10.13) becomes

$\frac{\partial}{\partial S_l}\left(D_l\frac{\partial C}{\partial S_l}\right) + \frac{\partial}{\partial S_t}\left(D_l\frac{\partial C}{\partial S_t}\right) - \frac{\partial}{\partial S_l}(\vec{v}_lC) = \frac{\partial C}{\partial t}$ (A10.13)

Equations (A10.8) through (A10.13) represent six forms of the advection-dispersion equation for solute transport in saturated porous media. Equation (A10.11) is identical to Eq. (9.3) in Section 9.2. The solution to any one of these equations will provide the solute concentration C as a function of space and time. It will take the form C(x, y, z, t) for Eqs. (A10.8) and (A10.9); C(l, t) for Eq. (A10.11); and C(S_l, S_t, t) for Eqs. (A10.12) and (A10.13). There are many well known analytical solutions to the simpler forms of the transport equation. Equation (9.5) in Section 9.2 is an analytical solution to Eq. (A10.11), and Eq. (9.6) is an analytical solution to Eq. (A10.9). In most field situations, two- or three-dimensional analyses are required. In addition, the velocities are seldom uniform and the dispersivities are usually variable in space. For these conditions, numerical methods adapted for digital computers must be used to obtain solutions.

The coefficient of dispersion in a three-dimensional, homogeneous, isotropic porous medium as expressed in Eq. (A10.12) is a second-order symmetrical tensor with nine components. D_l and D_t are the diagonal terms of the two-dimensional form. The directional properties of the dispersion coefficients are caused by the nature of the dispersive process in the longitudinal and transverse directions. If the medium is itself anisotropic, mathematical description of the dispersive process takes on much greater complexity. No analytical or numerical solutions are available for anisotropic systems. Only for isotropic media has the coefficient of dispersion been represented successfully by experimental methods.

It is possible to extend the transport equation to include the effects of retardation of solute transport through adsorption, chemical reaction, biological transformations, or radioactive decay. In this case the mass balance carried out on the elemental volume must include a source-sink term. For retardation due to adsorption, the transport equation in a homogeneous medium in a one-dimensional system along the direction of flow takes the form

$D_l\frac{\partial^2C}{\partial l^2} - \vec{v}_l\frac{\partial c}{\partial l} + \frac{\rho b}{n}\frac{\partial S}{\partial t} = \frac{\partial C}{\partial t}$ (A10.14)

where ρb is the bulk mass density of the porous medium, n the porosity, and S the mass of chemical constituent adsorbed on a unit mass of the solid part of the porous medium. The first term of Eq. (A10.14) is the dispersion term, the second is the advection term, and the third is the reaction term. This form of the reaction term is considered in Section 9.2 in connection with Eqs. (9.9) through (9.13). Analytical solutions to Eq. (A10.14) are presented by Codell and Schreiber (in press). Figure 9.12 presents some sample solutions carried out by Pickens and Lennox (1976) using numerical methods.