:70

Theory of Dispersion 

in a Granular Medium

GEOLOGICAL SURVEY PROFESSIONAL PAPER 411-1

OD

«!

a

i

£

O

i

I?

c MAY 2 6 1970

WRD

Theory of Dispersion 

in a Granular Medium

By AKIO OGATA

GEOLOGICAL SURVEY PROFESSIONAL PAPER 411-1

A review of the theoretical aspects of 

dispersion of fluid flowing through a 

porous material

UNITED STATES GOVERNMENT PRINTING OFFICE, WASHINGTON : 1970
UNITED STATES DEPARTMENT OF THE INTERIOR 

WALTER J. HICKEL, Secretary

GEOLOGICAL SURVEY 

William T. Pecora, Director

For sale by the Superintendent of Documents, U.S. Government Printing Office 

Washington, D.C. 20402 - Price 45 cents (paper cover)
CONTENTS

Symbols... .-.-.._._...__....___.__.._.___

Abstract. --.___-_____..___-___--.__---_______.____.

Introduction- __-___-___-_____--____--_______.______

Acknowledgments..__....__._______.._._______

Historical development-_________________________

Derivation of the basic field equations_________---___

Isotropic dispersion. ..^.........................

Anisotropic dispersion.__________________________

Some useful transformations.____________________

Dispersion equation in cylindrical and spherical 

coordinates. ___.__-_-_.___..._-____.__-_-_.

Page

IV

II

1

2

2

4

5

6

7

8

Page

Field equation for two-fluid systems... .----.-..-._ 19 

Adsorption_ _ _____________________________________ 10

Mathematical treatment of the dispersion equation. _ _ _ _ 12

Nature of the dispersion coefficient.___________________ 23

Evaluation of dispersion coefficient by analytical 

methods._-________.-_____---______----._-_._ 25

Evaluation of dispersion coefficient by experimental 

methods. _-.__ .__.__________.-_.__..__ 27

Summary....____..__,-._.._.__.--__. ._..___. 32

References cited..-.._.__._,-____________ 33

ILLUSTRATIONS

Page

FIGURE 1. Diagram illustrating mass balance in a cubic element in space.-_-_-___-_--___--._-__-_- 15

2. Schematic diagram of coastal aquifer system.______________________ ..._._.__. 9

3. Schematic diagram of transport in porous medium... __.__ __ 11

4. Definition sketch of one-dimensional system..___._____-____-_-___-_-.__.____-__..-,__ 12 

5-8. Graphs showing 

5. Concentration distribution; plane source at x=0 maintained at constant concentration. __ _._---__ ___.._. _._.____...__.____________________________ 14

6. Concentration distribution; region x<0 initially at constant concentration._______ 15

7. Steady-state concentration distribution; circular surface source at constant strength. . 16

8. Steady-state concentration distribution; circular surface source maintained at constant concentration._-__-__._-__.___-.__________-_.___-____.-___-_--_-.__ 17

9. Definition sketch of radial flow._______________________.._ _._ 17

10. Graph showing concentration distribution for radial-flow system._______________________ 19

11. Definition sketch of system depicting no dispersion in x direction.__-___---___--______-- 20

12. Approximate concentration profile due to lateral dispersion; radially symmetrical model._ 20

13. Approximate concentration profile due to lateral dispersion; two-dimensional model-....- 21

14. Graph showing comparison of examples 3, 4, and 6______.-_._.____.___.._______..___ 22

15. Concentration profile showing dispersion in an adsorbing medium with 5 S/dJ= bC... _______ 24

16. Concentration profile showing dispersion in an adsorbing medium with dS/dt=b(C mS),

£_->2C/dz8 = 0, anddS ^....................................................... 25

df d*

17. Sketch showing random path chosen by a fluid particle moving through the canal system. _ 26

18. Sketch showing typical setup for column experiment._____.__..____._-__.___.-___--__- 28

19. Plot of the Schmidt number versus the Reynolds number._____________________________ 30

20. Plot of the Peclet number versus the Reynolds number.-_._....--._.__.. _ __ 31

21. Graph showing example of comparision between theoretical and experimental results.___._ 32

m
SYMBOLS

Q total rate of flow 

u, v, w velocities in the x, y, and z

directions, respectively (the 

Darcy velocity divided by 

porosity) 

C Concentration in mass of solute

per unit volume of fluid 

C0 reference concentration 

C' ratio C/C0 

Di dispersion coefficient in the i=x,

y, z directions

