Transverse Diffusion in 

Saturated Isotropic 

Granular Media

GEOLOGICAL SURVEY PROFESSIONAL PAPER 411-B
Transverse Diffusion in 

Saturated Isotropic 

Granular Media

By AKIO OGATA

FLUID MOVEMENT IN EARTH MATERIALS

GEOLOGICAL SURVEY PROFESSIONAL PAPER 411-B

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

STEWART L. UDALL, Secretary

GEOLOGICAL SURVEY 

Thomas B. Nolan, Director

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

Washington 25, D.C. - Price 15 cents (paper cover)
CONTENTS

	Page 

Abstract ----.-----.------__-______.___-___-_-___--_-__-_- _-__ B-1

Introduction. ______________________________________________________________ 1

Acknowledgments. _ _-_________---___-________-_-__-____-__________------_--- 1

Two-dimensional dispersion._____--_________--_---_-__---____-____----__--_-- 1

Solution of the differential equation_________--_-_---_--___-----_------_------- 2

Solution of special cases of r equals 0 and r equals a______________-___-__-___--_- 4

Some properties of the hypergeometric series..--.___________-_____----_-_-_-- 5

General solution for computation of concentration distribution__________________ 5

Steady-state solution.-..__________________________ ___________-__-_--__ 6

Conclusion-.-....-.--.--_-__--------__--____-__________-----_____--_-_-_--- 6

References...----..------_-------__--___-_-___-______--__-_-----_---------- 8

ILLUSTRATIONS

FIGUBE 1. Mathematical model_____________________________________________ B-2

2. Plot of solution for various values of r/a_ _____-______--____-_-_----_- 7

3. Plot of solution for various values of a2/4Di________________________ 7

m

FLUID MOVEMENT IN EARTH MATERIALS

TRANSVERSE DIFFUSION IN SATURATED ISOTROPIC GRANULAR MEDIA

By AKIO OGATA

ABSTRACT

An analytical method of determining the ionic diffusion 

transverse to the direction of flow in granular media is presented. 

The basic assumption that the frontal zone created by longitudinal dispersion is fully developed and stationary before the 

transverse diffusion can cause spreading is made to simplify the 

mathematical analysis. Solutions were obtained in terms of 

both a series of Bessel's functions and a hypergeometric series. 

Numerical computations were made and results are presented 

in graphical form.

INTRODUCTION

The necessity for disposal of radioactive and other 

wastes and for more quantitative knowledge of microscopic flow in porous media has focused interest on 

dispersion phenomena occurring in saturated flow 

through porous media. Various laboratory techniques 

have been developed for determining the dispersion 

coefficient appearing in the analytical solutions. 

Emphasis, however, has been placed on studies of 

longitudinal dispersion for example, those of Beran 

(1957) and Day (1956) and by the writer (Dispersion in 

porous media, doctoral dissertation, Northwestern 

Univ., 1958). The data available on transverse 

diffusion are still extremely meager, owing to the 

difficulty in modeling transverse diffusion and 

measuring the concentration distribution once a model 

is constructed.

The Geological Survey in its ground-water research 

office at Phoenix, Ariz., has developed methods for 

measuring transverse diffusion, by means of radioactive 

isotopes and commercially available counters (Skibitzke 

and others, 1960). However, the studies to date have 

been directed toward feasibility of the methods, and 

quantitative data are not available. This paper will 

consider the mathematical aspects of an approximate 

description of the phenomena by dealing with a simplified 

mathematical model.

ACKNOWLEDGMENTS

The problem treated in this report was brought to the 

attention of the writer by H. E. Skibitzke, mathematician in charge of the Phoenix ground-water research 

office. The writer is greatly indebted to Mr. Skibitzke, 

whose initial development and ideas greatly facilitated 

the preparation of this report.

The author would like to express appreciation also to 

Prof. Richard Skalak, of Columbia University, for his 

thorough discussion and review of the paper. Professor 

Skalak points out that Goldstein, in his paper "Some 

Two-Dimensional Diffusion Problems with Circular 

Symmetry" (1932), obtained equation 27. However, 

his paper was not available to the writer at the time of 

