See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/13227170

Do Hydrodynamic Dispersion Coefficients Exist?

Article  in  Physical Review Letters · December 1996

DOI: 10.1103/PhysRevLett.77.4552 · Source: PubMed

CITATIONS READS
61 96

2 authors, including:

Daan Frenkel
University of Cambridge
640 PUBLICATIONS   42,017 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Modelling phoretic effects View project

DNA coated colloids View project

All content following this page was uploaded by Daan Frenkel on 30 December 2013.

The user has requested enhancement of the downloaded file.

VOLUME 77, NUMBER 22 PHYSICAL REVIEW LETTERS 25 NOVEMBER 1996

Do Hydrodynamic Dispersion Coefficients Exist?

C. P. Lowe1,2 and D. Frenkel1

1
FOM Institute for Atomic and Molecular Physics, Kruislaan 407, 1098 SJ Amsterdam, The Netherlands
2
Computational Physics, Faculty of Applied Physics, Delft University of Technology, Lorentzweg 1, 2628 CJ Delft, The Netherlands
(Received 22 August 1996)
We have calculated dispersion coefficients for tracer particles in a fluid flowing through a porous
medium consisting of randomly packed spheres. At high Péclet numbers, where the motion of the tracer
is determined largely by convection, we found evidence that the dispersion coefficient is diverging and
that the dispersion process is anomalous. A transient region of anomalous dispersion has been predicted
theoretically. However, our simulations suggest that, rather than being transient, this effect persists. We
argue that our findings are consistent with the available experimental data. [S0031-9007(96)01712-7]

PACS numbers: 47.55.Mh, 05.60. + w

Tracer dispersion is the simplest example of the flow of the velocity variations multiplied by a correlation time
an inhomogeneous two-component mixture. We assume ty . For packed spheres, ty is typically taken to be
that the two components, the “tracer” and the “solvent,” the time it takes to convect over a particle diameter.
have identical behavior but are distinguishable in the sense There is experimental evidence that, at least for a medium
that they are assigned different “colors.” The model may consisting of randomly packed spheres, this description is
be simple, but it is neither trivial nor irrelevant. For in- adequate for Péclet numbers ranging from 3 3 102 to 1 3
stance, the solvent can represent a pure fluid and the tracer 105 [4,6]. The dispersivity has a value of a ­ 1.8 6 0.4
a contaminant. If the solvent occupies the space within [6,7]. Because of its apparent success on the “laboratory”
some porous medium and is stationary, the tracer has a scale the phenomenological approach forms the basis for
diffusion coefficient somewhat less than its value in the most models of dispersion on much larger scales [8].
pure solvent (the latter we denote as D0 ). If, however, the Apart from the difficulty of predicting dispersivities,
solvent flows through the porous medium with some mean the physics of tracer dispersion at high Péclet numbers
flow velocity v̄, then, depending on the magnitude of v̄, is, it seems, well understood. Below, we argue that
the dispersion coefficient D (defined as the tracer diffusion this is not the case. Some evidence that there is a
coefficient in a frame of reference moving with the mean problem comes from experiment. According to the Taylor
flow velocity) can be very much greater than D0 . Under- view, the dispersion coefficient should reach an effectively
standing this phenomenon of “hydrodynamic dispersion” constant value in the time it takes a tracer to convect
[1] is obviously important for modeling pollutant transport over a correlation length of the porous medium (typically
in ground water. of the order of a few particle diameters). However,
Hydrodynamic dispersion arises when the variations in experimentally D is still found to increase on a much
the local fluid velocity spread the tracer over a larger longer time scale [9,10]. Theoretical approaches which
volume than one would expect by diffusion alone (this treat the solid-fluid interface more explicitly also suggest
is nicely illustrated in Ref. [2]). The Péclet number Pe that things might not be quite so simple [11–13]. Koch
is a measure of the relative importance of diffusion and and Brady [12] showed that for an isolated sphere (i.e.,
convection. It is defined as Pe ; U p l p yD0 , where U p in the dilute limit) correlations in the fluid velocity decay
and l p are, respectively, a characteristic velocity and a so slowly that the hydrodynamic dispersion coefficient
characteristic length scale. At high values of Pe , tracer diverges [12]. That is, the mean square displacement, Dstd,
transport over distances larger than l p is dominated by grows as Dstd , t lnstd, rather than linearly in time. Koch
convection, and dispersion is dominated by the variations and Brady, however, went on to assume that away from
in the local fluid velocity. the dilute limit the velocity field far from the surface of
With the notable exception of the Taylor-Aris result for a given sphere is given by the Brinkman equation. This
a fluid flowing through a tube [3], there are few analytic is a mean-field theory which predicts that correlations in
results for D. A phenomenological argument [4] is that, at the velocity decay rapidly at distances greater than the
high Péclet numbers, molecular diffusion in porous media square root of the permeability (the Brinkman length).
ceases to play a role and that dispersion is determined For dense beds of spheres this is a small fraction of a
solely by spatial variations in the fluid velocity. In this particle radius. This “screened” decay is rapid enough
limit DyD0 ­ aPe . The constant, a, characteristic of a for the dispersion coefficient to be finite. However, there
given medium, is known as the dispersivity. This result should exist a boundary layer in which the velocity field
is derived using the Taylor hypothesis [5] which states which gives rise to the divergence of D in the dilute
that the dispersion coefficient equals the covariance of limit persists. For the dispersion coefficient to reach its