F{ mass flux of dissolved component per unit area in the ith 

direction 

x, y, z space dimensions of Cartesian

system 

t time

A cross-sectional area 

r radial space coordinate 

K hydraulic conductivity 

k intrinsic permeability 

p pressure 

g gravitational constant 

5 flux per unit area 

/ porosity 

d diameter of the particle making

up the solid matrix 

y0 thickness of the aquifer 

8 concentration per unit volume

of solid material

b, m constants related to linear 

_ adsorption rates 

n number of steps in a random

walk analysis 

exp exponential

Laplace transform L{f(t) } 

Inverse Laplace transform

T 

X

Per 

Pe,

1 fC =db 4mjcC+tco

erf (a;) error function defined by erf (z)

erfc(x) 

J,(x) 

Y,(x)

/»(*) 

K,(x)

J(x, y)

complimentary error function

defined by 1 erf(z) 

Bessel function of the first kind

of order , 

Bessel function of the second

kind of order 

modified Bessel function of the

first kind of order , 

modified Bessel function of the

second kind of order ,

1-e-v f%-«/o[2foO M ]ett, or 

Jo

1 -XH f " e-*-Vo[2 (tf) M ]/i

[2(zt)"]dt

length of elementary canal 

density (also r/d) 

specific weight 

viscosity

kinematic viscosity 

stream function 

transformed coordinates as

specified 

ui\x 

D/ux

udlDr, radial Peclet number 

ua/Dr, axial Peclet number
FLUID MOVEMENT IN EARTH MATERIALS

THEORY OF DISPERSION IN GRANULAR MEDIUM

By AKIO OQATA

ABSTRACT

This report is a presentation of the theoretical aspects of the 

investigations and a few confirmatory laboratory studies of the 

dispersion process. The basic assumptions associated with the 

development of the dispersion equation, based on the applicability of a heuristic expression similar to Pick's law, are described. 

In general, these assumptions limit its application to a homogeneous medium in which the dispersion process is anisotropic. 

The present level of knowledge is made apparent through the 

presentation of several solutions for the differential equation 

that approximates the dispersion process. Adsorption and its 

effect on the fate of a contaminant is given only limited discussion because of the lack of formal studies, on a macroscopic 

scale, applicable to the hydrologic regimen. The equations 

approximating dispersion in a two-fluid system (salt-fresh 

water) are also presented.

INTRODUCTION

The rapidly growing need to dispose of radioactive 

waste products, and the increased rate of pollution of the 

ground-water resource by other chemical compounds, 

make it increasingly apparent that a much more detailed description of flow in porous media must be 

sought. Because of the dangers inherent in the presence 

of contaminants in the water supply, their fate and 

mode of travel downstream from their sources must be 

predicted. The fate of a contaminant depends on both 

the macroscopic and microscopic behavior of the fluid 

under the existing flow conditions and the physicochemical conditions within the environment of the 

granular material. The discussion to follow presents a 

simplified treatment of dispersion phenomenon, and a 

mathematical expression, in the form of a differential 

equation, is developed. As the medium becomes more 

complex, the flow description becomes indeterminate 

(Skibitzke and Robinson, 1963), since the description 

of the transport mechanism, owing to convection a.nd 

dispersion, requires first a correlation between properties of the porous medium and the flow conditions.

The present-day study of dispersion may be separated 

into three categories (C. V. Theis, written commun., 

Oct. 1967): 

1. The mathematical problem.

2. The laboratory problem for confirmation of (1) and 

its extension to include parameters that cannot 

be included in the mathematical model.

3. The field that is, the natural aquifer problem 

which introduces larger variations in parameters. 

Much progress has been made in segment (1) in recent 

years, and because of availability of electronic computers, complex mathematical equations can now be 

readily analysed. The problem (segment 3) is not 

clearly understood because of large variations in the 

parameters controlling the mixing mechanisms. Segment (2), on the other hand, is effectively used as a 

method to confirm mathematical analysis, but as yet, 

is limited in its use as a physical model of the real 

system. Bachmat (1967) pointed out that both the 

geometric and time-scale ratio for the model and 

prototype must be 1:1.

This paper is limited to the theory of dispersion or 

category (1) as described above. Because existing 

literature on the subject is extensive, a thorough review 

is beyond the scope of this paper.

The mechanisms associated with the transport of a 

dissolved substance in porous media are listed below. 

Some of these mechanisms are of little significance in 

natural ground-water bodies, but they are important 

in exchange systems that are used in the chemical 

industry. The mechanisms are as follows:

1. Molecular diffusion the transport of mass in its 

ionic or molecular state due to differences in 

concentration of a given species in space. The 

gross transport obeys Fick's first and second laws 

of diffusion (Bird and others, 1960, p. 502).

2. Mechanical dispersion the mixing mechanism that 

is present because of the variations in the microscopic velocity within each flow channel and from 

one channel to another. Microscopically, there is 

no mixing; however, if the average concentration 

of a given volume of fluid is considered, an apparent dilution or spreading is present. The gross 

transport, as in eddy diffusion, is given by an 

expression similar to Fick's law.

li
12 FLUID MOVEMENT IN EARTH MATERIALS

3. Eddy diffusion the mixing process that is due to 

the random fluctuation of fluid mass or the occurrence of eddies in the condition described as 

turbulent flow. In porous media these eddies exist 

only to a slight extent, if at all, because of the 

small size of the pores; but in media consisting of 

large grains, such as gravel, eddies may occur to 

some degree. In the more common types of media, 

however, the existence of low gradients tends to 

limit or preclude turbulence.

4. Mixing due to structural controls the apparent 

mixing mechanism due to large-scale structural 

variations in the granular material. These variations control the direction of movement of a given 

fluid particle; hence, if an average concentration 

is taken along a given plane parallel or transverse 

to the direction of flow, a large-scale mixing is 

observed. Whether this type of mixing can be 

treated in a manner analogous to that of mechanical dispersion is questionable, and, in any event, 

such a treatment would likely require a complete 

description of changes in the properties of the 

porous medium. The analysis would probably 

require the use of the Lagrangian concept in 

which each fluid particle must be traced along its 

path of flow.

5. Adsorption a process which differs from the others 

in that the amount of mass transported depends 

on the physicochemical interaction of the transported substance and the solid of the medium. 

The existence of an unbalanced force field causes 

migration of a liquid-borne contaminant from the 

liquid to the solid surface. In some instances the 

contaminant is permanently fixed on the surface, 

whereas in others it moves continuously from 

liquid to solid, or vice versa. Hence, the measurement of the liquid phase concentration shows a 

dimunition of the concentration. Also, the contaminant does not travel at the same rate as the 

fluid.

At present, there are basically two related methods of 

analysis that are useful in depicting dispersion. The 

first is a microscopic study of a fluid particle which 

employs the concepts of statistical analysis. The 

second is a microscopic analysis of dispersion or transport of mass which is expressed in terms of a heuristic 

law analogous to Fick's first law of diffusion. The 

microscopic analysis, which predicts the position of the 

fluid particle at any given time, is completely dependent 

on the mathematical model chosen to represent the 

medium. The most satisfactory model to date is the 

representation based on the assumption of pore channels located with complete randomness and multiply 

connected at some end points. On the other hand, the

heuristic approach requires that a partial differential 

equation similar to all transport systems be solved.

Acknowledgments

This review of the Theory of Dispersion was initiated 

while the author was associated with the Phoenix 

Research Office of the U.S. Geological Survey. The 

author gratefully acknowledges the detailed discussions 

and suggestions of R. H. Brown, H. H. Cooper, H. E. 

Skibitzke, C. V. Theis, and others in the review of the

manuscript.

Historical Development

Perhaps the first to recognize the occurrence of dispersion phenomena was Slichter (1905) in his investigation of the use of salt as a tracer in ground-water flow. 

He noted that at a point away from the source the salt 

concentration increased gradually rather than abruptly, 

as would be predicted by Darcy's expression. That is, a 

gradation of salt concentration in the fluid occurred as 

the front moved past a given point of measurement. 

Slichter explained this by noting that in flow through 

capillary tubes the velocity of the fluid varies across 

the cross section of each tube and that, because the soil 

complex is composed of a great number of these tubes, 

the sum of the deviation in each tube is likely to cause 

the mixing that he observed. Following the same line 

of reasoning, lateral dispersion dispersion perpendicular to the direction of the average flow was 

attributed to the repeated branching of the individual 

flow tubes within the solid matrix. This explanation is 

most commonly used today to describe the mechanical 

component of the dispersion phenomenon.

Earlier papers dealing with two-fluid systems have 

been concerned with the economical^ important 

problem of salt-water invasion in coastal aquifers. 

The first analytical descriptions and investigations 

were based on the supposition that the fluids could be 

treated as though they were immiscible. However, 

in the middle and later 1930's, a zone of diffusion 

which could not be explained by immiscible-fluid 

theory, was noted. This diffusion zone prompted 

recognition of the dispersion process as a possible cause 

for its occurrence. This phenomenon was treated 

theoretically by Kitagawa in 1934 when he made 

theoretical and experimental studies of dispersion in a 

porous system. The physical model he chose was that 

of dispersion of a salt tracer emitted from a line source 

within a horizontal flow field. The experimental model, 

like most laboratory models, was composed of sand. 

The theoretical analysis was based on the assumption 

that, at some distance from the source, the concentration distribution would conform with the normal 

probability curve.
THEORY OF DISPERSION IN A GRANULAR MEDIUM 13

Before Kitagawa's treatment became widely accepted, refinements of the displacement theory for a 

coastal aquifer came about largely through the efforts 

of Wentworth (1948). He postulated a "rinsing theory" 

based on the hypothesis that the pores of the medium 

act as mixing cells and that the ocean tides cause the 

alternate filling of the pore by salt and fresh water, 

thereby creating a partial mixing of the fresh and salt 

water. Wentworth derived an approximate expression 

for numerical computation.

Accelerated field investigations in the latter part of 

the 1950's and in the early 1960's indicated that instability of the salt water-fresh water front posed a 

threat to the water-supply systems in many coastal 

regions. The results of these studies were recorded in 

various publications; for example, Cooper (1959) and 

Kohout (1960, 1961) described the conditions in 

Florida with supplementary information on the Hawaiian Islands to provide evidence of a dispersion-induced 

circulation of sea water into the zone of diffusion and 

back to the sea.

The results of these investigations gave impetus to 

the development of analytical methods designed to 

predict the dispersive effects of ocean tides. Results of 

other analytical studies have been published by Carrier (1958), Jong (1959), Henry (1960), and Bear and 

Todd (1960). The limitations of these investigations 

lie in the incomplete knowledge of the flow system in 

heterogeneous coastal aquifers and in the necessary 

oversimplification of such a system. Hence, the theoretical results were correlated only with those from a 

laboratory model designed to conform with the mathematical model.

It was argued that, because of the discrepancy between the medium of the mathematical model and the 

real medium, the mathematical model cannot be expected to depict what actually took place. Because of 

this several investigators have preferred a statistical 

model to develop expressions for the dispersion coefficient. Two statistical models were used: the completely disordered flow of Scheidegger (1954) and the 

randomly oriented flow channels of Jong (1958) and 

Saffman (1959, 1960). Experiments on laboratory 

models, however, indicated a variance in the completely 

disordered model which could not be reconciled with 

the laboratory data. An important limitation of 

Scheidegger's model was that it did not include a 

directionally variable coefficient of dispersion. The 

assumption of random fluid paths did, however, provide 

for the directionally variable coefficient, because of the 

lateral components of the flow velocity. The difficulties 

that precluded a model description of the aquifer were 

apparent, and a variety of postulates were introduced

to simplify the flow regime so that an analysis was 

possible. The assumption that the medium can be 

described limited the analysis to an isotropic medium.

Most of the articles published assumed that the 

field equation, similar to Fick's second law of diffusion 

in a moving medium, was a good first-order approximation. Although the assumption that the transport 

equation consists of two modes the Darcy component 

and the dispersion component avoided the requirement that the individual flow paths bo known, it required that the average direction and magnitude of the 

flow be defined and that the dispersion coefficient be 

determined experimentally. Because the flow through a 

columnar model is readily defined, this type of analysis 

is especially useful in predicting the behavior of reactor 

columns or ion-exchange columns employed in the 

chemical industry. Dankwerts (1953) gave a comprehensive discussion of the analysis of a columnar model 

and described the types of flow found in experiments 

on packed columns. In most industrial applications, the 

flow rate in the exchange process is much higher than 

that of ground water in nature; therefore, the analysis 

may not be useful in the field of hydrology. The differential equation which describes the field is valid, but 

the nature and magnitude of the dispersion coefficient 

are altered.

An analysis which employs a "perfect mixture 

model" described, for example, by Aris and Amundson 

(1957), may also be used for the higher fluid velocities. 

This mathematical model assumed that each pore cell 

within the porous medium constituted a mixing 

chamber, within which the two fluids were perfectly 

mixed. The analysis differed from that of the low-flow 

model in that the laboratory experimental data indicate a concentration breakthrough curve which was 

more nearly approximated by the Poisson distribution.

In the middle and later 1950's there were many 

studies of dispersion in sand-filled columns, which were 

related to hydrology. Day (1956) investigated a saltwater replacement system and based his analysis on 

the diffusion equation and Scheidegger's disordered 

model to explain data obtained in the laboratory. 

Rifai and others (1956) published what was perhaps 

one of the most comprehensive laboratory and theoretical studies to that time. The theoretical part of 

their study employed both the capillary-tube and the 

disordered-medium models. As was expected, the two 

models gave identical results because of the applicability of the central limit theorem (Beran, 1957) in the 

disordered-medium model. In essence, the application 

of the central limit theorem indicated that the diffusion 

equation was valid. Both laboratory and theoretical 

results indicated that the dispersion coefficient in the
14 FLUID MOVEMENT IN EARTH MATERIALS

direction of flow was nearly proportional to the first 

power of the average velocity, which confirmed 

Scheidegger's prediction.

The results of studies by Rifai and others have been 

largely confirmed by many investigators; however, 

there are some indications that the dispersion coefficient is proportional to a power of the velocity which 

is slightly greater than one. It will be shown in a subsequent section that, if the flow similitude is specified by 

the Reynolds number and the mass transport by the 

Schmidt number, this power is very near one. Studies 

of friction losses in consolidated media (Muskat, 1946), 

indicated that the Reynolds number, unless modified 

to some extent, did not adequately specify the similitude conditions. Efforts to extend the Reynolds number 

concept to flow in other than isotropic unconsolidated 

material have not been successful.

Experimental evidence indicated that the magnitude 

of the dispersion coefficient was not equal in all directions. In fact, the very nature of the flow system makes 

it highly unlikely that the magnitudes would be the 

same. Hence, Jong (1958) and Saffman (1959, 1960) 

developed a model of random flow paths to compute 

the magnitude of the dispersion in directions both 

parallel and transverse to the flow. The mere existence 

of individual grains in porous material gave rise to 

microscopic velocity components transverse to the 

direction of the average flow, thereby creating the 

lateral movement of the contaminant and the velocity 

dependency of the dispersion coefficient.

This concept was extended by Bear (1961) in a paper 

on the tensor form of the dispersion coefficient. He 

verified the tensor characteristic of the coefficient 

experimentally by observing the dispersion of a tracer 

from a point source. Although there was little doubt 

that the coefficient was a tensor, Bear's experiment did 

not correct for the spread due to the volume of the 

tracer injected. On the basis of experience in other 

fields, such as heat flow and molecular diffusion, it can 

be assumed that the generalized transport coefficient 

may be characterized as a tensor.

The magnitude of the dispersion coefficient was 

dependent not only on the fluid velocity but also on the 

characteristics of the porous medium. Many experiments reported in the literature show this to be true. 

Among these were experiments by Orlob and Radhakrishna (1958), who explored the effect of entrapped 

air in a granular material, and by Biggar and Nielsen 

(1960), who examined sands and other types of naturally occurring soils to determine the effect of inhomogeneity and variations of particle sizes. Also, 

Skibitzke and Robinson (1963) obtained qualitative 

laboratory results that demonstrated the importance of 

heterogeneity of the medium.

To date, both theoretical and experimental investigations have been confined primarily to unidirectional flow. 

Although the experimental data indicated that dispersion was directionally dependent, they were not sufficient to enable a statement as to how the ratio of the 

magnitudes of the two components of the coefficient can 

be predicted in advance of an experiment. To depict 

transverse dispersion, Ogata (1961) developed a simplified mathematical model for steady flow, principally 

to analyze experimental results given in a later report 

(Ogata, 1964a; also Skibitzke, 1964). Harleman and 

Rumer (1962) used a similar mathematical model to 

correlate results of their experiments. Many results 

published in chemical journals will not be mentioned 

because they involved flow rates higher than those of 

interest here. Articles have appeared also in petroleum 

journals, but most of these articles are reviews of 

previous works generally with additional experimental 

data to supplement results of others so that empirical 

expressions can be developed.

DERIVATION OF THE BASIC FIELD EQUATIONS

The equation derived in this section is a statement of 

the law of conservation of mass. For simplicity, the 

development will be in terms of Cartesian (x, y, z) 

coordinates. The field equations can then be obtained 

readily in terms of other coordinate systems through 

known transformation.

The description of the mechanism of mass transport 

through a complex porous medium necessitates the 

adoption of a fictitious model such that average conditions can be more easily expressed. The first general 

assumption is that Darcy's law applies. This is corollary 

to the assumption that a piston-type flow occurs. In 

other words, the model is such that all particles of 

water move at the same rate of speed through the 

porous medium. Because of this assumption, which 

neglects microscopic variations in fluid velocity, a 

refinement must be employed to approximate the real 

system. For the dispersion mechanism, this refinement 

involves a transport component that is analogous to 

Fick's first law. The justification for the use of this 

type of expression is given by Taylor (1953) in his 

study of laminar flow through pipes.

The dispersion coefficient that appears hi the dispersion component is assumed to be approximated by a 

constant, whereas it generally is to some extent dependent on the concentration and the type of dissolved or suspended substance. The development of 

the differential equation, however, is based on the 

assumption that the concentration is sufficiently small 

that the dispersion coefficient is independent of concentration. No chemical reaction is assumed to occur 

between the solid and the liquid phases. This supposes
THEORY OF DISPERSION IN A GRANULAR MEDIUM

that within the fluid system no loss or addition of 

matter can take place.

ISOTROPIC DISPERSION

To establish the mathematical statement of the 

conservation of mass, consider a cubic element hi the 

three-dimensional Cartesian space, as represented in 

figure 1. The two modes of transport of a fluid within

bx

FIGURE 1. Mass balance in a cubic element in space.

a porous medium have been described as (a) convection 

and (b) dispersion. In mathematical terms these modes 

for direction x (say) :are expressed:

Transport by convection='wC'eLl

dC Transport by dispersion=D,,. ^- dA>

where dA is an elemental cross-section area of the 

cubic element. Assuming that these two components 

may be superposed, the total amount of material 

transported parallel to any given direction is obtained 

by summing the convective and dispersive transports. 

Thus, if Fx represents the total amount of mass per 

unit cross-sectional area transported in the x direction 

per unit tune, then

where / is the porosity of the medium. The negative 

sign before the dispersive term indicates that the 

contaminant moves toward the zone of lower fluid 

concentration.

370-837 O 70 2

Similarly, expressions in the other two directions are

written

Fv=vC-JD dC

From figure 1, the total amount of solute entering a 

cubic element hi space is

F* dzdy+Fv dzdx+F, dxdy. 

The total amount leaving the cubic element is

dy\ dzdx

where the partial terms indicate the spatial change of 

the fluid mass hi the specified direction. Accordingly, 

the difference in the amount entering and leaving the 

cubic element is

Because there is no loss in the mass of the liquid, the 

difference between the amount entering and leaving 

must be equal to the amount of mass accumulated 

within the element. The rate of mass change is represented mathematically by the expression

The equating of the difference of outflow and inflow to 

the amount of dissolved substance accumulated within 

the cubic element leads to the relationship

dF, dF, .^ ~ (1)

Equation 1 is a mathematical statement of the law 

of conservation of mass under the conditions stipulated. 

Vectorially, the equation is written

V-F=JdC/dt.

Substituting the expressions for Fx, Fv, and Ft and 

assuming that Dx=Dv=De=D=constant and since 

V-q=0, equation 1 becomes

dt (2)

j a a a a

where-37is the operator-rr-fw-5 \-v-~ 1-10-5 Equation 2 dt r dt ox oy oz
16 FLUID MOVEMENT IN EAETH MATERIALS

applies to both diffusion and heat transfer within a 

moving medium. Note that the terms u, v, and w as 

used in equation 2 now represent average velocities or 

the Darcy velocity divided by the porosity.

Because the characteristics of the dispersion coefficient other than its directional nature are not well 

understood, equation 2 was developed on the assumption that D was a constant. If D is a function of concentration, that is, if D=/(<7), equation 2 is written

Equation 3 is nonlinear for D=/(<7) and can be solved 

only for simple expressions of j(O), generally by 

numerical methods. These types of problems are encountered in the analysis of soil moisture and in some 

other types of diffusive processes. Excellent discussions 

of the numerical methods applicable in obtaining solutions for the one-dimensional equations are presented 

by Crank (1956) and also in many articles relating to 

soil-moisture movement.

ANISOTROPIC DISPERSION

All experiments conducted in two dimensions indicate 

that the magnitudes of the dispersion coefficient depend on the reference direction, with the larger value 

oriented in the direction parallel to the flow. To include 

this directional dependency in the equations requires 

that the coefficient D be represented by a more general 

mathematical construct known as a tensor. The complex properties of tensors are better described in texts 

on advanced mathematics and need not be included 

here.

In generalizing the dispersion coefficient into a tensor, 

Bear (1961) assumed an analogy with a tensor representing elastic deformation or tensor of the fourth 

rank. In three dimensions this tensor has 81 components, but because of symmetry the tensor has 21 

nonzero components (Jeffreys, 1961, p. 78). Scheidegger 

(1961), in a verification of Bear's tensor, assumed dispersion to be analogous to turbulence rather than to 

elasticity. Because turbulence is not symmetrical to 

the same extent as elasticity, its tensor consists of 

36 rather than 21 nonzero components. For unidirectional flow in an isotropic medium, however, the 

results of both Bear's and Scheidegger's analyses show 

that the dispersion coefficient is described completely 

by longitudinal and transverse components in other 

words, by a tensor composed of only two components.

The difficulties inherent in the application of the 

tensor to evaluate mass transport arise from difficulties 

in measuring the various components. Thus, as hi heat 

flow or diffusion, it is generally necessary to -assume

that the dispersion coefficient is characterized by three 

independent components parallel to the chosen reference axes. Under this assumption the dispersion tensor 

is a second-rank tensor consisting of nine components. 

Jost (1960, p. 4), in discussing diffusion systems, 

argued that the existence of components other than 

those making up the main diagonal of the matrix is 

doubtful. However, this argument cannot readily be 

extended into the study of dispersion in granular 

media. Thus, for simplicity it is assumed that the 

dispersion coefficient is approximated by a tensor 

composed of, at most, nine components.

Using the standard notation for second-order tensors, 

the dispersion component of the transport equation 

can be expressed

 .Gi =DiJ i j=12tB. (4)

In other words, the three components of mass transport 

are written

dC

and the dispersion tensor Dtf can be represented by the 

matrix

The advantage of the tensor notation is that it 

provides a shorthand method of describing, in general, 

the physical phenomena. Equation 4 describes the 

most useful case of dispersion in a porous medium and 

indicates that the three components of mass transport 

are determined not only by the direction of the maximum concentration gradient but also by the gradients 

at right angles to this direction. In an anisotropic 

medium, because the dispersion constant is dependent 

on the fluid and medium properties, the elements 

forming the matrix can be both time and space dependent. The dispersion tensor can be described by a 

single number, as in isotropic dispersion, or by a 

matrix which reflects the complexity of the porous 

medium itself. Only for isotropic media has the coefficient of dispersion been described successfully by 

experimental methods.

If the medium is isotropic and the Darcy velocity is 

constant, the transport equation in the x direction is 

Fx=uC Qx ; where Qx is the dispersion component.
THEORY OF DISPERSION IN A GRANULAR MEDIUM 17

The substitution of the transport equation into the 

continuity equation results in the field equation

t d*C 

dxdy (5)

Equation 5 describes the field distribution for a system 

of anisotropic mass transport. If Dn=DM =Da& =constant and Z^=0 for i^j, equation 5 is identical with 

equation 2. This results from the fact that equation 2 

was developed under the assumption that D was a 

tensor of zero rank.

Equation 5 is known as a quadric, and by use of 

standard transformations it can be reduced to the form 

of equation 3. This transformation involves rotating 

the coordinate axes so that the reference axes parallel 

the principal axes of dispersion. Recent experimental 

and analytical studies point to the important fact that 

hi isotropic and homogeneous media the principal axes 

of dispersion are oriented parallel and transverse to 

the mean direction of regional flow. This indicates that 

for such media the mass transport system can be 

defined by two characteristic dispersion components 

that are specified when the mean direction of regional 

flow is known.

Assuming that the principal axes can be defined, the 

dispersion tensor can be transformed so that only the 

elements of the major diagonal remain, all others being 

zero. The matrix representation of the tensor then 

becomes

'DX Q 0

0 Dy 0

0 0 D.

In unidirectional flow, symmetry about the mean flow 

line exists so that DV=DZ .

The possible representation of the dispersion coefficient in the manner depicted by the matrix is significant 

because of the internal complexity of the porous environment. Actually, the dispersion of a fluid flowing 

through a homogeneous porous medium is the result of 

the averaging process necessary to describe the flow 

of ground water. In essence, this signifies that provided 

all the possible microscopic changes in the magnitude 

and the direction of the pore channels can be expressed 

quantitatively, the components of the dispersion tensor 

can be computed by statistical means.

SOME USEFUL TRANSFORMATIONS

The diffusion, or heat-flow, equation where grad

g<7=0 has been studied extensively and published in 

detail by Carslaw and Jaeger (1959), Crank (1956), 

and Jost (1960). For various types of bourdaryvalued problems it is expedient to transform the field 

equation so that the convective terms do not appear 

explicitly.

Two methods can be used to remove the explicit 

convective terms for problems in which the boundary 

conditions involve media of infinite extent. These 

transformations require the introduction of a new 

coordinate system or a change in the independent 

variables. First consider the transformation of a 

standard-rectilinear coordinate system into a new 

system called moving coordinates. The dispersion equation including the directional properties of coefficient 

D\a

dC -£7-\-U . dC 3 \-V . dC. 3 \-W dC-~-=

dt dx dy dz

 ^- dz2

5-0

oy

(6)

A transformation similar to one commonly used in the 

study of wave equations is given by the expressions

£=x ut, r}=y vt, £=z wt, and r=t.

Using the chain rule of partial differentiation, the 

above expressions are substituted in equation 6. For 

example, the transformation of dC/dx is given by

dx dx

Two similar expressions for rj and £ coordinates can be 

obtained in the same manner. The tune differentiation 

follows the same rule and is given by

dC_dCd![ ,dCdr_dC_ dC 

dt~d^drdrdt dr U dl-'

A similar method is utilized to determine the second 

differentiation with respect to the x, y, and z coordinates. Substitution in equation 6 reduces it to the form

dC 

6T Q 2

(7)

which is analogous to the heat conduction or diffusion 

equation. The transformation of coordinates did not 

involve rotation of the coordinate axes; hence, the 

dispersion coefficient D< is still oriented in the original 

x, y, and z or the new £, 77, and f directions.

Introduction of the above transformation in effect 

creates a system based on coordinates that move in
18 FLUID MOVEMENT IN EARTH MATERIALS

space at a rate dependent on the magnitudes of the 

average fluid-flow velocities. Physically, this may be 

pictured as an observer moving along with a velocity ~g 

while watching the process of dispersion. The choice of 

this new coordinate system is possible because the 

transport due to convection is a conservative process; 

hence, stopping the fluid movement to simplify the 

analysis of the dispersion process is permissible. 

However, the transformation does pose difficulties as, 

for example, hi problems where certain conditions are 

specified at an initial point. In such problems the 

transformation tends to increase the difficulty of 

inserting the specified conditions in the solution of the 

differential equation. To illustrate, consider the problem 

where the concentration C is specified at the point 

x=y=z=Q; say, C(Q,Q,Q,f) = C0, where <70 is a constant. 

It is noted that when x, y, z=Q; £, TJ, f= qt, where 

q=u, v, w, respectively. The problem now becomes one 

of describing the concentration for a system in which 

the concentration is specified at points hi space depending on the time. This tends to increase the complexity 

of the mathematical analysis, unless there is symmetrical distribution of concentration about the source 

of fluid at 3=0.

The second transformation method involves the 

substitution of a new dependent variable. In this 

instance, it is assumed that the solution of the problem 

can be expressed functionally by

u

_w_ "2A

where T(x, y, z, t) is an unspecified function. Assuming 

that, generally, the solution of the diffusion equation is 

of this nature, the next step is to determine the function T(x, y, z, t). The solution must satisfy equation 6; 

hence, substitution gives the expression

~*\7=:L'f

32r

(8)

which is the diffusion equation. The boundary conditions, as given in the preceding paragraph, must be 

written as a time-dependent variable at a specified 

point in space. Because F0 (say) is specified at a given 

point in space, an analysis using equation 8 introduces 

no further difficulties.

In general, the transformations leading to equations 7 

and 8 work well for dispersion in media of infinite extent. In the event of a finite boundary the concentration is specified along x, y} z=L some other transformation of equation 6 must be utilized to circumvent 

the specification of concentration at infinity.

Usually, mathematical solutions are readily found 

for the dispersion equation describing isotropic conditions; that is, for Dx=Dy=Dg=D. Thus, introducing 

the independent variable

*'=*A/ . 2/'=2/VrT> 2'=ZV7T

and substituting these independent variables in equation 6 gives

where

i, j, fc=unit vectors hi the x, y, and z directions. To 

reduce this relation to the form of equations 7 and 8, 

the transformations previously discussed can be used.

DISPERSION EQUATION IN CYLINDRICAL AND 

SPHERICAL COORDINATES

The field equation as previously written is in rectangular or Cartesian coordinates. For some problems the 

cylindrical or the spherical coordinates may be more 

suitable for analysis. The fundamental equations 

governing the dispersion process are readily derived in 

much the same manner as before by stating the mass 

balance within a reference space as defined by the given 

coordinate system. Similarly, coordinate transformations are readily accomplished because the interrelationships between any given coordinate systems 

are known.

In cylindrical coordinates the isotropic dispersion 

equation is

dc BO ac do

And in spherical coordinates the equation is 

dC . dC. /l\atf, / 1 \dC

de

(10)
THEORY OF DISPERSION IN A GRANULAR MEDIUM 19

FIELD EQUATION FOR TWO-FLUID SYSTEMS

In the development of the dispersion equation the 

density and viscosity of the contaminant fluid was 

assumed to be the same as the surrounding original or 

uncontaminated fluid. In cases where high concentrations of a contaminant occur, as in an aquifer subject 

to salt-water encroachment in coastal regions, this 

assumption is no longer valid. In such an environment 

both the viscosity and the specific weight of the contaminated fluid change. Addition of the gravity force 

is now necessary for analyzing the concentration distribution along the interface between the contaminated 

and uncontaminated fluids. The gravity force exerts so 

large an influence on the flow field that the position of a 

particle of contaminated fluid at any given time is 

largely predetermined by this force. That is, the 

description of the velocity field is dependent on the 

dynamic equilibrium position of the interface between 

the two fluids. Thus, definition of the dispersion process 

in a two-fluid system requires complete knowledge 

of both the fluid properties and the flow field.

When no chemical reaction between the fluids and 

the solids take place, the dispersion equation is the 

same as previously described with the provision that 

velocity is space dependent. Hence, the field equation 

is written

dC (11)

Equation 11 indicates that the field equation is completely specified provided the flux g is determined. 

Although conditions depart somewhat from those 

studied by Darcy, it is assumed that the flow equation 

may be generalized such that

2=- (V#- K f*, (12)

The schematic diagram of a two-fluid model under 

dynamic equilibrium conditions is shown in figure 2. 

The description of the transport of salt requires the 

complete solution of equation 11 subjected to equation 

12. Inasmuch as equation 11 is nonlinear it is virtually 

unsolvable without the use of numerical methods or 

various simplifications. For example, Jong (1959) 

defined the flow system in terms of curvilinear coordinates oriented with respect to the assumed position of 

the interface, thereby enabling the use of a twodimensional system. Paralleling the classical hydrodynamics, a stream function is introduced into the 

analysis. The pertinent components of the stream 

function are defined as

Land surface

Sea level

Fresh water

Zone of diffusion 

Interface

FIGURE 2. Schematic diagram of coastal aquifer system.

where s and n denote the directions parallel and normal 

to the interface, respectively, as shown in figure 2. The 

components of flow in the specified reference frame 

(eq 12) are

_- 

ds y ds

(13)

where k= - Accordingly, substituting the stream

pg 

function in equation 13, the continuity equation may

be written

The terms in the brackets are included to show the 

effect of viscosity M which is dependent on the salt 

concentration in the fluid. The right side of the equation 

describes the effects of changes in density and viscosity. 

Equation 14 is analogous in form to a vorticity equation that frequently appears in the potential theory in 

hydrodynamics. This vorticity equation has been 

studied in classical hydrodynamics, and methods of 

solution for various boundary conditions are described 

in a number of texts; for example, Milne-Thomson 

(1961, p. 564). Equation 14 is also substantially the 

same equation derived by Henry (1960) in evaluating 

the two-fluid system. Henry used the same basic flow 

equation, expressed the specific weight 7 in terms of the 

salt concentration (7, and assumed that the viscosity
no FLUID MOVEMENT IN EARTH MATERIALS

change due to that concentration is small. Hence, the 

equation he derived is in the form

(15) Q dz'

where ki=k ( ), y0 is the aquifer thickness, Q is the 

\ Po /

aquifer discharge, C' is the ratio between the concentration of salts in the aquifer water and the sea water 

C/CB, and xf is the dimensionless length x/y0 . Equations 

14 and 15 differ only in that the first specifies possible 

effects of the viscosity changes.

To simplify the analysis, Jong (1959) dealt with the 

vorticity equation (eq 14), excluding the terms in 

brackets since the change in viscosity is small. He 

showed that the analysis can be carried out by replacing 

real fluids with fictitious fluids and by introducing 

singularities parallel to the interface of the system. 

The characteristics of the singularities are determined 

by the actual observed velocity distribution; in other 

words, the singularities are chosen so that the original 

velocity distribution in the two fluids is reproduced. The 

method, thus, makes possible computation of the position of the interface between two fluids of different 

densities. Because the equations are derived for steadystate conditions, the final position of the interface can 

be computed for any given flow conditions.

Equation 15 describes the flow conditions, in terms 

of a stream function, for the change in specific weight 

of the fluid with the change in salt content. However, 

the transport of the salt is also dependent on mechanical 

dispersion; thus, equation 15 specifies only a part of 

the problem. In other words, equation 1.5 describes the 

flow field within the transition zone containing fluids of 

different densities, and the transport of mass within 

this zone needs to be described by the dispersion 

equation. The solution of the vorticity equation given 

by Jong makes it possible to compute the location of 

the interface for any specified time, from which the 

stream function can be determined for the whole flow 

field. When unsteady conditions prevail, these solutions 

may still be used as long as the changes in the flow 

field are small. The transport system, however, must 

be described subsequent to description of flow in order 

to determine the dispersion effects on the frontal zone.

vSince by definition q*=6V/dy and qv= 6V/dx, substitution into the two-dimensional steady-state transport or dispersion equation leads to

Q

dy' 6V (16)

Equation 16 is thus an approximation of the steadystate mass transport process which can be completely

specified, provided * is specified by equation 15. Henry 

(1960) has developed solutions for two special instances 

of the transport equation and has evaluated them 

numerically.

ADSORPTION

In the discussion and development of dispersion 

equations to this point, no means of mass transport has 

been considered other than mechanical dispersion and 

convective transport. In the isotropic systems usually 

studied in the laboratory, the spread of tracers in the 

fluid is generally observed to be small, which indicates 

that the magnitude of the dispersive term is small. All 

soil complexes display an additional component of mass 

transfer generally termed "adsorption." Adsorption is 

a chemical reaction in which mass transfer occurs by 

actual removal of the dissolved substance from the 

liquid phase owing to mutual attraction between the 

substance and the solid phase of the porous medium. 

Generally speaking, in natural flow systems this process 

is of major magnitude with the fate of any contaminant 

introduced into the soil or the ground-water system 

dependent to a large extent on the capacity of the solid 

matrix material to adsorb the dissolved component. 

This is recognized by most workers in the field of radioisotope disposal; however, the major research is confined to surface chemistry which takes into account the 

microscopic aspects of the adsorption process.

Vermeulen and Hiester (1952), in a macroscopic investigation of ion-exchange columns, specified that the 

ratio of the tracer velocity to the main body of fluid 

velocity under equilibrium conditions may be expressed

by

UT 1 iLf Sp.

U

where

S= trace component adsorbed per gram of solid

matrix,

p=bulk density of solid matrix, and 

UT=velocity of tracer.

This expression was established empirically and has 

been used to some extent to determine rate of movement of radioactive isotopes.

To develop the transport equation which includes 

adsorption, it is assumed that the rate of adsorption 

is the third component in the dispersion equation. This 

assumption simply implies that the chemical reaction 

or the rate of adsorption can be superimposed on the 

dispersion and convection system previously described 

by equation 3. The adsorption component in the physical model can be pictured as a mathematical sink or 

source with its strength dependent on the difference of 

concentration of the liquid phase and the solid phase. 

A simplified model is illustrated in figure 3. Here, it is
THEORY OF DISPERSION IN A GRANULAR MEDIUM 111

Total inflow

Solid 

phase

Liquid 

phase

dx

Convection and 

dispersion components

(From Bird and others, 1960)

FIGURE 3. Schematic diagram of transport In porous medium.

assumed that the adsorption component is in one direction only.

Superposition of the three components adsorption, 

dispersion, and convection allows the writing of the 

continuity or mass-balance equation as follows:

(17)

where S is defined simply as the concentration of the 

solid phase. The function 8 has the same dimension 

as C.

The analysis of dispersion in adsorbing media that 

is, the use of equation 17 depends on the representation of S in terms of C, where it is assumed that the 

relationship £)$/()<=/(<?) holds. To describe this phenomenon completely requires a macroscopic study of 

the energy state at any equilibrium condition as prescribed by thermodynamic principles. Because of the 

complexity of the subject, especially when the kinetics 

of transfer rates is included, no detailed discussion is 

presented. Excellent discussions are included in Glasstone (1961) and Adamson (1960). This discussion is 

limited to the use of simplified expressions that have 

been applied in studies of reaction columns.

One of the expressions most frequently used in the 

analysis of ion exchange is

(18)

where b and ra are constants. As indicated by Crank 

(1956), this equation describes a first-order reversible 

reaction and is applicable to some ion-exchange processes in fluids flowing in porous substances. Note that, 

when the term mS is much smaller than C, the expression reduces to a form similar to Langmuir's equilibrium 

isotherm for low concentrations or

(19)

Equation 19 is applicable to a first-order irreversible 

reaction wherein the soil particle may be thought to 

act as a mathematical sink.

Extending the expression given as equation 18 a 

little further, it may be written

where n< 1. Note that for the special case of n= 1, this 

is equation 18. The need to determine three unknown 

parameters makes the above equation difficult to 

evaluate even under the rather idealized conditions 

encountered in a laboratory. Again assuming that 

mS^C, the equation may be written in the simpler 

form

%j-=bCn > (20) d\

which is virtually the expression for Freundlich's 

equilibrium isotherm.

The three types of reaction rates represented by 

equations 19-21 are discussed by Crank (1956, p. 121) 

in reference to a diffusion system involving various 

boundary conditions. The solution of the simultaneous 

equations represented by equations 3 (for a constant 

D), 18, and 19 can be determined by operational 

methods. A numerical method, however, must be 

used to solve equation 20. When equation 19 is used 

to represent the reaction rate, direct use can be made of 

solutions for various boundary conditions in heat flow 

reported extensively by, for example, Carslaw and 

Jaeger (1959). The reversible reaction system, equation 

20 for diffusion in a plane sheet, has been analytically 

determined by Crank (1956, p. 132) for the boundary 

conditions

£=(7=0, -a<x<a,t=Q

x, x=±a,

The solution, however, is too complex to be reproduced. 

Ogata (1964b) obtained a solution for a semi-infinite 

medium, which will be presented as one of the mathematical examples.

Despite the fact that the adsorption process is 

generally ill-defined, it is an important factor in the 

control of mass transport within a porous medium. In
112 FLUID MOVEMENT IN EARTH MATERIALS

most cases the fate of the contaminant in a groundwater reservoir is largely dependent on the amount 

transferred between the liquid and solid phases. 

Perhaps the most realistic macroscopic appraisal 

available is in the results of reaction-tower experiments. 

The analytical difficulties that appear whenever the 

adsorption isotherm is nonlinear are formidable, 

however, not insurmountable.

MATHEMATICAL TREATMENT OF THE 

DISPERSION EQUATION

The analytical treatment of dispersion phenomena 

within a porous matrix, based on the preceding theoretical development, presupposes that the dispersion 

coefficient is a known parameter. Thus, the full definition of the physical system depends on accumulation 

of data concerning this specific parameter. It will 

become apparent in the following discussion that only 

the simpler dispersion systems are treated primarily 

because of the need for laboratory confirmation. The 

discussion is, therefore, not intended to be a complete 

presentation of all known solutions of the differential 

equation. Also, this section presents methodology 

considered useful in the analysis of the dispersion 

equation.

The analytical methods, namely, Laplace transforms, 

the product solution, and the point source solution are 

the principal tools used in the study of the differential 

equation. These techniques are discussed to some 

extent in the work of Carslaw and Jaeger (1959) and in 

many texts in advanced calculus.

All the examples to follow are limited to flow in a 

semi-infinite medium. The use of more complex boundaries adds little to the study of the fundamental 

characteristics of the mass transport system because 

of the complexity of the final expression. The principal 

objective is to investigate initially the nature of dispersion. Hence, the convective term div ~^C will 

be assumed to be given by u dC/dx-, in other words, 

the flow field is assumed constant and unidirectional. 

Space dependency of the dispersion coefficient is not 

considered; variation is limited to directional dependency only. In systems involving an adsorptive solid 

medium, the problem becomes highly complex and 

thus no analysis has yet been attempted for anything 

more complex than one-dimensional mass transport.

The examples given are intended only as samples of 

solutions available and may not be physically realistic. 

This is especially true in examples 3 and 4 where a 

system involving a disk source was chosen to depict a 

two-dimensional case. Example 3 actually may not 

be physically realistic; however, it may possibly occur 

when the velocity u is small. That is, example 3 implies 

that a concentration gradient occurs at z=0.

Numerous solutions available in the published literature are presented in the following pages. Whenever 

possible, a reference is given so that the complete 

development can be consulted by those who are 

interested. In addition a limited number of functions, 

which the author arranged to have computed, are 

presented in graphical form as appropriate.

The first two examples considered are for dispersion 

within.a semi-infinite medium in a unidirectional flow 

field. It is postulated that the concentration is specified 

at z=0 and at the initial time f=0. The system is 

defined schematically in figure 4. For the specified field 

conditions, the field equation (eq 6), omitting the 

subscript x, becomes simply

D n 6 5-5-='^ 2 0 dC.dC 5 T-5T'

dx2 dx ot

(21)

The subscript x has been dropped in equation 21, since 

dispersion occurs only in the x direction.

Example 1. The plane a*=0 is maintained at a concentration C= CQ .

The boundary conditions are given by

0z0=0, x>Q

0(00,0=0,

The method best suited is the Laplace transfonnat

u average areal flow

C specified at x = 0

FIGURE 4. Definition sketch of one-dimensional system.
THEORY OF DISPERSION IN A GRANULAR MEDIUM 113

tion; substitution of the transform leads to the 

subsidiary system given by the equation

dC. -

The solution of this subsidiary equation subjected 

to the restraints due to the boundary condition is

ux

The inverse transform as given by Erdelyi and others 

(1954,p.246) is readily determined as

/ot N

(21a)

Example 2. The region z<0, initially at concentration

The boundary conditions of interest are specified

by

C(x, 0) = C0, z<0, 

C(x, 0)=0,x

Again, various analytical methods are available; 

the simplest and most straightforward method 

involves the use of a plane-source solution as the 

fundamental solution. The concentration at any 

point in space (see, for example, Eifai and others, 

1956, p. 54) is given by the expression

Summing over all possible values of £, gives

= erfc

where £=x ut.

(21b)

The graphical representation of the solutions given 

for examples 1 and 2 appears in figures 5 and 6. Note 

that because a constant concentration is maintained at 

x==0, example 1 exhibits an asymmetrical concentration distribution about x ut. However, it has been 

shown that the conditions of flow through porous 

media are such that the second term in equation 21 a is 

negligible (Ogata and Banks, 1961). Hence, the two 

systems exhibit no difference at some distance from the 

source.

370-837 O .70 3

Consider now, the system that is described within a 

two-dimensional frame of reference because of symmetry about a given axis. The definition plane is given 

in figure 4. Because of the mathematical complexities, 

this system is illustrated using circular surface sources 

maintained at constant strength and constant concentration. The results are sufficiently complex that in a 

discussion on evaluation of the dispersion coefficient 

(p. 130-132) an approximation is used to simplify the 

computations.

Assuming again a unidirectional ground-water flow 

field with parameters Dx and u oriented in the x direction and DT in the r direction normal to x, the field 

equation for the radially symmetrical system becomes

Equation 22 is equation 9 when symmetry exist about 

the axis x; hi other words, the variation in concentration with respect to the angle 6 is zero.

Example 3. A disk source, 0<r<a located at x=Q 

emitting at constant rate F=FQ , 

The boundary conditions are

 Fo at x=Q, Q<r<a, t> 

=Q at z=0,

Initially, the fluid in the porous medium is at zero

concentration.

The solution of this system is readily determined 

by use of the point source. Using the method described by Carslaw and Jaeger (1959, p. 266), the 

expression applicable for an instantaneous disk 

source is

-^) f" e- 

*Ux t/ Jo

where £=x ut. For the special case DT =DZ=D and 

u=0, the expression is the same as equation 9, page 

260, in Carslaw and Jaeger (1959).

To obtain an expression for the continuous disk 

source, the above equation is integrated between the 

limits 0 to some time t. After some algebraic manipulation the integral equation is written

C^ f

 f o Un Jo

where

u ux 4DXDT
114 FLUID MOVEMENT IN EARTH MATERIALS

0.001

FIGURE 5. Concentration distribution; plane source at i=0 maintained at constant concentration.
THEORY OF DISPERSION IN A GRANULAR MEDIUM 115

-2.0

(x-ut) 

FIQUEE 6. Concentration distribution; region x<0 initially at constant concentration.

The above expression may be written in two 

alternative ways, or

C

where

and

l=_ *1/2Jof e

-\dr

erfc|j--)9Vl+72X2l

-exp

For the steady-state condition, that is, letting 

t >°° in the above expression, the equation obtained 

is

The steady-state solution, although simpler, cannot 

be integrated and must be evaluated numerically. 

Along the axis r=0, the integral equation can be 

expressed in terms of an elementary function or,

uC ,

Example 4- A disk surface source 0<r<a located 

at z=0 and maintained at a constant concentration 

C0.

This example considers the steady state only since 

the solution is not available for the transient system. 

The boundary conditions for this specific problem 

are

C(Q, r)=0, r>a 

C(a>,r)=0, r>0. 

C(x, co)=0, z>0.

To analyze this problem, it is assumed that the 

solution is represented by the product of two functions of single variables, that is,

C=R(r)X(x),
116 FLUID MOVEMENT IN EARTH MATERIALS

where R(r) is a function of r only, and X(x) is a 

function of x only. Substitution for C in equation 22, 

noting that for steady state dC/dt=Q, the partial 

differential equation is reduced to the subsidiary 

equations

where f2 is an unspecified constant.

The solutions of the subsidiary equation are readily 

obtained. The solutions are

)=_42 exp [2a

Because of the condition specifying C be zero as 

x »oo, the positive sign before the radical in the 

second equation is discarded. The solution is thus

=r exp

The method of determining ^l(f) is described fully 

by Sneddon (1951, p. 59), the result is

-£=a f°

ko Jo

exp [2a (lA similar solution for DT=DX D and u=0 is

given by Carslaw and Jaeger (1959, p. 214).

The steady-state solutions to examples 3 and 4 are 

almost identical, except for the factor l/yT+VX2- The 

two solutions, presented as examples 3 and 4, are represented .graphically in figures 7 and 8. The numerical 

computations for the integrals in examples 3 and 4 were 

accomplished by integrating under a curve with a polar 

planimeter. In both examples the function could not be 

integrated into a readily calculable form. However, the 

integrals can be evaluated by numerical methods for 

various values of rj=r/a, which may be used in evaluating experimental data.

The previous examples considered dispersion in flow 

systems described by a single-space variable and by 

two-space variables. Examples 3 and 4 considered a 

radially symmetrical system with flow occurring along 

the cylindrical axis. Consider now, a line sink or a line 

source at r=0 in a cylindrical reference frame. This is 

generally termed "a radial flow system." The definition 

of the geometrical configuration is presented in figure 9, 

and the applicable differential equation follows.

Writing the expression for the conservation of mass, 

with reference to the definition sketch (fig. 9), the 

differential equation is

d /n ac\.z>ao do dc - D -- - --=-' d /n a- ( D

dr\

5-

or/ or or dt

,OQ ,

(23)

i . i i i i i i i ___ i ___ i ___ i ___ i _____ i ___ i

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

FIQUBE 7. Steady-state concentration distribution; circular surface source of constant strength.
THEORY OF DISPERSION IN A GRANULAR MEDIUM 117

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

FIGURE 8. Steady-state concentration distribution; circular surface source maintained at constant concentration.

Here, the concentration and the dispersion coefficient 

do not vary with the angle 8 nor with the distance z 

which is normal to r. Equation 23 is written to allow for 

a possible space dependency of the dispersion coefficient 

D in the direction r. The basis of this space dependency 

comes from experimental and statistical findings that 

D is directly proportional to the velocity of the fluid 

stream. Since it is readily seen that the velocity is

FIGURE 9. Definition sketch of radial flow.

inversely proportional to the distance from the source 

at f=0, this same relationship is assumed to hold for 

the dispersion coefficient.

The expression for fluid velocity may be formulated 

provided the quantity of flow or flux is specified at some 

value of r=r0) generally near the source. Since flow is 

assumed steady, continuity requires that the velocity 

at any point is given by the expression u u<,rQfr, where 

u0 is the average velocity at r r0 . Substituting the 

variable velocity into equation 23 leads to the expression

dr

_ 

dr~ dt (24)

Experimental evidence indicates that the dispersion 

coefficient D in a homogeneous medium varies very 

nearly with the first power of velocity; thus, D can be 

written

where Dm is a constant related to the solid and liquid 

properties. Substituting the preceding expression into 

equation 24 gives

6V V

D/dC

r2 dr

l_dC} 

DQ dt'
118 FLUID MOVEMENT IN EARTH MATERIALS

where a=u0r0 . Accordingly, collecting the terms leads 

to the equation

(25) _a_ dC=j^ dC 

6V2 D0 6V DQ dt'

Example 5. Dispersion in a fluid flowing radially from 

a line source at r=0.

The boundary conditions to be satisfied are

tf(r,0)=0, r>rQ 

#(00, 0=0, t>0.

The method of analysis because of convenience and 

simplicity in determining the coefficient of the solution is the Laplace transform. Equation 25 transforms to

_ == 

dr2 D0 dr D0 '

The solution of the transformed problem is

where z=r a/4D0p, z0=r0 a/4D0p and K^(x) is 

the modified Bessel function of the second kind of 

order %.

The inversion of the Laplace transformed system 

has been obtained by the use of the hi version theorem 

(Ogata, 1958, p. 88). The final expression is given by

where

and JK(X), YX(X) are Bessel functions of the first and

second kind of order %, respectively.

The solution hi its present form has not been evaluated; however, a finite difference scheme has been applied to the differential equation. The integral in its 

present form is cumbersome to evaluate because of the 

appearance of the Bessel function in both numerator

and denominator. Shamir and Harleman (1967) have 

obtained an approximate solution of equation 25, 

numerically. Experimental analysis of this problem 

has been done by Lau, Kaufman, and Todd (1959); 

however, certain fundamental simplifications were required for the comparison of experimental with theoretical results. A result of the finite difference computation for the value D0/u0r0=l'XlQ~2 is presented in 

figure 10.

The solutions for dispersion hi one-dimensional, 

two-dimensional, and the special case of radial-flow 

systems have been presented. Continuation of this 

analytical trend would suggest that dispersion in a 

three-dimensional system should be considered next. 

However, the relative difficulties in coping with a 

three-dimensional ground-water system renders the 

attempt hardly worthwhile. Thus, it seems more useful 

in the final two examples of dispersion without chemical 

reaction to consider approximate solutions for twodirectional dispersion. In the two-dimensional field 

equations (eq 22) it is assumed that, if the term 

Dxd?C/dx2 is small compared to the other terms, the system can be described as consisting of flow and dispersion 

that are orthogonal to each other. This is signified in 

the definition sketch presented as figure 11.

For dispersion process from a disk source, the field 

equation is

£ (26)

For the plane source hi which y>0 initially at a constant concentration, the appropriate field equation

is

dC dC dx' (27)

Equation 27 is identical with the diffusion through a wall, 

and equation 26 is identical with diffusion from a 

cylinder moving hi the direction x at the rate u.

Example 6. Approximate radially symmetrical model 

of lateral dispersion.

The approximation of lateral dispersion, equation 

27, in which the justification of its use is given, was 

treated by Ogata (1961). In previous discussion on 

transformation of the coordinates, it was pointed 

out that the coordinate transformation £=x utm&y 

be used to eliminate the convective term. Equation 

26 may thus be written

^VoWrdr/ dt

Although C depends on x, this need not be specified 

in the boundary conditions since x does not appear
THEORY OF DISPERSION IN A GRANULAR MEDIUM 119

1.0

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

EXPLANATION

1 (/ -/ )= 1.75 feet

2 (r-y0 )=2.75feet

3 ( f-f ) = 3.50 feet

4 ( f-r0 ) = 4.00 feet

0.68 0.76 0.84 0.92 1.00 

<W/ r(r-r0 )

1.08 1.16 1.24

FIGURE 10. Concentration distribution for radial-flow system.

explicitly in the differential equation. Hence, the 

boundary conditions are stated

C(r,0)=0,r>a, 

(d<7/dr)r=0=0, t>o.

The solution of the above differential equation for 

the specified boundary condition is obtained by 

Laplace transformation. The subsidiary equations 

for r>a and r<a are

cPX , 1 d\ ,-

f r>a dr* r dr * ' ' '

- C where \=C-\ - and g*=p/D. The two regions are

described by identical equations when X is used.

Particular solutions valid for the region of interest 

with the specified constraints are

-Q fl 

PL

^(qa^gr) "1

c=^ /i(go)g0(gr)

p

The inverse of the above expression is given by 

Carslaw and Jaeger (1959, p. 346). The solution for 

all values of r is given by

c r° 7^=a 

ko Jo

exp ( 

From equation 27, it is immediately apparent that 

for steady-state conditions the solution can be obtained

T»

by replacing t with -. The steady-state conditions are 

thus computed from the same expression.

Example 7. Approximate steady-state model of lateral 

dispersion in a unidirectional flow field.
120 FLUID MOVEMENT IN EARTH MATERIALS

r or y

Fr = Drb C/6r

FIGURE 11. Definition sketch of system depicting no dispersion in X direction.

The model considered is dispersion in a fluid moving 

in the direction x with the plane ?/=0 initially at a 

constant concentration C0 . This special case was 

treated by Harleman and Rumer (1962, p. 24). 

Equation 27 is the appropriate differential equation 

and the boundary conditions are

,y)=0, <=°<2/<0 

=0, x>QThe application of Laplace transform- in equation 

27 leads to the subsidiary equation

where C'=C/C0 .

The solution of this subsidiary equation is readily 

determined as

Hence, ushig the complex inversion theorem, the 

solution obtained is,

The expression for example 6 has been computed

numerically for various values of

a'u 

4DTx

The expression

for example 7, on the other hand, is a standard tabulated function; hence, numerical values are available. 

These two expressions are presented graphically as 

figures 12 and 13. Because of the simpler nature of 

these expressions, as compared with those obtained in 

examples 3 and 4, their use, whenever applicable, is 

favored for computing-numerical values of the dispersion coefficient for any given experiment. Figure 14

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1-1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0

FIGURE 12. Approximate concentration profile due to lateral dispersion; radially symmetrical model.
THEORY OF DISPERSION IN A GRANULAR MEDIUM 121

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1,5 1.6 1.7

FIGURE 13. Approximate concentration profile due to lateral dispersion; two-dimensional model.

shows the curves for examples 3, 4, and 6 for two 

values of 4Drx/ua?. As expected for higher values of 

iDrX/ua2 , the difference between the three functions 

become appreciable. For a first approximation, however, the use of the solution obtained in examples 6 

or 7 would be simpler and of sufficient accuracy.

The following examples of mathematical analysis of 

dispersion consider an adsorbing medium. The system 

will be limited to a semi-infinite medium. The adsorption rate is assumed to be linear or more specifically 

of the form dS/dt=b(C mS). For flow in a medium 

with finite boundaries, the reader is referred to several 

examples given by Crank (1956) as the complete solutions are too complex to be reproduced here.

When mass transfer due to adsorption is superimposed 

on the convective and dispersive transport, the continuity consideration in a one-dimensional flow system 

leads to the differential equation

(28)

The tune rate of change in the solid-phase concentration 

was previously indicated to be approximated by the 

linear equation

~=b(C-mS). (29)

In the examples to follow, equations 28 and 29, with 

suitable modifications, are used to describe the concentration distribution within an adsorbing medium. The 

flow occurs in a semi-infinite medium; hence, for all 

following examples the boundary conditions are

C(Q,t) =<?<>, *>0, 

<7(z,0) =0, x>0,

0(00, £)=0, t>Q,

S(z,0) =0, z>0.

Example 8. Dispersion in an adsorbing medium; 

adsorption rate given by dS/dt=bC, (m=O).

The solution of simultaneous equations 28 and 29 

(with m=0) for the prescribed conditions can be 

determined by assuming that it may be written in the 

form

C(x,t)=T(x,t) exp -

Substituting the above expression into equation 33 

leads to an equation with T as the dependent variable, 

or
122 FLUID MOVEMENT IN EARTH MATERIALS

exp(Dr x\2/ua2 ) F(\)d\

4Dr x/ua =2.0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7

FIGURE 14. Comparison of examples 3, 4, and 6.

For the function T, the condition at z=0 is given by

The general solution of the equation above is given 

by Carslaw and Jaeger (1959, p. 63). For this specific 

case the solution is

=o exp

The integral above is a well-tabulated function called 

the error function, or

-= exp [(l-M)ux/2D] erfc lwhere

Example 9. Dispersion in an adsorbing medium; 

adsorption rate given by dS/dt=k(C mS) with 

^0 and dC/dt«dS/dt.

exp [(l+M)ux/2D] erfc

The differential equation for the specified condition 

is reduced to udC/dx dS/dt. Using the Laplace 

transform, the subsidiary equations obtained are

u<fc_=_ ^

dx * 

pS=bC-bmS.

The solution satisfying these simultaneous equations 

is

C 1 TT=- exp

Co p

L r~ bx

\p+km/_\

The inverse of the above expression is available in 

tables of the inverse transform. The solution obtained by Hougen and Marshall (1947) is

Q

7~-=l exp (

fb

 bmt) I 

Jo

exp

Example 10. Dispersion in an adsorbing medium; 

adsorption rate given by dS/dt=b(C mS), with 

dS/dt»dC/dt.

The postulated conditions reflect primarily the 

mathematical viewpoint of looking toward a simplified equation. Again, the use of the Laplace integral
THEORY OF DISPERSION IN A GRANULAR MEDIUM 123

transform reduces the partial differential equations 

to the subsidiary expressions

dC 

p'S=b(C-mN}.

The solution of these simultaneous equations, 

obtained by classical methods, is

6XP

The inverse is obtained by standard means and is 

described by Ogata (1964b). The expression obtained

is

O 2i& I

Co -y/TT JO 

 \TftJ 0

exp (£2 -<

exp -£-

where a=ux/2~D, p=Dba2/u2 , and 3(x,y) is the 

Goldstein J function which has been defined in 

example 3.

Example 11. Dispersion in an adsorbing medium; 

adsorption rate given by dS/dt=b(C mS), with

The equation considered is again a vjariation of 

the complete field equations 28 and 29 with the 

condition that dispersion in the direction of flow is 

small. The method of analysis again involves application of the Laplace integral transform. The 

subsidiary equations for this case are

u

pS=b(C-mS).

The solution of these simultaneous equations 

can be determined by elementary means, and the 

resulting expression is

JLl

Co p

f f» 

L (p+bm) 

The solution as given by Ogata (1964b) is

a r / r\~\ rte/w r / 7 r\ I £=l-exp -bmtt--) e-*I0 \ Jbm(t--)\ \d\ Co *L \ u/JJo U L\ \ uj J

Note that the preceding solution is the same as that 

for example 9 when (t x/u) is replaced by t.

Example 12. Dispersion in an adsorbing medium; 

adsorption rate given by dS/dt=b(C-mS).

This final example employs equations 29 and 30. 

The subsidiary equations, obtained by use of the 

Laplace transform, are

pS=b(C-mN).

The particular solution of the transformed equations 

for the specified boundary conditions is given by the 

 the expression,

The process of inversion is described by Ogata 

(19646), and the final expression is written

Jl

-= |_,_ exp (-r2V/W

where c=6 72/£2, T=bm(t-y2/?), a=ux/2D, y2=x2/4D

and J(c, T) is the Goldstein J function as defined

previously.

Note that this solution is identical with the solution 

in example 10 when the term (t y2/^2) is replaced by t.

Because the development of the solutions for examples 10-12 is relatively recent, their numerical evaluation has not been completed. The solution for example

8 can be readily computed, whereas the solution for

9 has been used in analyzing reactor columns. Figures 

15 and 16 are graphical representations of the expressions for these two examples.

The number of conditions, other than the onedimensional system, that have been investigated are 

relatively few. Because of the similarity between equations describing dispersion in ground water and those 

describing diffusion and heat-conduction systems, 

numerous analytical solutions are available. And with 

the increased use of digital computers, a variety of more 

complex mathematical problems can now be solved 

numerically.

NATURE OF THE DISPERSION COEFFICIENT

In the analytical development of the macroscopic 

theory of dispersion it is assumed that the characteristics 

of the dispersion coefficient are known, and that it is 

virtually a constant. In reality, the nature of this 

coefficient is unknown and is the subject of considerable
124 FLUID MOVEMENT IN EARTH MATERIALS

0.5 10.0

FIOUBE 15. Concentration profile of dispersion in an adsorbing medium with dS/dt=bC.

investigation. Extension of the theory describing dispersion in a porous medium, or the ability to use 

dispersion equations for systems other than isotropy, is 

entirely dependent on the evaluation of this coefficient 

by either analytical or laboratory means. Further, this 

critical evaluation depends entirely on the ability to 

describe analytically the physical properties of the 

porous medium. This stage of investigation has not been 

reached as yet; thus, the discussion to follow is limited 

to evaluation of the dispersion process in isotropic 

material.

Many experimental and analytical papers published 

in journals show the dependency of the magnitude of 

the dispersion coefficient on the flow rate and the effective gram size diameter of the porous medium. This 

dependency comes about because, provided the sand is 

uniform, the geometric characteristics are described by 

a single parameter the sand-grain diameter. But, 

whenever the materials that make up the matrix depart 

from isotropic conditions, the representation of the 

medium by a so-called characteristic length becomes a 

difficult task, as experienced in attempts to generalize 

the concept of a friction factor in fluid flow through 

porous media (Carman, 1956). Generalization of the 

results of laboratory and analytical studies cannot be

realized unless a porous matrix can be described in its 

generality, including parameters such as the grain 

distribution and degree of cementation. This need is 

reflected in all problems involving the granular medium 

and in some quarters gave rise to such correction terms 

as "tortuosity." However, the only route opened at 

present is to assume that the Reynolds number establishes flow similitude and the Schmidt number establishes similitude in the diffusion. Then, the correlation 

of these two similarity criteria makes it possible to write 

an analytical expression for the functional relationship 

between the flow and diffusion processes.

Characteristics of the dispersion coefficient for conditions encountered in the laboratory are generally 

studied through statistical evaluation, or by correlation 

of the experimental data with results of a mathematical 

model. The statistical model is important in that to a 

limited extent it does point out the parameters that 

control the magnitude of the coefficient. Its weakness 

lies in the need to reduce the porous matrix to an extremely simplified structure in order to render the 

problem mathematically tractable. The use of the 

macroscopic model is, however, a little more difficult 

since the nature of the coefficient is only specified as a 

constant. Thus, any correlation of these coefficients with
THEORY OF DISPERSION IN A GRANULAR MEDIUM 125

FIQTIBK 16.- Concentration profile of dispersion in an adsorbing medium with dS/df ^(C-mS), and T^'r 

various parameters must be accomplished by methods 

of regression, unless similarity can be established between the various flow conditions in the laboratory and 

the field.

Numerous experiments indicate that the use of the 

Reynolds number and Schmidt number are useful in 

establishing results comparable to those obtained by 

statistical analysis, provided the porous matrix is isotropic. However, in natural soils, or when air entrapment occurs in flows, the dispersion coefficient cannot 

be readily correlated to an effective diameter or a 

characteristic length. Thus, the expression for the D is 

written Dt=iDmu, where Dm is an experimental constant 

that depends on characteristics of the porous medium.

The statistical development of the dispersion theory 

is of interest because it reveals the dependence of the 

magnitude of the dispersion coefficient on the characteristics of the flow and the porous medium. Stringent 

conditions, however, must be postulated to simplify the 

model so that the mathematical concept can be utilized. 

As is generally true in most mathematical development, 

the mathematical model adopted has little resemblance

to the actual porous matrix. The basic theory of the 

statistical analysis is described in four papers published 

in separate scientific journals. These papers were written 

by Scheidegger (1954), Jong (1958), and Saffman (1959, 

1960). Saffman's and Jong's models are similar in 

nature; however, Scheidegger, who introduced the concept of statistical hydrodynamics, chose a mathematical model that assumed random motion of the fluid 

particles.

The results of both Jong and Saffman are reproduced. 

The porous medium model of Jong is reproduced as 

figure 17. The pore channels are assumed to be capillary 

tubes multiply connected at some end points and uniformly distributed in all directions in space. The medium 

is assumed to be made up of spherical particles that 

form polyhedral cavities connected by triangular capillary tubes. The velocity within a tube varies according 

to the position of the fluid particle within the canal; 

however, for simplicity it is postulated that the mean 

velocity is represented by the rate at which the fluid 

particle containing, say, a contaminant or tracer, passes 

from point to point. The multiple connections at the 

tube-end points represent a limited choice of direction 

confronting an arriving particle of fluid containing the 

contaminant or tracer. This postulate, to some degree,
126 FLUID MOVEMENT IN EARTH MATERIALS

(From Jong, 1958)

FIGURE 17. Random path chosen by a fluid particle moving 

through the canal system.

indicates randomness of the particle path, and together 

with the fact that the pressure gradient is proportional 

to cos 6 (fig. 17), makes it possible to trace or determine 

the probable location of the fluid particle at any given 

time.

The complete analysis presented by Jong (1958) is not 

repeated here inasmuch as it is basically a justification 

for the use of a given probability function. Of immediate 

interest, however, is the end result which follows and 

indicates the functional dependency of the dispersion 

coefficient on the flow field and the pore characteristics. 

The use of average velocity to characterize the flow 

makes it possible to characterize the pore system by a 

length parameter defined here as the elementary

length I. The final results of Jong's analysis are the 

statements that the dispersion coefficients are related 

to the flow and media parameters according to the 

expressions

DT =Bu Z/16 

Dx=u I -ln 7)/6. (30)

In equation 30, the symbols are defined as follows:

u= average velocity,

/=length of elementary canal,

A =f unction of Z0 where Z0 is the distance 

along the z direction at which the maximum number of particles are located at 

time T0 , 

In 7=Eulers constant = 0.577.

The elementary length of the porous medium, and 

the residence time of the tagged fluid therein, constitute 

the variables of the statistical analysis. Although the 

elementary length may be used to describe the medium, 

the magnitudes of the dispersion coefficient cannot be 

computed directly unless there is a relating measurable 

parameter associated with the characteristic length. 

Jong reports that for his experimental model, in which 

the medium was composed of glass beads, the calculated 

value of I was one-third of the mean grain diameter. 

For conditions other than isotropy, it is expected that 

the value of I will depend on the grain-size distribution 

and the mode of packing.

Staffman (1959), in his first paper on the theory of 

dispersion, used a statistical model identical with that 

of Jong's. The basis of his analysis is the visualization 

of the fluid path as a random walk process where the 

length, the direction, and the duration of each step 

are random variables. The expressions for radial and 

axial dispersion derived from this mathematical model 

are

(31)

where

or

The symbols in equation 31 and the related expressions 

are defined as follows:
THEORY OF DISPERSION IN A GRANULAR MEDIUM 127

£=time, 

<o=estimate of time for appreciable diffusion to

take place along the pore, 

u=average velocity in the x direction,

Z=pore length, and 

n=number of steps in the random walk analysis.

Equations 30 and 31 indicate a slightly different result 

for Dx , but the radial diffusion coefficients DT derived 

by these two methods are identical

Saffman (1960), in a later paper, investigated dispersion phenomena by using a Lagrangian correlation 

function in an attempt to compute values of dispersion 

coefficients for a more generalized condition in terms 

of the dimensionless number u l/D, where D is the 

molecular diffusivity. The results of the analysis are 

much more complex, and the expressions are not 

reproduced here. Saffman states that comparisons of 

these results with experimental analysis show a large 

discrepancy, probably attributable to nonuniform 

packing of the granular material within the experimental column and possible channeling along the 

boundary between container and porous material. Or 

this discrepancy may be the consequence of oversimplifying the existing mathematical model in order to 

render the analysis tractable. The conclusion that can 

be drawn is that too much weight is being placed on 

the single parameter I the elementary canal length. 

Differing somewhat from Jong, Saffman calculated the 

value of the coefficient on the assumption that the 

magnitude of I equals the mean diameter of the sand 

grains.

In an attempt to reconcile some of the discrepancies 

noted, especially transverse to the flow direction, 

Saffman (1960) pointed out that in a strictly twodimensional system, as depicted by the mathematical 

model, two of the basic assumptions may not represent 

good approximations. Physically, because of the prevailing steady-state laminar flow, there can be no 

intertwining of the streamlines except in a threedimensional flow system. Hence, the variation of the 

width of the stream tubes is on the order of the diameter 

of the sand grain or less. Thus, the assumptions of 

statistical independence of each streamline, and the 

uncorrelated nature of the mean displacement transverse to the flow, cannot be strictly valid. Using 

statistical models these investigations find that at 

some time after the introduction of the contaminant, 

or at some distance from the contaminating source, 

dispersion is described by equation 6. This conclusion 

based on the assumption that the central limit theorem 

is applicable was first stated by Beran (1957).

It is apparent that the choice of the mathematical 

model in the statistical analysis is limited entirely by

the mathematical difficulties, and thus cannot in any 

way represent a true porous medium. The departure 

from the physical features of the actual granular material may be the primary reason for the poor correlation 

between calculated and observed values. These difficulties are inherent because of limited knowledge in 

the statistics of sedimentation and other geologic 

processes. This differs to a large degree from other 

transport systems in which the kinetic theory of gases 

or the statistical mechanics is applicable. In systems 

where these theories are applicable, meaningful statistical averages can be derived from the mathematical 

model, making possible the computation of the relationship between the coefficient and the various parameters. Certain limitations are imposed by the chosen 

mathematical model, but they are not so severe as in 

the analysis of flow through a porous matrix.

The generalization of the expressions for the dispersion coefficient cannot be readily achieved at present. 

However, the statistical studies may be used as a 

guide in an attempt to formulate empirical expressions. 

Articles such as that of Perkins and Johnston (1963) 

give a number of empirical expressions that have been 

used to describe systems of specific interest.

EVALUATION OF DISPERSION COEFFICIENT BY 

EXPERIMENTAL METHODS

Analysis in the preceding section indicates that the 

magnitude of the coefficient is dependent on the 

velocity of fluid flow and the elementary canal length 

used to describe the porous material. The validity of the 

dispersion equation is also established, although no 

proof was given. The magnitude of the dispersion coefficient is obtained by correlation of the solution of 

equation 6 with experimental data. This section discusses the use of the dispersion equation as the tool in 

the experimental determination of the magnitude of the 

dispersion coefficient. The difference between the two 

methods is that the statistical method requires expermental determination of the medium characteristic, 

whereas the macroscopic method requires evaluation of 

the coefficient itself.

Many solutions have been presented in a previous 

section to indicate the nature and extent of theoretical 

investigations of transport phenomena. These solutions 

were obtained for boundary conditions which involved 

semi-infinite porous media, readily amenable to representation by laboratory test models. The verification of 

theoretical results has generally been confined to the 

analysis of one-dimensional models that are easier to 

build and test. These models are most appropriate for 

experimentally establishing the magnitude of the longitudinal dispersion coefficient. Accordingly, the quanti-
128 FLUID MOVEMENT IN EARTH MATERIALS

fication of the coefficient is based on the conditions 

specified by

Urn G(x,t)=Q,t>0. 

X »oo

The first statement indicates that over a defined plane, 

say z=0, the concentration is to be maintained at a 

known and feasible value. The second condition is 

merely a statement of the concentration of the fluid 

which initially fills the porous matrix. This condition 

can be established to be either zero or any known value <70 . 

For the conditions depicted, the solution of the differential equation is as given in example 1. Most studies 

indicate that the second term in that equation is generally small, so that the solution presented as example 2 

gives a good approximation (Ogata and Banks, 1961). 

Thus, for purpose of correlation, the expression to be 

used (see eq. 26b) is

_ .e ft 2 erfc

where T=ut/x and X=Dfux.

The laboratory model constructed to conform with 

the specified boundary conditions consists of a tubular 

column filled with unconsolidated sand. Except for 

the contact at the wall of the tube, where some flow 

channeling may occur, the model satisfies the conditions specified for a semi-infinite porous medium. 

The tracer chosen depends on the nature of the investigation. Common salt has proven the easiest to 

measure and is used whenever density is not a critical 

factor. Because of its electrolytic characteristic, in-place 

measurement of the salt concentration can be obtained 

by means of an electric probe. At points where large 

enough samples of fluid may be taken, however, chemical methods represent an economical and rapid means 

of determining the concentration at any specified tune. 

In either event the requirement is that measurements 

be taken at some distance from the tracer source, at a 

location and in a manner which does not disturb the 

flow. A typical setup for a columnar model is presented 

in fibure 18.

Since the solution is available, the magnitude of D 

can be rapidly computed. The simplest means of calculating this result is to establish the value of the ratio 

X=8/<t>, where 0=(l-T)/2 JXT and 6=(l-T)/2 

. For any value of C/C0, <t> can be determined from

the tables of the error function. By definition, 

8=i/X<l>; thus, since T=ut/x and T is known for any

L-B1

Screen 

6

 =P

1 Tracer reservoir 

2 Solute mixing and constant 

head tank 

3 Tap water constant head tank 

4 Overflow tanks 

5 Petcocks 

6 Sand column 

7 Downstream overflow 

8 Overflow control valve 

9 Sampler tubes

7 

-^ -Perforated steel plate

9 r^

- O-~ Vnlves

FIGURE 18. Typical setup for column experiment.
THEORY OF DISPERSION IN A GRANULAR MEDIUM 129

given series of tests, the values of 6 and <f> are available. 

The values of these two parameters are plotted on 

rectangular coordinates that associate 0 and <t> with a 

given value of C/C0 . The resulting plot can be analyzed 

by a straight-line correlation that relates the <£ and 6 

values. The slope of the straight line is readily determined and by noting that the average slope is 

I) QI^=JX = -\I , the value of D can be computed. V ux

This method affords a simple means for computing 

the value of the dispersion coefficient directly from 

any given series of tests.

The results obtained from any such analysis, however, 

yield simply a number applicable to a particular series 

of tests. A generalized expression for the results of all 

such tests cannot be realized unless relating parameters 

that hold for all conditions are identified. One method 

for obtaining these parameters is by rewriting equation 

6 in dimensionless form. For this purpose let,

x=x/d, r=ut/d, and C'=C/C0

where d is defined simply as a characteristic length. 

Substituting the dimensionless variables in equation 

22 results in the expression

D dC' ^

ud 3x2 dr

Accordingly for the postulated conditions, the generalization of the test results depend on the single parameter 

D/ud. The rewriting of the field equation in its dimensionless form permits identification of a parameter 

which, in effect, is the similarity criterion for the 

process described by the equation.

The resulting parameter, however, is not a familiar 

one in fluid mechanics literature. To resolve this, consider the flow mechanism in a porous matrix. Dimensional analysis like that described by Muskat (1946, p. 

56) indicates that the flow is characterized by a parameter known as a Reynolds number, which appears frequently in the hydraulics. Since flow in a porous 

medium is virtually flow in capillary tubes, it is not 

surprising that the Reynolds number describes the 

similarity conditions. Assuming the medium is described by a characteristic length d, the Reynolds number may be written, ud/p, where d for a homogeneous 

medium is either the pore channel diameter or the grain 

diameter and v is the kinematic viscosity of the fluid. 

Accordingly, the parameter Djud may be rewritten as

 /Re. The parameter D/v is the familiar dimensionless

group appearing in the study of isothermal diffusion 

and is called the Schmidt number (Bird, and others,

1960, p. 512). The coefficient of the differential equation 

in terms of these two dimensionless groups may be 

written

D/ud= l/(Sc) (Re) = 1/Pe,

where Sc is the Schmidt number; Re, the Reynolds 

number; and Pe, the Peclet number. Thus, based on the 

various assumptions made in the development of the 

differential equation, the dispersion process in a porous 

medium is characterized by the Schmidt number and 

the Reynolds number.

For purpose of generalization, the magnitudes of the 

dispersion coefficient resulting from computation of 

experimental observations are correlated by plotting 

the Schmidt number against the Reynolds number. 

An illustration of this correlation method is given as 

figure 19. The plot is on a log-log scale and produces a 

straight line having a slope of nearly unity over the 

range of Reynolds numbers less than one.

Another common practice in evaluating an exchange 

system is to correlate the Peclet number, udjD, with 

the Reynolds number. It can be seen that because 

Pe=(Sc)(Re), if similarity in the dispersion phase is 

represented by the Schmidt number, the correlation 

wDl result in a horizontal line. An illustration of a plot 

of the Peclet number and Reynolds number is given in 

figure 20.

These are the two methods of correlation most 

commonly used hi analyzing mass transport in a homogeneous granular porous medium. A majority of investigations indicate that a plot of a Reynolds number 

versus a Schmidt number (fig. 19) results in a straight 

line of unit slope, thus, indicating a direct relationship 

of the dispersion coefficient to the velocity and the 

characteristic length the average grain diameter or

Daud.

The works of Jong and Saffman, however, point out 

that the proportionality constant includes the residence 

time of the fluid; hence, the coefficient is generally 

written D=Dmu, where Dm is an unspecified constant. 

Scheidegger (1960, p. 259) further indicated that the 

relationship should be written D=Dmun, where the 

value n is either 1 or 2, depending on the velocity dis^ 

tribution taken to represent flow in the pore channel. 

In more recent analysis, HarJeman and Rumer (1962) 

obtained values of Dm =0.027 and n=l.W for glass 

beads of 0.39 mm average grain diameter, and Z)m=0.09 