writing. Professor Skalak points out also that there 

is a shorter method of derivation of equation 9 by use of 

the Hankel transform. The author utilized the Laplace 

transformation because this method has a wider range 

of applicability and is a more familiar method of 

approach.

TWO-DIMENSIONAL DISPERSION

The differential equation that describes the twodimensional dispersion of a contaminant hi saturated 

granular media, where radial symmetry exists, is

(D

where C= concentration of contaminants hi lb/ft3 

r=radial distance, in feet, measured from

center of contaminant filament 

£=longitudinal distance, in feet, measured

in the direction of fluid flow 

£=time, in seconds

Z>r=radial diffusion coefficient, in ft2/sec 

Z>z=longitudinal dispersion coefficient, in ft2/

sec 

u=average fluid velocity, hi ft/sec

Equation 1 is readily obtained from the law of conservation of mass. Owing to difficulty in solving the 

general case, it is necessary to make some simplifying 

assumptions that are physically realistic.

B-l

584498 61-
B-2 FLUID MOVEMENT IN EARTH MATERIALS

The physical nature of the dispersion process is shown 

in figure 1. The system may be divided into two regions, the frontal zone and the cylindrical filament of 

contaminant upstream from the frontal zone. In the

Contaminant 

supply

Fluid 

reservoir

Granular media.

Fluid velocity u - constant

Circular source

Filament of 

contaminant

h-H

Frontal zone 

of dispersion

FIGURE l. Physical setting of mathematical model.

frontal zone the dispersion, or spreading, of the contaminant is due to longitudinal plus transverse dispersion. However, in the region upstream from this 

frontal zone the spreading is due principally to transverse diffusion and convection because of the virtual 

nonexistence of a concentration gradient in the direction 

of flow. Thus, if the frontal zone is not considered,

&C £*zr~2 =0 and equation (1) reduces to

A|/r ac\ ac ac

ror\or/ ox ot

In discussing the frontal zone, consider a moving 

coordinate system that is, let %=x ut and r=t. 

Substituting this transformation into equation (1) gives

_

Physically, the above equation indicates that the 

observer is moving at a velocity, u.

Experimental and analytical development of the 

dispersion in a one-dimensional system indicates that 

the length of the frontal zone remains constant, once 

it is fully developed. Von Kosenberg (1956) indicates 

that the zone in which the concentration varies by a 

certain percentage depends on the velocity and the 

dispersion coefficient. In the system considered here, 

in which no mass exchange occurs between the liquid 

and the solid phase, this front will progress through 

the medium with the velocity of the fluid.

Because of the existence of another interface in the 

radial direction, the two-dimensional system is not the

same as the one-dimensional system. But, if it is 

assumed that the establishment of this front takes 

place at a relatively rapid rate in comparison with 

transverse diffusion, then the tune rate of change of 

concentration depends only on the radial dispersion. 

Thus, the differential equation describing the process is

Substituting %=x ut, t=r in equation 2 and resubstituting t for r gives

Ur o f OO \_OC/ (3)

Equation 3 is the fundamental differential equation 

describing the dispersion process throughout the region 

for long times. This is the differential equation found 

in various radial-diffusion or heat-flow problems, whose 

solutions are known for various boundary conditions. 

The boundary conditions for the case to be considered 

are

(2) £=0, (7=0, r

(3) r=0,= (4)

It is evident that an initial-value problem needs to be 

solved. Because the longitudinal-dispersion coefficient, 

DX) no longer appears, the subscript r will be dropped 

from the radial-diffusion coefficient.

SOLUTION OF THE DIFFERENTIAL EQUATION

Solution of equation 3 may be obtained by employing 

the Laplace transform; that is, assume that there exists 

a function

, t)dt

If equation 3 is multiplied by e~pt and integrated with 

respect to t within the limits 0 to <» , equation 3 may 

be written

(5a)

_ _ n 

It should be noted that by setting C \ - in 5a

the two equations become equivalent.
TRANSVERSE DIFFUSION IN SATURATED ISOTROPIC GRANULAR MEDIA

Hence, equation 5 may also be written

B-3

p. D'