4552 0031-9007y96y77(22)y4552(4)$10.00 © 1996 The American Physical Society

VOLUME 77, NUMBER 22 PHYSICAL REVIEW LETTERS 25 NOVEMBER 1996

asymptotic value there must be time for the tracer to diffuse volume fraction of 0.45, were generated by a Monte Carlo
into (or out of) this layer, the width of which scales as simulation of a hard sphere fluid. This volume fraction
P 21y3
e d (d is the sphere diameter). On shorter time scales is somewhat lower than that of the beds used experimen-
the anomalous dispersion should persist [13]. This would tally. Salles et al. [19] reported only a weak dependence
have three important consequences. First, the dependence of the dispersivity on the volume fraction (at least at high
of D on the Péclet number becomes DyD0 , Pe ln Pe . densities). We also found this to be the case so we ex-
Second, there should be a time for which the dispersion pect to be able to make a meaningful comparison with ex-
is anomalous. Third, rather than becoming irrelevant at periment. At the faces of the simulation box we applied
high Péclet numbers, the molecular diffusion coefficient periodic boundary conditions. However, in all cases we
actually determines the time scale on which the dispersion were careful to eliminate the effects of the periodic bound-
coefficient converges. In this Letter we describe the results ary conditions on the VACF by comparing data for various
of computer simulations of hydrodynamic dispersion in system sizes. At high Péclet numbers it was only possible
random arrays of spheres. We restrict ourselves to the to exclude the influence of the periodic boundary condi-
following questions: How valid is the Taylor hypothesis? tions up until the time required for particles to convect, on
And can we find evidence for the transient anomalous average, one quarter of the box length. This in itself sug-
dispersion predicted by Koch and Brady [12]? gests that correlations in the fluid velocity persist beyond
In our simulations we calculated the velocity autocor- the Brinkman length. The simulations were performed in
relation function (VACF) for tagged particle motion in a a rectangular box with lengths 28 particle diameters (paral-
flowing fluid. The VACF, which we denote Cy std, was lel to the flow direction) and 10 particle diameters (perpen-
calculated along the flow direction and in a frame of ref- dicular to the flow direction), a particle diameter being 5 or
erence moving with the velocity, v̄. We therefore have 9 lattice units. The lattice spacing corresponds to just less
Cy std ­ kfys0d 2 ȳg fystd 2 ȳgl, where ystd is the com- than the Brinkman length for the particles of diameter 5
ponent of the particle velocity parallel to v̄. It is conve- and just more than the Brinkman length for particles of di-
nient to introduce
Rt a time-dependent dispersion coefficient, ameter 9. However, the dispersion coefficients calculated
Dstd ; 0 dt 0 Cy st 0 d. The dispersion coefficient, if it ex- using the different sizes of sphere differed by only a few
ists, is given as the limit of Dstd for t ! `. percent and both gave the same time dependence for Dstd.
The fluid was described by a lattice Boltzmann model. Figure 1 shows a plot of the dispersion coefficient as a
This is a preaveraged version of a lattice gas [14], in which function of the Péclet number. The Péclet number was
the state of the fluid is specified by the average number defined with the velocity equal to the mean fluid velocity
of particles, ni sr, td, with velocity ci , at each site r. The and the length equal to the particle diameter. We have also
time evolution of the distribution functions is described plotted experimental values measured on similar systems
by the discretized analog of the Boltzmann equation [15]. by Pfannkuch [6]. Both the simulations and experiments
Tagged particle motion in the lattice Boltzmann system were conducted at low Reynolds numbers. By using the
can be treated as follows. The probability that a particle moment propagation technique the errors in the VACF,
moves with a velocity ci after a collision, which we denote
as pi srd, is given by ni srdyrsrd, where rsrd is the total
number of particles at the node (here we consider a fluid
which is in a steady state so there is no time dependence).
Once we know pi srd it is possible to calculate the VACF
averaged over all possible tagged particle trajectories using
the moment propagation method [16]. This has been
described elsewhere [17].
Our calculations proceeded as follows. Having gener-
ated an appropriate solid geometry we applied a constant
force density to the lattice Boltzmann fluid. Once the flow
reached a steady state, we calculated Dstd. We tested the
scheme against the few known analytic results and found
it to reproduce accurately, even with a relatively crude lat-
tice, both the time dependence and the value of the dis-
persion coefficient. We also calculated Dstd for a simple
cubic array of spheres where we can compare with the
experimental results of Gunn and Pryce [18]. We found FIG. 1. Hydrodynamic dispersion coefficients, D, as a func-
somewhat better agreement than that reported by Salles tion of the Péclet number, Pe . The dotted line is a guide to
the eye through the (uncorrected) simulation results (crosses).
et al. [19]. A detailed description of this will be given The diamonds are the experimental results of Pfannkuch
elsewhere; here we are primarily interested in the randomly [6]. The solid line shows the simulation results corrected by
packed spheres. Configurations of spheres, with a solid extrapolating the “tail” effect to experimental time scales.

4553
VOLUME 77, NUMBER 22 PHYSICAL REVIEW LETTERS 25 NOVEMBER 1996

for a given configuration, are essentially zero. We do, are two further points to note from Fig. 1. First, at high
however, need to average over a number of different Péclet numbers we are starting to miss a significant part of
configurations of spheres which introduces a statistical the integral by only integrating up to tcut . This means
error. The number of configurations we used depended that our values for D are underestimated. Second, the
on the Péclet number, ranging from a minimum of 5 (at decay of the VACF is not the exponential characteristic
the lowest Péclet numbers) to a maximum of 24 (at the of hydrodynamic dispersion in simple geometries. In
highest). For the values shown in Fig. 1 the error due to particular, the initial rate of decay is rather slow. A
ensemble averaging is at most 5%. A second source of detailed analysis shows that it is in fact proportional
error comes from integrating Cy std up to a time tcut rather to D0 —even at very short times molecular diffusion is
than over all times. At low Péclet numbers (Pe # 80) playing a role.
the VACF has decayed to less than 2% of its initial value We now move to our second question, can we find
by t ­ tcut . The error incurred by truncating the integral evidence for transient anomalous dispersion? In Fig. 3 we
should therefore be small (unless, of course, the integral have plotted the VACF multiplied by t. If the dispersion
diverges). This truncation of the integral becomes more is anomalous the VACF decays as 1yt and the plot should
problematic at high Péclet numbers (see Fig. 2). However, show a plateau. For Péclet numbers Pe $ 480 we see
we should bear in mind that the experimental dispersivities that the VACF is indeed decaying as 1yt. The existence
are also measured in a finite time interval. Although our of this region of anomalous dispersion is in agreement
measurements of D are obtained by truncating the time with the theoretical prediction of Koch and Brady [12,13].
integral at what is, by experimental standards, a short time, However, we find no evidence that this 1yt decay is only
Fig. 1 shows that our numerical values for D agree rather transient. The Koch and Brady prediction is that the 1yt
well with experiment. Our value for the dispersivity at decay persists only for a time proportional to P e1y3 , after
the highest Péclet number we studied (Pe ­ 1280) is 1.35. which the decay should be more rapid. We find that at
This is towards the lower end of the experimental range, Pe ­ 320 the decay is more rapid than 1yt, but appears
a ­ 1.8 6 0.4. The dispersivity, a, is not yet constant to approach a 1yt decay at longer times. However, when
but is still increasing slightly with the Péclet number. increasing the Péclet number by only a factor of 1.5 (to
This is despite the fact that we are well into the regime Pe ­ 480) we find that 1yt decay persists right up until
where dispersion is commonly assumed to be dominated tcut . At higher Péclet numbers we find evidence for only
by convection. 1yt decay. The fact that the anomalous dispersion persists
The VACF at higher Péclet numbers is shown in Fig. 2. once it appears seems more consistent with it being the
The dimensionless time, t, is the time divided by the true asymptotic behavior, rather than a transient effect that
average time required for the tracer to convect a particle scales weakly with the Péclet number.
diameter. If the Taylor hypothesis is correct then the There is a way to exclude the possibility that the disper-
VACF should become independent of D0 at high Péclet sion coefficient is diverging. Our simulations cover a time
numbers; i.e., the curves shown in Fig. 2 should collapse scale that, experimentally, is very short. We can, however,
onto a single curve. On going from Pe ­ 640 to Pe ­ take our data and extrapolate the value of the dispersion
1280 the curves are, however, still distinguishable. There coefficient to experimental time scales, assuming that the

FIG. 3. The velocity autocorrelation function, Cy std, at vari-

FIG. 2. The velocity autocorrelation function, Cy std, at vari- ous (high) Péclet numbers, multiplied by the reduced time t
ous (high) Péclet numbers, plotted as a function of the reduced and plotted as a function of t. A plateau in this plot means
time t. that Cy std is decaying as 1yt.

4554
VOLUME 77, NUMBER 22 PHYSICAL REVIEW LETTERS 25 NOVEMBER 1996

1yt decay persists. If this procedure leads to an incon- is part of the scientific program of FOM and is supported
sistency, we could conclude that the 1yt decay must dis- by the Nederlandse Organisatie voor Wetenschappelijk
appear on longer time scales than those accessible in the Onderzoek (NWO). Computer time on the CRAY-C98y
simulation. In the experiments performed by Pfannkuch 4256 at SARA was made available by the Stichting Na-
[6] the (dimensionless) cutoff time, tcut , varied from 400 tionale Computer Faciliteiten (Foundation for National
to 1500. Our extrapolated values of D for tcut ­ 400 Computing Facilities).
are also shown in Fig. 1. The important point to note
is that the assumption that the 1yt decay persists (and
hence that D diverges) leads to estimates for D that are,
at least, consistent with experiment (and even marginally
better than if we simply ignore the 1yt decay beyond tcut ). [1] C. S. Slichter, Water Supply Papers 14R, U.S. Geological
Our data suggest that, in the range 320 # Pe # 1280, the Survey (1905).
tail effect leads to an increase of about 25% in the ap- [2] N. S. Martys, Phys. Rev. E 50, 335 (1994).
parent dispersivity on going from tcut ­ 400 to tcut ­ [3] G. I. Taylor, Proc. R. Soc. London A 219, 186 (1953).
1500. Pfannkuch’s data in this range gives a ­ 1.5 6 0.2 [4] J. J. Fried and M. A. Combarnous, Adv. Hydr. Sci. 7, 169
(1971).
(tcut ­ 400) and 2.2 6 0.1 (tcut ­ 1500), an increase of
[5] G. I. Taylor, Proc. London Math. Soc. 120, 196 (1921).
some 50%. This is subject to appreciable statistical error, [6] H. D. Pfannkuch, Rev. Inst. Fr. Petrol. 18, 215 (1963).
although not as much as Fig. 1 might, at first sight, suggest [7] N. E. Rifaı̈, W. J. Kaufman, and D. K. Todd, Sanitary Eng.
(a is obtained by averaging over all data points within the Report No. 3, I.E.R. series 90, Berkley, CA (1956).
specified range of Péclet numbers; this reduces the error in [8] D. A. Chin, J. Hydr. Eng. 112, 591 (1985).
a as compared to the error in the individual data points). [9] N. Han, J. Bhakta, and R. G. Carbonell, AIChE J. 31, 277
We can at least conclude that this trend is in agreement (1985).
with the simulation. [10] E. Charlaix, J. P. Hulin, and T. J. Plona, Phys. Fluids 30,
In conclusion, our simulations suggest that at high Péclet 1690 (1987).
numbers there is at least a region of anomalous dispersion. [11] P. G. Saffman, J. Fluid Mech. 6, 321 (1959).
This means that the Taylor hypothesis is incorrect and that [12] D. L. Koch and J. F. Brady, J. Fluid Mech. 154, 399
(1985).
the modified relation predicted by Koch and Brady is more
[13] D. L. Koch and J. F. Brady, Chem. Eng. Sci. 42, 1377
appropriate. However, our results are more consistent with (1987).
the conclusion that the hydrodynamic dispersion coeffi- [14] U. Frisch, B. Hasslacher, and Y. Pomeau, Phys. Rev. Lett.
cient actually diverges. A detailed experimental study of 56, 1505 (1986).
the limiting behavior of Dstd would clearly be desirable. [15] A. J. C. Ladd, J. Fluid Mech. 271, 285 (1994).
We stress that there is no reason why a dispersion coef- [16] D. Frenkel and M. H. Ernst, Phys. Rev. Lett. 63, 2165
ficient must exist—there are several limiting cases where (1989).
it does not (in a two-dimensional fluid [17], in a turbulent [17] C. P. Lowe and D. Frenkel, Physica (Amsterdam) 220A,
fluid [20], and even, in the dilute limit, for hydrodynamic 251 (1995).
dispersion itself [12]). Moreover, there is experimental [18] D. J. Gunn and C. Pryce, Trans. Inst. Chem. Eng. 21, 583
evidence that the dispersion coefficient often does not (1969).
[19] J. Salles, J.-F. Thovert, R. Delannay, L. Prevors, J.-L.
reach a constant in time. Explanations that have been put
Auriault, and P. M. Adler, Phys. Fluids A 5, 2348 (1993).
forward to account for this observation are that the medium [20] J. K. Klafter, M. F. Shlesinger, and G. Zumofen, Phys.
has a fractal structure [21] or that the dispersion mecha- Today 49, No. 2, 33 (1996).
nism is locally nonlinear [22]. Our results provide a much [21] T. A. Hewett, presented at the 61st Annual Technical Con-
simpler explanation. ference of Petroleum Engineers, New Orleans, Louisiana,
We thank Professor D. Koch and Ignacio Pagonabar- 1996.
raga for their comments. The work of the FOM Institute [22] J. P. Pascal, Physica (Amsterdam) 223A, 99 (1996).

4555

View publication stats