A Solution of the 

Differential Equation of

Longitudinal Dispersion 

in Porous Media

GEOLOGICAL SURVEY PROFESSIONAL PAPER 411-A
A Solution of the 

Differential Equation of 

Longitudinal Dispersion 

in Porous Media

By AKIO OGATA and R. B. BANKS

FLUID MOVEMENT IN EARTH MATERIALS

GEOLOGICAL SURVEY PROFESSIONAL PAPER 411-A

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 Documents5 U.S. Government Printing Office 

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

ILLUSTRATIONS

Abstract-___-___________-______-______________-------__--------__------_- A-l

Introduction_____________________________________________________________ 1

Basic equation and solution__________________________.___________________ 1

Evaluation of the integral solution________________________________________ 3

Consideration of error introduced in neglecting the second term of equation 8____ 4

Conclusion_____________________________________________________________ 7

References ______________________________________ ____ 7

Page 

FIGURE 1. Plot of equation 13____________________________________________ A-5

2. Comparison of theoretical and experimental results__________--_-__ 5

3. Plot of equation 18 -_-______________-____--_-_----__--_-----_- 6

m

FLUID MOVEMENT IN EARTH MATERIALS

A SOLUTION OF THE DIFFERENTIAL EQUATION OF LONGITUDINAL DISPERSION IN

POROUS MEDIA

By AKIO OGATA and K. B. BANKS

ABSTRACT

Published papers indicate that most investigators use the 

coordinate transformation (x ut) in order to solve the equation 

tor dispersion of a moving fluid in porous media. Further, the 

boundary conditions O=0 at x=«> and 0=00 at x= «> for 

Z>0 are used, which results in a symmetrical concentration distribution. This paper presents a solution of the differential 

equation that avoids this transformation, thus giving rise to an 

asymmetrical concentration distribution. It is then shown that 

this solution approaches that given by symmetrical boundary 

conditions, provided the dispersion coefficient D is small and 

the region near the source is not considered.

INTRODUCTION

In recent years considerable interest and attention 

have been directed to dispersion phenomena in flow 

through porous media. Scheidegger (1954), deJong 

(1958), and Day (1956) have presented statistical 

means to establish the concentration distribution and 

the dispersion coefficient.

A more direct method is presented here for solving 

the differential equation governing the process of dispersion. It is assumed that the porous medium is 

homogeneous and isotropic and that no mass transfer 

occurs between the solid and liquid phases. It is 

assumed also that the solute transport, across any fixed 

plane, due to microscopic velocity variations in the 

flow tubes, may be quantitatively expressed as the 

product of a dispersion coefficient and the concentration 

gradient. The flow in the medium is assumed to be 

unidirectional and the average velocity is taken to be 

constant throughout the length of the flow field.

BASIC EQUATION AND SOLUTION

Because mass is conserved, tl^e governing differential 

equation is determined to be

d<7

(1)

v

where

D=dispersion coefficient 

C= concentration of solute in the fluid

u=average velocity of fluid or superficial velocity/ 

porosity of medium

x=coordinate parallel to flow 

y,z coordinates normal to flow 

2=time.

In the event that mass transfer takes place between 

the liquid and solid phases, the differential equation 

becomes