Above is Bessel's equation, whose solution is as 

follows (Hildebrand, 1954, p. 167)

where /<,(«) =modified Bessel function of the

first kind of order zero 

K0 (a)=modified Bessel function of the

second Mnd of order zero 

A,B=arbitrary constants

Since K0 remains finite as r >oo ; this solution is valid; 

that is, .4=0 for r^>a. In the region r<X /o remains /de\ finite as r >0, and also since I ~>T J must be zero, 70 

is the solution, or 5=0. Hence,

r<a

Atr=a Ci=Ci; thus

In addition, the mass rate of transfer across r=a is 

given by

giving

 BKi(aq)=AI1 (qa)

From the above two equations the coefficients A and B 

may be determined; that is,

P

p K0(ga)Ii(qa)+Ki(qa)I0 (q[)

Substituting values for A and B into equation 6 gives 

the Laplace solution of equation 5,

P a) +K0 (qa)Ii(ga)

(7)

"I

To obtain the solution of equation 3 for C as a function 

of r and t it is necessary to apply the Bromwich inversion 

theorem,

1 C+i<* _

e"C(r,p)dp

fC+i =^-.

61*% Jc+i

The inversion of equation 7, however, has been evaluated by Carslaw and Jaeger (1948, p. 285); the 

result for this specific case is

-Deft Ji(aa)[J0 (ar)4>(a)]da 

' «[02 («)] 

J0 (ar)Ji (aa)da C2 (r,t)=-

TT

where

0(a) = Jl (aa)Y0 (aa) J0 (dot) YI (aa)

Further (Watson, 1948, page 77),

(8)

_^2

TTZ

Thus, because v=Q,

2 

iraa

Substituting the value for 0(a), the result is

^=0^=0=0^

Jo

(9)

where

J\(a} and J0 («)=Bessel function of the first kind of

order one and zero, respectively 

a=radius of contaminant filament

By direct substitution, the above can be shown to 

satisfy the original differential equation. It is now 

desired that equation 9 be shown to satisfy the 

boundary conditions. For t 0 the solution reduces to 

the following (Watson, ] 948, p. *06):

 =af e

CQ Jo

Ji(aa)Jo(ar)da

 a*

"0 for

 -

2a

 I/a

(10)

The condition at r=a is due to an averaging process 

such that the concentration curve becomes a continuous 

function.

In any numerical integration process difficulties are 

inherent; thus, if possible, equation 9 should be written 

hi terms of tabulated functions to facilitate computation. However, before going into the general case it
B-4 FLUID MOVEMENT IN EARTH MATERIALS

may be advantageous to consider special cases for r=0 

and r=a. For these two values of r, equation 9 may be 

readily integrated, thereby simplifying the analysis considerably. The diffusion coefficient may be readily 

obtained from experimental data by using these two 

expressions.

SOLUTION OF SPECIAL CASES OF r=0 AND r=a

Consider first the concentration variation along r=0. 

Because J0 (x) *l as x-*Q the equation above becomes 

simply

c r

7T=<H

^o Jo

,-Dtc? (11)

Integration of the above is given by Watson (J 948, p. 

394), the result being

But

\ VT )=t/ Sinh a;

1TX

Hence, letting 7=37^* the above may be written

oUt

C 

Further,

^-=2 exp ( 7) Sinh (7) 

Sinh (»)=i (c* O

Thus, the concentration distribution at r=Q is given by 

gr=l-exp (-2T) (13)

Note that this expression may be obtained directly from 

the more general solution, equation 24.

However, if the radius of filament a is large enough 

that the concentration of the centerline does not vary 

sufficiently in the experimental model, it would be 

advantageous to choose some other value of r. In this 

instance the best value would be r=a, provided time 

is large. Substituting for r in equation 9 gives

C f 77=a I e-

^o Jo

Dta2 Ji(aa)J0 (aa)da

The integration of the above is given by Bateman 's 

Manuscript Project (1953b, p. 50), expressed in terms 

of a hypergeometric series; that is,

CL where z=jji' Some properties of the hypergeometric

series will be presented in a latter part of the paper. 

Equation 14 may be written as a series

or

C 1

. . . An zn)