and n=1.18 for glass beads of 0.965 rnm average grain 

diameter. Because the available data are limited, it is 

difficult to confirm or refute the validity of this type of 

expression. The limited data, coupled with ease of 

computation, popularize the use of n=l in analyzing 

dispersion in a homogeneous medium.
130 FLUID MOVEMENT IN EARTH MATERIALS

10

1.0

0.1

0.01

0.001

Sc =0.08/Re

EXPLANATION 

O Sand 1, d = 0.00371 ft 

O Sand 2, d = 0.00199 ft 

C Sand 3, d = 0.00135 ft

-O- Glass beads, d = 0.0116 ft 

(j) Glass beads, d = 0.00306 ft

Data fromt^ Dankwerts (1953)

- Day (1956) 

| Jong (1958)

-%-Kramers and Alberda (1953) 

J^f Banks (personal commun., 1958)

Sc --= 1 /Re

-o

I I I I I I___I I I I I I I I

0.01 0.1 1-° 

Re = ud/v

10 100 1000

FIGURE 19. Plot of the Schmidt number versus the Reynolds number.

In the section headed "Field equation for anisotropic 

dispersion," the directional properties of the dispersion 

coefficient were discussed. Laboratory observations and 

analytical studies, like those of Bear (1961), indicate 

that the magnitude of the coefficient is different in at 