_ 5(7 d(7 &F

where F is the concentration of the solute in the solid 

phase.

The specific problem considered is that of a semiinfinite medium having a plane source at z=0. Hence 

equation 1 becomes

Initially, saturated flow of fluid of concentration, 

(7=0, takes place in the medium. At t Q, the concentration of the plane source is instantaneously changed 

to (7=(70 . Thus, the appropriate boundary conditions 

are

<7(co,£)=0; *>

The problem then is to characterize the concentration 

as a function of x and t.

To reduce equation 1 to a more familiar form, let

(4)

A-l

586211 61 2
A-2 FLUID MOVEMENT IN EARTH MATERIALS

Substituting equation 4 into equation 1 gives

The boundary conditions transform to

= Co exp (

It is thus required that equation 5 be solved for a timedependent influx of fluid at 2=0.

The solution of equation 5 may be obtained readily 

by use of Duhamel's theorem (Carslaw and Jaeger, 

1947, p. 19):

If C=F(x,y,z,t) is the solution of the diffusion equation for semi-infinite media in which the initial concentration is zero and its surface is maintained at concentration unity, then the solution of the problem in which 

the surface is maintained at temperature <j>(t) is

F(x,y,z,t \)d\.

This theorem is used principally for heat conduction 

problems, but the above has been specialized to fit 

this specific case of interest.

Consider now the problem in which initial concentration is zero and the boundary is maintained at 

concentration unity. The boundary conditions are

r(z,0)=0;z>0

This problem is readily solved by application of the 

Laplace transform which is defined as

/»

Jo

T(x,p)= e-ptT(x,t)dt. 

Jo

Hence, if equation 5 is multiplied by e~pt and integrated 

term by term it is reduced to an ordinary differential 

equation

(6) dx*~D' __ 

The solution of equation 6 is

where

The boundary condition as x >«> requires that B=0 

and boundary condition at x=0 requires that A=l/p, 

thus the particular solution of the Laplace transformed 

equation is

T=i e -v. 

P

The inversion of the above function is given in any 

table of Laplace transforms (for example, Carslaw and 

Jaeger, p. 380). The result is

Utilizing DuhameFs theorem, the solution of the 

problem with initial concentration zero and the timedependent surface condition at a; =0 is

Since e~^ is a continuous function, it is possible to 

differentiate under the integral, which gives

AA r

 JZ&J- e~n G

Thus

Letting

\=x/2^D(t T)

the solution may be written

(8)

Since <j>(t) = C0 exp (uH/^D) the particular solution of 

the problem may be written,

(9)

i where e a2-JDi
DIFFERENTIAL EQUATION OF LONGITUDINAL DISPERSION IN POROUS MEDIA 

EVALUATION OF THE INTEGRAL SOLUTION

A-3

The integration of the first term of equation 9 gives

(Pierce, 1956, p. 68)

For convenience the second integral may be expressed 

in terms of error function (Horenstein, 1945), because 

this function is well tabulated. 

Noting that

The second integral of equation 9 may be written

exp [-H)>}-

Since the method of reducing integral to a tabulated 

function is the same for both integrals in the right side 

of equation 10, only the first term is considered. Let 

z=e/\ and adding and subtracting

the integral may be expressed

Further, let

in the first term of the above equation, then

exp

J"

Similar evaluation of the second integral of equation 

10 gives

Again substituting 

result is

-*""£"" -

$=-t z into the first term, the

Noting that

substitution into equation 10 gives

J

oo / t 

J« 

 a 

a a

Thus, equation 9 may be expressed

2(70

-e* f (11)

However, by definition,

erfc

Writing equation 11 in terms of the error functions 

2

erfc «++e~2i erfc «-

Thus, substituting into equation 4 the solution is

g-i [erfc (a-i)+e- erfc («+£} (12)
A-4

Resubstituting for e and a gives

 £=i/ erfc e (£=«*Yf<P e erfc TtG

FLUID MOVEMENT IN EARTH MATERIALS

be reduced to

which may be written in terms of dimensionless 

parameters,

where l-=ut/x and i)=D/ux.

Where boundaries are symmetrical the solution of 

the problem is given by the first term of equation 13. 

This symmetrical system was considered by Dankwerts 

(1953) and Day (1956), utilizing different analytical 

methods. The second term in equation 13 is thus due 

to the asymmetric boundary imposed in the more 

general problem. However, it should be noted also 

that if a point a great distance away from the source 

is considered, then it is possible to approximate the 

boundary condition by C( °° ,f) = O0 , which leads 

to a symmetrical solution.

A plot on logarithmic probability graph of the above 

solution is given in figure 1 for various values of the 

dimensionless group ri=D/ux. The figure shows that 

as 77 becomes small the concentration distribution 

becomes nearly symmetrical about the value £=1. 

However, for large values of t\ asymmetrical concentration distributions become noticeable. This indicates 

that for large values of D or small values of distance x 

the contribution of the second term in equation 13 

becomes significant as £ approaches unity.

Experimental results present further evidence (for 

example, Orlob, 1958; Ogata, Dispersion in Porous 

Media, doctoral dissertation, Northwestern Univ., 

1958) that the distribution is symmetrical for values 

of x chosen some distance from the source. An 

example of experimental break-through curves obtained 

for dispersion in a cylindrical vertical column is shown 

as figure 2. The theoretical curve was obtained by 

neglecting the second term of equation 13.

CONSIDERATION OF ERROR INTRODUCED IN 

NEGLECTING THE SECOND TERM OF EQUATION 8

Experimental data obtained give strong indication 

that in the region of flow that is of particular interest 

it is necessary to consider only the first term of equation 

13. Owing to complexity of the overall problem of 

determining the error, it would facilitate analysis to 

determine the value of £ at which the function e17 " erfc 

/l+i\ is a maximum. This then will enable the determination of the value of rj at which equation 13 may

(14)

without introducing errors in excess of experimental 

errors.

The necessary condition that the function /0?,£) is a 

stationary point is given by

(15)

To determine whether the function is either maximum 

or minimum at a given point the sufficient conditions 

are given by

>0

O2f O27

(b) Maxima, p-p <0, consequently ^ <0

d2/ d2/

(c) Minima, ^ >0, consequently -^ >0 (16)

(Irving and Mullineux, 1959, pp. 183-187)

Further, if 16 (a) is greater than zero, the stationary 

point is called a saddle point. 

Let

Differentiating the function

where e2=(l £) 2/4£7/. From the above expression it 

can be seen that £=1 and £=<» are the stationary 

points of the function.

The second differentials can readily be obtained by 

direct methods. The results are

'[I

and

also

_

2 4
DIFFERENTIAL EQUATION OF LONGITUDINAL DISPERSION IN POROUS MEDIA A-5

FIGURE 1. Plot of equation 13.

2.0 -2.0

FiotJBE 2. Comparison of theoretical and experimental results.
A-6 FLUID MOVEMENT IN EARTH MATERIALS

At point £=1, the following expressions are obtained:

d2/

{=1

-8/2

=0.

It can be shown that

drj

(17)

<0

for all values of 17 by numerical consideration of the 

expansion of the complementary error function.

Accordingly, the function el/1> erfc ( T I satisfies

conditions 16 (a) and 16(b), indicating that maxima of 

the function occurs at {=1. The point at infinity is 

not considered here since both terms of equation 13 

approach zero as £ approaches infinity, which is indicative of a minimum condition.

There are no general analytical means available by 

which it is possible to obtain a general expression to 

determine the error involved in neglecting the second 

term of equation 13. Accordingly, consideration will 

be given in obtaining a reasonable numerical value of ?? 

for which the second term in equation 13 may be 

neglected. Consider equations 13 and 14 at £=1. Note 

that equation 14 reduces the value of J2. 1 while 

equation 13 gives C"o 2

where X= 1

(18)

The function e? erfc X is tabulated in

Carslaw and Jaeger (1948) up to the value X=3.0. For 

large values of X or small values of i\ the function may 

be approximated by

erfc X=-= -

'2X3 "22Xfi

Hence, equation 18 may be readily computed.

A semilogarithmic plot of equation 18 is given in

FIOCBE 3. Plot of equation 18.
DIFFERENTIAL EQUATION OF LONGITUDINAL DISPERSION IN POROUS MEDIA A-7

figure 3. It indicates that for values of i7<0.002 a 

maximum error of less than 3 percent is introduced by 

neglecting the second term of equation 13.

In all experiments reported, the dispersion coefficient 

D ranged from approximately 10~2 to 10~* cm2/sec. 

Further, it has been established that D is proportional 

to the velocity, hence the relationship D Dmu may be 

written, where Dm is the proportionality constant which 

is believed to be dependent on the media. Accordingly, 

since ii=Dlux, i\ may be expressed as i\ Dmtx. This 

then indicates that for £>5.DOTX102 the second term of 

equation 13 becomes negligible. Orlob and Radhakrishna (1958) obtained values of Dm ranging from 0.09 

cm to 2.79 cm. Using these values, measurements must 

be obtained at values of x greater than 45 cm or 1395 

cm. However, if an error of 5 percent is permitted, the 

above values are reduced by a factor of 4, thus x must 

be greater than 10 cm or 350 cm.

CONCLUSION

Consideration of the governing differential equation 

for dispersion in flow through porous media gives rise 

to a solution that is not symmetrical about x=ut for 

large values of 17. Experimental evidence, however, 

reveals that D is small. This indicates that, unless the 

region close to the source is considered, the concentration distribution is approximately symmetrical. TheoC 1

retically, 7? -> only as r?-»0; however, only errors of GO ^

the order of magnitude of experimental errors are 

introduced in the ordinary experiments if a symmetrical 

solution is assumed instead of the actual asymmetrical

one.

REFERENCES

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

solids: Oxford Univ. Press, 386 p. 

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

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

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

advancing through saturated sand: Am. Geophys. Union

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

deJong, G. de J., 1958, Longitudinal and transverse diffusion in

granular deposits: Am. Geophys. Union Trans., v. 39,

p. 67-74. 

Horenstein, W., 1945, On certain integrals in the theory of heat

conduction: Appl. Math. Quart., v. 3, p. 183-187. 

Irving, J., and Mullineux, N., 1959, Mathematics in physics and

engineering: Academic Press, 883 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, p. 648-659. 

Pierce, B. O., and Foster, R. M., 1956, A short table of integrals:

Boston, Mass., Ginn and Co., 189 p. 

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

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

O