where values of An are positive coefficients. Also 

note that

Since (l) n n\; (2) =(?&-{-1)! equation 14 further 

reduces to

1 /3N

Co

2V2

nlnl

n=*l

The above expression is the series representation of 

the confluent hypergeometric function; thus it may be 

written

where

However, since a relationship exists between the 

confluent hypergeometric function and the Bessel 

function, it would be advantageous to write equation 

15 in terms of Bessel functions. Bateman (1953a, p. 

265) gives the relationship

; 1+2v;

where r(i>+l) is the 

70(x)=/o( »); thus

gamma function. Further,

(16)

The concentration distribution along r=a is readily 

computed by means of equation 16.

The function /0 (a?) has been tabulated extensively; 

thus there is no difficulty in computing equation 16. 

As stated previously, equations 13 and 16 provide a 

means of computing the radial-diffusion coefficient from 

experimental data.
TRANSVERSE DIFFUSION IN SATURATED ISOTROPIC GRANULAR MEDIA 

SOME PROPERTIES OF THE HYPERGEOMETRIC SERIES

The hypergeometric function occurs as solution of a 

linear second-order differential equation, which has at 

most three singularities at 0, c°, and 1. In this paper 

the Pochhammer notation is used. By definition, the 

generalized hypergeometric series is as follows (Watson, 

1948, p. 100; Bateman Manuscript Project, 1953a, 

p. 56):

; pi, p2 ,

pe ',

where

r? n

(a+n-l);(a) 0 =l

It should be noted that if p=2 and q=l the above 

reduces to

(18) _f -i/ N «=o n\(p) n

Note further that if al = I and <x2 =P in equation 18 this 

series becomes an elementary geometric series, or

n=0

Owing to extreme difficulties in dealing with the 

generalized hypergeometric function, only the convergence of equation 18 will be stated. For a more 

general discussion of the hypergeometric series the 

reader is referred to the Bateman Manuscript Project 

(1953a, p. 56).

A hypergeometric series of the type given by equation 18 is generally convergent for values of |2|<OSeries converge for x 1 only if p a1 «2>0 and for 

x= l only if p ai a2 +l>0 (Hildebrand, 1954, 

p. 179).

GENERAL SOLUTION FOR COMPUTATION OF 

CONCENTRATION DISTRIBUTION

As stated previously, although equation 9 is the 

required solution of the problem, the numerical evaluation of the integral is extremely difficult. Thus, it 

would be advantageous to rewrite it in such a way as to 

facilitate computation. Equation 9 may be written in 

a series of Bessel functions or hypergeometric functions. 

Because tables of Bessel functions were readily available, this series was used. However, both of the series 

will be presented.

The product of Bessel functions generally may be 

written

= l An Jn (aO)Jn(ar) 

»=o

by use of the recurrence formula

Explicitly,

B-5

(19)

(f)©'

«=2

(f)wn l -Jn (aa)Jn (<xr) (20)

Substituting equation 19 into equation 9 and using the 

known integral

= exp

(21)

Equation 9 may be readily integrated. The result of 

the integration is

Q~~ C

where £=

_?!±^i /

±Dt J \

, a and = '

1=2 m J (22)

By use of the recurrence relationship for the modified 

Bessel function,

(23) _2v T 

The second term in equation 22 may be written

co («) w ~2 co ^m oo

TO=2 ffl' -I 7n=0 "*"i -^ TO=I

This on substituting into the original equation gives

b r?w/w (£) (24)

The above result was obtained by Goldstein (1953) in 

the study of exchange processes in fixed columns. 

Equation 24 may be written in another form, however :

£=exp f

OQ L m=l

(25)

Equation 25 also may be written in its alternate form 

(Goldstein, 1953)

m=0 ^

(26)
B-6 FLUID MOVEMENT IN EARTH MATERIALS

Note that for ?j=->l, <1; thus equation 26 converges more rapidly than 25 in the region ->1. The T

converse is true for equation 25; that is, it would be 

advantageous to use equation 25 for computation in

the region -<1. It should be noted further that 

r

equations 13 and 16 may be obtained directly from 

equations 25 and 26.

As stated previously, equation 9 may be written in 

terms of the hypergeometric series. There exist the 

relationships (Watson, 1948, p. 148)

J0(ar)Ji(aa)

/fta\ ^-^, 

= -7T S

( m 1 m;!;-

for i7 2==-r >1, and ri

da

m=0 mlml

Substituting the above relationships into equation 9 

and using the known integral

/.'

«/ o

2pl/2UH

the expressions obtained are

(27)

and

(28)

where V X\=Tfn' 

4JJt

Hypergeometric functions appearing as coefficients 

of the power series are polynomials of m.

Although both equations 27 and 28 and equations 25 

and 26 may be readily computed, equations 25 and 26 

were used because of their rapid convergence. A plot

f C a f . , f a . . . oi 7=7- versus -7-5=77 for vanous values of rt=- is given in GO 4Z?£ r &

C a figure 2, and a plot of -~- versus - for various values of

 J=T: is given in figure 3.

STEADY-STATE SOLUTION

Owing to the nature of the experiments conducted 

(Skibitzke and others, 1960), it would be advantageous 

to assume that a steady state is reached. Accordingly,

^v VV

consider -^-7=0 in equation 2, which gives

r -^--^- -Tr or \r or/ ax

The boundary conditions become

<7(r,0)=0

/00, (29)

The conditions given are the same as those given for the

np

unsteady-state problem. Further, by letting T=-,

u

equations 29 and 3 are identical. Thus, the solution fo r 

steady-state condition may b 

in equation 25, 26, 27, or 28.

steady-state condition may be obtained by letting t -

u

CONCLUSION

In attacking the mathematical problem of transverse 

diffusion, an important assumption was made: that the 

front due to the longitudinal dispersion is established 

rapidly as compared to the transverse diffusion. If the 

diffusion transverse to the direction of flow includes 

mechanical dispersion, this assumption would be erroneous. However, since a homogeneous medium is 

assumed, the assumption that diffusion is wholly due 

to molecular agitation seems to be valid. Although 

there have been no quantitative studies, qualitative 

study using dye tracers in various types of models 

seems to bear this out.

Up to this time, investigation has been directed 

primarily toward evaluating experimental and analytical methods. The approximate method presented leads 

to what is believed to be a realistic solution that can 

be computed. When laboratory experiments are completed the resulting data, in conjunction with the 

solution presented, will furnish a means of determining 

the magnitude of the coefficient of transverse diffusion.
TRANSVERSE DIFFUSION IN SATURATED ISOTROPIC GRANULAR MEDIA B-7

FIGTJEE 2. Plot of solution for various values of r/o.

0.9 -

0 O.I 0.2 0.3 0.4 0.5 0.6 0,7 0.8 0.9 3 |.4 1.5 1-6 1.7 1.8 1.9 2.0

FIGTJEE 3. Plot of solution for various values of (H/4ZW.
B-8 FLUID MOVEMENT IN EARTH MATERIALS

REFERENCES

Bateman Manuscript Project, 1953a, Higher transcendental

functions, v. 1: New York, McGraw-Hill, 302 p. 

 1953b, Higher transcendental functions, v. 2: New

York, McGraw-Hill, 396 p. 

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

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

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

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

Day, P. R., 1956, Dispersion of moving salt-water boundary

advancing through saturated sand: Am. Geophys. Union

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

Goldstein, S., 1932, Some two-dimensional diffusion problems

with circular symmetry: London Math. Soc. Proc., Ser. 2,

34, p. 51-88.

Goldstein, S., 1953, On the mathematics of exchange processes in 

fixed columns: Royal Soc. (London) Proc., v. 219, p. 151-185.

Hildebrand, F. B., 1954, Advanced calculus for engineers: New 

York, Prentice-Hall, 594 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. 

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

Von Rosenberg, D. U., 1956, Mechanics of steady-state singlephase fluid displacements from porous media: Am. Inst. 

Chem. Eng. Jour., v. 2, p. 55.

Watson, G. N., 1948, A treatise on the theory of Bessel functions; 

Cambridge Univ. Press, 804 p.

o