least two directions, parallel to and transverse to the 

direction of flow.

The few available laboratory data for transverse 

diffusion are generalized by a plot of the Peclet number 

versus the Reynolds number. A general review of the 

published data relating to oil production is presented 

by Perkins and Johnston (1963). However, the specific 

area of interest, as displayed by the analytical and 

experimental data, is that related to flow conditions for 

which the Reynolds number >1. In a ground-water 

regime, this represents a range of flow that is not normally observed.

For the lower range of the Reynolds number, representative of ground-water flow conditions, the analysis 

uses the two-dimensional dispersion equation with the 

constant dispersion coefficient developed previously. 

The most practical flow system for studying dispersion, 

because of ready representation in the laboratory, is

the axially symmetrical system described by the 

equation

1 dC\ . d2C dC dC

Again, DT and Dx represent the radial and longitudinal 

dispersion coefficients, respectively, and the fluid flow 

is in the direction parallel to the x axis.

To determine the grouping of the parameters involved 

in the process of mass transport, the above equation is 

rewritten in dimensionless variables. To accomplish this, 

put x=xd, r=pd, and t r d/u where d is a representative or characteristic length describing the porous 

medium. Substitution of these dimensionless variables 

leads to the expression

_ _ Pe dp2 P dp dx2 ~

dC' 

