Appendix X

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 \delta(nC)/\delta x = n\delta C/\delta 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{\delta C}{\delta 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 \alpha_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{\delta C}{\delta 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{\delta C}{\delta y} (A10.5)

F_z = \vec{v}_znC - nD_z\frac{\delta C}{\delta 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{\delta F_x}{\delta x}dx\right)dz \hspace{1mm} dy \left(F_y + \frac{\delta F_y}{\delta y}dy\right)dz \hspace{1mm} dx + \left(F_z + \frac{\delta F_z}{\delta 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{\delta F_x}{\delta x} + \frac{\delta F_y}{\delta y} + \frac{\delta F_z}{\delta 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{\delta C}{\delta t} dx \hspace{1mm} dy \hspace{1mm} dz

The complete conservation of mass expression, therefore, becomes

\frac{\delta F_x}{\delta x} + \frac{\delta F_y}{\delta y} + \frac{\delta F_z}{\delta z} = n\frac{\delta C}{\delta 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{\delta}{\delta x}\left(D_x\frac{\delta C}{\delta x}\right) + \frac{\delta}{\delta y}\left(D_y\frac{\delta C}{\delta y}\right) + \frac{\delta}{\delta z}\left(D_z\frac{\delta C}{\delta z}\right)\bigg] - \bigg[\frac{\delta}{\delta x}(\vec{v}_xC) + \frac{\delta}{\delta y}(\vec{v}_yC) + \frac{\delta}{\delta z}(\vec{v}_zC)\bigg] = \frac{\delta C}{\delta 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{\delta^2C}{\delta x^2} + D_y\frac{\delta^2C}{\delta y^2} +  D_z\frac{\delta^2C}{\delta z^2} \bigg] -  \bigg[ \vec{v}_x\frac{\delta C}{\delta x} + \vec{y}_x\frac{\delta C}{\delta y} + \vec{v}_z\frac{\delta C}{\delta z} \bigg] = \frac{\delta C}{\delta t} (A10.9)

In one dimension,

D_x\frac{\delta^2C}{\delta x^2} - \vec{v}_x\frac{\delta C}{\delta x} = \frac{\delta C}{\delta 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{\delta^2C}{\delta l^2} - \vec{v}_l\frac{\delta C}{\delta l} = \frac{\delta C}{\delta 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{\delta^2C}{\delta S_l^2} + D_t\frac{\delta^2C}{\delta S_t^2} - \vec{v}_l\frac{\delta C}{\delta S_l} = \frac{\delta C}{\delta 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{\delta}{\delta S_l}\left(D_l\frac{\delta C}{\delta S_l}\right) +  \frac{\delta}{\delta S_t}\left(D_l\frac{\delta C}{\delta S_t}\right) - \frac{\delta}{\delta S_l}(\vec{v}_lC) = \frac{\delta C}{\delta 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{\delta^2C}{\delta l^2} - \vec{v}_l\frac{\delta c}{\delta l} + \frac{\rho b}{n}\frac{\delta S}{\delta t} = \frac{\delta C}{\delta t} (A10.14)

where \rho_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.