dx

dCf 

dr

where PeT=ud/DT is the radial Peclet number and 

¥ex=ud/Dx is the axial Peclet number.

The Schmidt number, for directions r and x, and the 

Reynolds number characterize the dispersion and convection, respectively, for a radially symmetrical system.
THEORY OF DISPERSION IN A GRANULAR MEDIUM 131

1.0

-O-

-o

f

II

0.11 C

0.1 f)

0.08

o

O

O

O

EXPLANATION 

O Sand 1 

3 Sand 2 

 Sand 3

-O- Glass beads, large 

(J) Glass beads, small

Data from - 

^ Dankwerts (1953)

- - Kramers and Alberda (1953) 

P Banks (personal commun.,1958)

0.01 i i i i

1.0 1.0 10 100

Re = ud/v 

FIGURE 20. Plot of the Pedet number versus the Reynolds number.

Published articles, for example, Kramers and Alberda 

(1953), indicate that the longitudinal and radial Peclet 

numbers are constant. Therefore, both radial and 

longitudinal dispersion coefficients may be expressed 

explicitly by the relationship D=Dmu, where Dm is a 

constant for any given granular material.

Because of difficulties in building a two-dimensional 

physical model, the radial dispersion is investigated by 

considering a system where the longitudinal dispersion 

is small. One approximate analysis for steady-state 

condition is described by Ogata (1961) and may be 

used to calculate the magnitude of the radial dispersion 

coefficient from any suite of experimental data. This 

approximate equation is given in example 6. The complete solution, however, as previously stated, is too 

complex for rapid computation of the dispersion 

coefficient; hence, the expression

C_

CQ =1i ~exp

applicable for r=0 may be used.

Evaluation of the experiment conducted in Phoenix 

by the Geological Survey (Skibitzke and others,

1961) was based on the preceding expression. Because 

of some variation in the experimental data, only the 

average value of I>r=lX10~5 cm2/sec was determined 

to be significant. The velocity of the fluid in the experimental model was approximately 4X10"4 cm/sec, and 

the computed value of D compares well with results 

obtained by others and summarized by Perkins and 

Johnston (1963). It is also interesting to note that this 

value of the dispersion coefficient is near that of 

molecular diffusion of salt in water.

A second model of interest is the one used by Harleman and Rumer (1962). A two-dimensional dispersion 

model was analyzed both theoretically and experimentally. The mathematics of the mass transport 

system is described in example 7 where the solution was 

presented as

T^ofl+erf

Computed values of the dispersion coefficient, 

based on the observed and the approximate mathematical models, show a range of Dv from 9X10~5 to 

6.5X10"4 cm2/sec for values of the fluid velocity rang-
132 FLUID MOVEMENT IN EARTH MATERIALS

ing from 0.012 to 0.28 cm/sec. By correlation of the 

dispersion coefficient with the Reynolds number, the 

proposed expression for the dispersion coefficient 

transverse to the direction of flow was

Using this expression with the results of the Phoenix 

experiment, a value of 0.96 X10~5 cm2/sec was obtained 

for Dv, which shows good agreement.

SUMMARY

The theory of dispersion of miscible fluids has been 

developed extensively for fluid flow within a homogeneous porous medium. In the development, a fluid 

containing a given dissolved substance is described 

as moving according to Darcy's law, provided density 

differences between that fluid and the natural host 

fluid are negligible. However, the transportation of 

the dissolved substance has both a convection component, as described by Darcy, and a dispersion component. In all investigations the dispersion component 

is considered to be primarily an effect of the microscopic 

velocity variation within the pore channel and, hence,

is simply a perturbation from the described, average 

flow. There are other mechanisms involved hi the transport of the dissolved component, but their magnitude 

is small for normal ground-water flows, excepting 

adsorption. This mechanism of dispersion is treated to 

be analogous to eddy diffusion; thus, an expression 

similar to Fickian law of diffusion is used in the 

mathematical model to describe the dispersion 

component.

Evaluation of dispersion theory, based on the stipulated assumptions, entails laboratory investigations 

that are necessarily simple. Laboratory analysis of a 

linear flow system can be readily carried out and is 

especially useful in evaluating the one-dimensional progress of a slug of tracer. Such laboratory investigation 

is described in the section on determination of the 

dispersion coefficient. The results comparing theoretical and actual concentration distributions in the frontal 

zone of a linear flow system are typified in figure 21. 

For such a flow system there is ample evidence of this 

nature which indicates that the theory of dispersion 

hi a homogenous medium is a good first approximation.

On the basis of findings from experimental evaluation 

of the theory of dispersion, the application of the

Experimental - glass beads

-1.5 -2.0 

(From Ogata, 1958)

FIGURE 21. Example of comparison between theoretical and experimental results. From Ogata, 1958.
THEORY OF DISPERSION IN A GRANULAR MEDIUM 133

theory is progressing to more complex flow systems. 

A good example is the derivation of an analytical expression representing the mixing of two miscible fluids 

of different densities. This effort is directed particularly 

toward analyzing the conditions encountered in coastal 

aquifers where observations reveal the development of 

a diffused zone between the fresh- and salt-water zones. 

The theoretical treatment is basically that of a onedimensional dispersion system in which the coordinates 

are oriented with the position of the interface between 

the two fluids.

The field applicability of the analysis for a unidirectional transport system is limited. For problems involving any injection of contaminant, the dispersion 

phenomena take place in at least two-dimensional 

space. Accordingly, a description of the transport component perpendicular to the direction of flow is imperative. The approach to this description is to assume that 

the generalization can be accomplished by assigning to 

the dispersion coefficient the properties of a tensor. 

Justification of this method is based on experimental 

evidence which indicates that the magniture of the 

dispersion coefficient for a homogeneous medium is 

different in the directions transverse and parallel to 

the line of flow. Because the generalization depends on 

the ability to identify separately all given causes of 

spreading of the contaminant, such as the effect of 

the injection source, the tensor components are not yet 

quantified.

Additional mechanisms, such as mass transfer from 

the liquid to the solid phase due to adsorption, are of 

fundamental importance whenever it is required that 

the fate of a given substance in a ground-water body 

be determined. The adsorption process may be described by the historical evaluation of the mechanism 

involved, and its effect on the transport of the given 

substance through the field of flow. In many instances 

of transport through a porous medium, the large solid 

surface area available makes the adsorption phenomenon the most important feature. Description of this 

phenomenon involves a chemical process; hence, the 

thermodynamics of the mass transfer must be determined for detailed analysis. At present, for lack of 

better data, approximation of the adsorption rate is 

used in the mathematical analysis. There is, as yet, no 

experimental evidence to indicate the validity of these 

expressions used in the mathematical model.

In all attempts to describe the flow conditions in a 

porous matrix, it has been found that the generalization 

of the expression for flow is formidable piimarily 

because of the inability to represent statistically the 

physical properties of the given medium in relation to 

the fluid flow (Skibitzke and others, 1961). Thus, the 

limitations placed on the analytical theories of dispersion are apparent. The generalized description of the 

transport of mass within any heterogeneous porous 

medium is unattainable by any presently known 

methods. However, the probelms encountered in a 

real situation do not require complete and rigorous 

formulation; usually, an engineering estimate of the 

maximum possible effect under the observed conditions 

is acceptable.

The investigations of mass-transport phenomena have 

proceeded along both analytical and experimental lines. 

The experimental methods have been used to verify the 

analytical expressions obtained by solving the differential equation of mass transport. Qualitative data 

defining the phenomenological process have also been 

obtained in the laboratory to demonstrate the nature 

of the flow. The analytical developments have proceeded 

along established lines to obtain expressions for predicting the possible concentration of a given dissolved 

substance under various specified conditions. As with 

all analyses, their applicability is limited by the ability 

to describe the parameters that control the behavior 

mechanisms; hence, only the most simple linear flow 

systems are described.

Within a limited range, the theories of dispersion 

have been verified by the results of experimental 

analysis. However, these expressions are only useful 

whenever ground-water investigations can be brought 

into the laboratory; for example, when the use of tracers 

is needed in establishing certain details of the fluid flow. 

Extension of such investigations into more complex 

porous media depends on successfully directing the 

research in dispersion theory so that the requisite 

studies of the controlling parameters are possible.

REFERENCES CITED

Adamson, A. W., 1960, Physical chemistry of surfaces: New

York, Interscience Publishers, 629 p. 

Aris, R., and Amundson, N. R., 1957, Some remarks on logitudinal mixing or diffusion in fixed beds: Am. Inst. Chem.

Eng. Jour., v. 3, no. 2. p. 280-281. 

Bachmat, Y., 1967, On the similitude of dispersion phenomena

in homogeneous and isotropic porous medium: Water

Resources Research, v. 3, no. 4, p. 1079-1083. 

Bear, Jacob, 1961, On the tensor form of dispersion in porous

media: Jour. Geophys. Research, v. 66, no. 4, p. 1185-1198. 

Bear, Jacob, and Todd, D. K., 1960, The transition zone between

fresh and salt waters in coastal aquifers: California Univ.,

Water Resource Center Contribution 29, 156 p. 

Beran, M. J., 1957, Dispersion of soluble matter in flow through

granular media: Jour. Chem. Physics, v. 27, no. 1, p.

270-274. 

Biggar, J. W., and Nielsen, D. R., 1960, Diffusion effects in

miscible displacement occurring in saturated and unsaturated porous materials: Jour. Geophys. Research, v. 65,

no. 9, p. 2887-2896. 

Bird, R. B., Stewart, W. E., and Lighfoot, E. N., 1960, Transport

phenomena: New York, John Wiley & Sons, 780 p.
134 FLUID MOVEMENT IN EARTH MATERIALS

Carman, P. C., 1956, Flow of gases through porous media: New

York, Academic Press Inc. 

Carrier, G. F., 1958, The mixing of ground water and sea water

in permeable subsoils: Jour. Fluid Mech., v. 4, no. 5, p.

479-488. 

Carslaw, H. S., and Jaeger, J. C., 1959, Conduction of heat in

solids: London, Oxford Univ. Press, 510 p. 

Cooper, H. H., Jr., 1959, A hypothesis concerning the dynamic

balance of fresh water and salt water in a coastal aquifer:

Jour. Geophys. Research, v. 64, no. 4 p. 461-467. 

Crank, J., 1956, The mathematics of diffusion: London, Oxford

Univ. Press, 347 p. 

Dankwerts, P. V., 1953, Continuous flow system distribution of

residence times: Chem. Eng. Sci., v. 2, no. 1, p. 1-13. 

Day, P. R., 195(5, Dispersion of moving salt-water boundary

advancing through saturated sand: Am. Geophys. Union

Trans., v. 37, no. 5, p. 595-601. 

ErdeUyi, A., Magnus, W., Oberhettinger, F., and Tricomi, F. G.,

1954, Tables of integral transforms, v. 1: New York,

McGraw-Hill Book Co., 391 p. 

Glasstone, Samuel, 1961, Textbook of physical chemistry:

Princeton, N. J., D. Van Nostrand Co., 1320 p. 

Harleman, D. R. F., and Rumer, R. R., Jr., 1962, The dynamics

of salt-water intrusion in porous media: Massachussetts

Inst. Technology Hydrodynamics Lab. Rept. 55, 125 p. 

Henry, H. R., 1960, Salt intrusion into coastal aquifers: Internat.

Assoc. Sci. Hydrology, Comm. Ground Waters Pub. 52, p.

478-487. 

Hougen, 0. A., and Marshall, W. R., Jr., 1947, Adsorption from

a fluid stream flowing through a stationary granular bed:

Chem. Eng. Process, v. 43, no. 4, p. 197. 

Jeffreys, Harold, 1961, Cartesian tensors: London, Cambridge

Univ. Press, 92 p.

Jong, G. de Josselin de, 1958, Longitudinal and transverse diffusion in granular deposits: Am. Geophys. Union Trans., v.

39, no. 1, p. 67-74.

 1959, Vortex theory of multiple phase flow through porous 

media: California Univ. Water Resource Center Contribution 23, 80 p. /

Jost, W., 1960, Diffusion in/solids, liquids, gases: New York, 

Academic Press Inc., 55o p.

Kitagawa, K., 1934, Sur le dispersement et l'e"cart moyen de 

l'6conlement des eaux sonterraines, I. Experiences avec un 

modele le laboratoire: Kyoto Imperial Univ., Coll. Sci. Mem. 

Ser. A, v. 17, p. 37-42.

Kohout, F. A., 1960, Cyclic flow of salt water in the Biscayne 

aquifer of southeastern Florida: Jour. Geophys. Research, v. 

65, no. 7, p. 2133-2141.

 1961, Fluctuation of groundwater levels caused by dispersion of salts: Jour. Geophys. Research, v. 66, no. 8, p. 

2429-2434.

Kramers, H., and Alberda, G., 1953, Frequency response analysis 

of continuous flow system: Chem. Eng. Sci., v. 2, no. 4, p. 

173-181.

Lau, Leung-Kai, Kaufman, W. J., and Todd, D. K., 1959, 

Dispersion of a water tracer in radial laminar flow through 

homogeneous porous media: Calif. Univ., Inst. Eng. 

Research Ser. 93, Issue No. 5, 78 p.

Mime-Thomson, L. M., 1961, Theoretical hydrodynamics: New 

York, MacMillan Co., 660 p.

Muskat, Morris, 1946, The flow of homogeneous fluids through 

porous media: Ann Arbor, Michigan, J. W. Edwards, Inc., 

749 p.

Ogata, Akio, 1958, Dispersion in porous media: Evanston, 111., 

Northwestern Univ., Ph. D. dissert.

 1961, Transverse diffusion in saturated isotropic granular 

media: U.S. Geol. Survey Prof. Paper 411-B, 8 p.

 1964a, The spread of a dye stream in an isotropic granular 

medium: U.S. Geol. Survey Prof. Paper 411-G, 11 p.

1964b, Mathematics of dispersion with linear adsorption

isotherm: U.S. Geol. Survey Prof. Paper 411-H, 9 p.

Ogata, Akio, and Banks, R. B., 1961, A solution of the differential equation of longitudinal dispersion in porous media: 

U.S. Geol. Survey Prof. Paper 411-A, 7 p.

Orlob, G. T., and Radhakrishna, G. N., 1958, The effects of 

entrapped gases on the hydraulic characteristics of porous 

media: Am. Geophys. Union Trans., v. 39, no. 4, p. 648-659.

Perkins, T. K., and Johnston, O. C., 1963, A review of diffusion 

and dispersion in porous media: Soc. Petroleum Engineers 

Jour., v. 3, no. 1, p. 70-84.

Rifai, M. N. E., Kaufman, W. J., and Todd, D. K., 1956, 

Dispersion phenomena in laminar flow through porous 

media: Calif. Univ., Inst. Eng. Research Ser. 93, Issue 

No. 2, 157 p.

Saflfman, P. G., 1959, A theory of dispersion in a porous medium: 

Jour. Fluid Mech., v. 6, p. 321-349.

 1960, Dispersion due to molecular diffusion and macroscopic mixing in flow through a network of capillaries: Jour. 

Fluid Mech., v. 7, p. 194-208.

Scheidegger, A. E., 1954, Statistical hydrodynamics in porous 

media: Appl. Physics Jour., v. 25, p. 994-1001.

 1960, The physics of flow through porous media: Canada, 

Univ. Toronto Press, 313 p.

1961, General theory of dispersion in porous media: Jour.

Geophys. Research, v. 66, no. 10, p. 3273-3278. 

Shamir, U. Y., and Harleman, D. R. F., 1967, Numerical

solutions for dispersion in porous media: Water Resources

Research, v. 3, no. 2, p. 557-581.

Skibitzke, H. E., 1964, Extending Darcy's concept of groundwater motion: U.S. Geol. Survey Prof. Paper 411-F, 6 p. 

Skibitzke, H. E., Chapman, H. T., Robinson, G. M., and

McCullough, R. A., 1961, Radio-tracer techniques for the

study of flow in saturated porous material: Internat. Jour.

Appl. Radiation and Isotopes, v. 10, no. 1, p. 38-46. 

Skibitzke, H. E., and Robinson, G. M., 1963, Dispersion of

ground water flowing through heterogeneous materials: U.S.

Geol. Survey Prof. Paper 386-B, 3 p. 

Slichter, C. S., 1905, Field measurements of the rate of movement

of underground water: U.S. Geol. Survey Water-Supply

Paper 140, p. 9-85.

Sneddon, I. N., 1951, Fourier Transforms: New York, McGrawHill Book Co., 542 p. 

Taylor, Geoffrey, 1953, Dispersion of soluble matter in solvent

flowing slowly through a tube: Royal Soc. London Proc.,

ser. A, v. 219, p. 187-203. 

Vermeulen, Theodore, and Hiester, N. K., 1952, Ion exchange

chromatography of trace elements: Ind. Eng. Chemistry,

v. 44, p. 636-651. 

Wentworth, C. K., 1948, Growth of the Ghyben-Herzberg

transition zone under a rinsing hyopthesis: Am. Geophys.

Union Trans., v. 29, p. 97-98.

U. S. GOVERNMENT PRINTING OFFICE : 1970 O - 370-837