Dispersion and Reservoir Heterogeneity

Atul Arya, * SPE, U. of Texas

Tom A. Hewett, SPE, Chevron Oil Field Research Co.
Ronald G. Larson, AT&T Bell Laboratories
Larry W. Lake, SPE, U. of Texas

Summary. Macroscopic dispersion is the mixing, on the scale of several hundreds of grain diameters, at a point in a permeable
medium that is free of boundary effects. Megascopic dispersion is the one-dimensional (lD) dispersion derived by averaging across
an entire cross section. This work investigates how both dispersions vary with heterogeneity, aspect ratio, diffusion coefficient, and
autocorrelation. The theoretical results are compared to existing field and laboratory data and to existing theories for limiting cases.
The degree of autocorrelation in the medium determines whether or not megascopic dispersivity (dispersion coefficient divided
by velocity) is uniquely defined. Large correlation distances (with respect to the medium dimensions) imply a dispersivity that
grows with distance traveled. Small correlation distances imply a dispersivity that is eventually stabilized at some constant value.
This value is related to the heterogeneity of the medium. On the field scale, diffusion is insignificant, but on a laboratory scale, it
can stabilize the dispersivity even if the medium is correlated. Macroscopic dispersivity is sensitive to diffusion in both the
laboratory and field scale. It is smaller than or equal to megascopic dispersivity, also in conformance with experimental data, and
comparable to laboratory-measured dispersivity.

Introduction
Dispersion is mixing caused by variations (heterogeneity) in the recovering mixture in the pores. This happens in developed miscible
velocity within each flow channel and from one channel to another. floods or through the generation of an optimal salinity in micellar/
Molecular diffusion is the transport of mass because of spatial con- polymer floods.
centration differences. Dispersion and diffusion in permeable media Megascopic and macroscopic do not necessarily coincide with
play an important role in miscible displacement, where channeling field and laboratory scales. The dispersivity measured in a labora-
and/or fingering of the displacing fluid occurs. We examine the tory experiment is megascopic, but it differs from a field-scale
interrelationship between heterogeneity and diffusion in this paper. megascopic dispersivity because of a smaller heterogeneity, a gener-
ally smaller aspect ratio, and a greater influence oflateral bounda-
Purpose and Scope ries. We shall see below that megascopic dispersivities measured
In this paper, mixing resulting from macroscopic variations in the on laboratory displacements are comparable to the macroscopic dis-
permeability of the medium is of primary interest. These heter- persivity in a field-scale displacement.
ogeneities calise fluctuations in the velocities of individual fluid ele- Both macroscopic and megascopic dispersions in a permeable
ments. To focus on the mixing resulting from heterogeneity and medium result from contrasts in hydraulic conductivity. Generally,
diffusion, the process under consideration is one in which one fluid the larger the contrast between the elements that make up the medi-
displaces another fully miscible fluid. The investigative tool is nu- um, the larger the dispersion. In these cases, mixing is produced
merical simulation of first-contact-miscible, equal-density, constant- by permeable-medium nonidealities that are responsible for changes
mobility displacements in two-dimensional (2D), randomly heter- in the direction and velocity of flow. For most geologic systems
ogeneous (RH) flow fields. of interest, the most significant dispersion will be generated by this
A primary purpose of this work is to investigate the behavior mechanism.
of longitudinal dispersivity in field-scale miscible displacements.
We place special emphasis on analyzing systems with large heter- Classifying Mixing
ogeneity (Dykstra-Parsons 1 coefficient, VDP =0.6 to 0.8), which
We classify mixing by the time behavior of the dimensionless mixing
is the norm in oilfield cores. 2 These values are considerably higher
or transition zone, IlxD' IlxD is the distance, normalized by L,
than those generally investigated in the past. 3 We also investigate
over which the cross-sectionally averaged displacing fluid concen-
the effects of diffusion and system aspect ratios [system length
tration changes a fixed amount (usually from 0.1 to 0.9). A gener-
parallel to bulk flow divided by length perpendicular to flow (Lib)]
al form of the lD linear material-balance equation can be used to
on the behavior of longitudinal dispersivity. classify different types of mixing zones:
Definitions
Macroscopic dispersion is the mlxmg, on the scale of several ac ac a2 c
l/J-+uF(C)- -D(C)- =S(C), ................. (1)
hundreds of grain diameters, at a point in a permeable medium that at ax ax 2
is free of boundary effects. Megascopic dispersion is the lD dis-
persion derived by averaging across an entire cross section.
The megascopic description is typically the size of a gridblock where
in a numerical reservoir simulation model 4 and contains many F( C) = flux function,
macroscopic elements. The megascopic dispersivity, Cl.ME, is im- D(C) = diffusion function, and
portant in field-scale simulation of EOR processes because it con- S(C) = capacitance or source function.
trols volumetric sweep efficiency within a gridblock.
A macroscopic description describes a permeable medium in terms For illustration, we take D(C) constant = l/JKe, and S(C) =k'C,
of average properties and their variations at scales much larger than where k' is a rate constant. In dimensionless form, Eq. 1 becomes
pores. 5 In EOR, the macroscopic dispersivity, Cl.MA, controls oil
recovery by determining the rate of formation of an effective oil-
k 2
'Now at Intera Technologies Inc.
!.£+F(C) ac _(l/J e)a C = k'L C, ............... (2)
Copyright 1988 Society of Petroleum Engineers atD aXD uL aXD u

SPE Reservoir Engineering, February 1988 139

 ......
-QI
Ie.
z
......
_
t!)
z .....
>< Q
fI ux Induced
f1 ux Induced
-
J:~
(I) .....
x
", ...\ .............. '"., ...................... ,... ,..... . Di.pe,..ive
111111111111

(I):t:
UJI-
-'t!)
zz
QUJ
Capacitive
- -'
(l)UJ
11'1'lllfll"""'IIIIIIIIII""lltlllllllll'l

Zz
~Q
_ N
Q Capacitive

OU1ENSIONLESS TIME (to)

DIMENSIONLESS TIME (to)

Fig. 1A-Classifying mixing by Axo behavior. Fig. 1B-Classifying mixing by N Pe1 behavior.

where tion to immiscible displacements, or in any displacement that is

locally bypassing oil (channeling or viscous fingering). Capacitive
uL mixing occurs because of mass transfer from locally stagnant regions
- - = N pe , in the medium or because of nonequilibrium chemical reaction. Dis-
<J>Ke persive mixing is caused by diffusion or mechanical mixing in mis-
cible displacemcnts or by capillary pressure in immiscible
k'L displacements.
- = NOa , In oil reservoirs, mixing between fluids is caused by a variety
u of mechanisms. For example, in a miscible flood, mixing can be
a result of mass transfer by extraction, vaporization, condens&tion,
dissolution, and channeling. Therefore, it is difficult to categorize
mixing in field-scale displacements exclusively into one of the three
types described. A combination of these causes 1:J.xD and disper-
sivity to behave differently than described earlier. However, these
basic mechanisms help understand behavior in more complex het-
erogeneous reservoirs.

Experimental Measurement of Dispersivities

N pe = Peclet number,
Field-Scale Dispersivities. The best summary of measured mega-
NOa = Damkohler number, scopic dispersivity, OIME' is by Lallemand-Barres and Peaudccerf.p
tD = dimensionless time, and Fig. 2 shows all the field data collected in this study and the meas-
XD = dimensionless distance. ured OIME's for eight different formation types plotted vs. the
length over which the dispersivities were measured. There is. con-
In what follows, we use OIME and NPel interchangeably. They are siderable scatter in this log-log plot in dispersivity measured at the
related by same length. The trend clearly points, however, toward increasing
OIME with measurement distance. Pickens and Grisak 7 also present
a comprehensive summary of OIME' These values, also plotted in
--"'--"'-- ............................... (3) Fig. 2, follow the same trend.
L uL The trend in Fig. 2 can bc qualitatively explained by the obser-
vation that the "scale" of heterogeneity increases as the sample
A flux-induced mixing zonc occurs when NPel =NOa =0. In this volume increases. The trend in OIME may also be a consequence
case, 1:J.xD is a linear function of tD and is independent of L. In of streamline effects. If streamlines diverge or converge in some
a capacitive mixing zone, NPel =0 and F is constant. Hcre, 1:J.xD manner, ID model analyses may yield an apparent decrease or in-
approaches constancy with respect to tD, and its magnitude is in- crease, respectively, in the interpreted value of megascopic dis-
versely proportional to L. In a dispersive mixing zone, NDa =0 and persivity. This is unlikely to be a factor in Fig. 2 because most
F is constant. 1:J.xD grows in proportion to the square root of time. measurements were point samples from sources producing at low
1:J.xD also varies with L -'/2. Fig. IA shows schematically the var- rate and hence sampled only a few streamlines.
iation of mixing zones with respect to tD'
We can also use dispersivity itself as an indicator of the type of Laboratory Corefloods. One important criterion for selecting ex-
mixing occurring during a displacement-i.e., we match the solu- perimental data to present here was that the mobility ratio of the
tion of Eq. 2 with NDa =0 and F= I to the actual response of the floods be I or slightly less than I. The da~Jrom numerous miscible
system. When we do, this dispersivity is no lopger constant in gener- flood experiments are also plotted in F~. 2.8
al. In the case of flux-induced mixing, dispersivity grows linearly The trend including the laboratory data is also linear but with
with tD, and this growth is independent of L. In a purely capacitive a larger slope (al\ valucs for distances less than 5 m [16 ft] are lab-
mixing zone, dispersivity decreases with tD and reaches an asymp- oratory values). The laboratory data thcmselves do not show a trend.
totic lower limit at large t D' It depends on L, with small systems The log-log least-squares fits to the Fig. 2 arc for field data,
having larger dispersivities. In a dispersive mixing zone, the dis-
persivity is invariant with respect to time. Fig. 18 shows the varia-
tion of dispersivity with tD for the three types of mixing zones. OIME=0.229Lo.755, ............................... (4a)
Except for the purely dispersive case, the dispersivities should be
called apparent because they are not actu&lly a measure of dispersive and for all data,
mixing.
The three types of mixing are manifestations of a particular mech-
anism. Flux-induced mixing occurs in the Buckley-Leverett solu- O'ME=O.044LI.13 . ................................ (4b)

140 SPE Reservoir Engineering, February 1988

 We will show later that displacements in an RH permeable medi- 1000.000
um are generally non-Fickian (DiME is a function of time or mean
D
travel distance).

Description of Numerical Model. The numerical model is an in- 100.000

- ~/' - ,/
compressible, 2D, finite-difference simulator that solves the
convection-diffusion (CD) equation, along with a continuity equa-
tion, in the x-yplane: 10.000
c
--.--
• c /
/

ac a (uxC-</>D o-ac) +-
</>-+-
at ax
a (uvC-</>Do-ac) =0, ... (5)
ax ay' ay 5
in
a:
UJ
I-
UJ

U.I
1.000
t/-
-Ir.
- !.II;~r- ,.
c
10

au, au,.
I:
'1::1 v/~J c- - lal1emand-Barres
0.100
- + - ' =0, PIckens lie Grlsak
.................................. (6a)
V
D
ax ay .-/
and •
0.0101 ;<!
c
• Lab data
- All data

.
kj ap .
Uj=---' }=x,y . .............................. (6b)
Jl aj
0.001
I· (
0.1
•
•••
1.0 10.0 100.0
.... Field data

1000.0 10000.
The model uses the following initial and boundary conditions: DISTANCE (METERS)
C(X,O) =0, ....................................... (7) Fig. 2-Field and laboratory dispersivities.

ac
-(L,t)=O, ...................................... (8) where C is the concentration of the displacing agent injected at unit
ax concentration and

and 2 roo 2
erfc(x)=--J e- U du
";;x
(UXC-</>D O ac) =u ......... .................... (9)
ax x=o is the complimentary error function. Eq. 10 is semi-infinite because
the condition (Eq. 8) has been replaced by C(oo,t)~O.
In thcse equations (Eqs. 4 through 9), u is the cross-sectional With the second term omitted, the solution (Eq. 10) has been wide-
average of uX ' These equations are reduced into dimensionless ly used for analyzing finite corefloods and tracer tests. II This ap-
form, and a finite-difference analog for the partial-differential equa- proximation is good only when tD is large and/or Nrel is small,
tions is developed. 8 Eq. 5 contains only a diffusion coefficient, or under circumstances where the inlet boundary appears as though
Do, which is the effective binary molecular diffusion coefficient it were a large distance from the displacing front for most of the
divided by the product of the porosity and formation resistivity fac- flood.
tor. 9 Upon nondimensioning, this becomes the dimensionless In the system, however, modeled boundary effects significantly
diffusion coefficient D=</>Do/uL. Except where noted, all the dis- influence the concentration profiles, even for homogeneous-
persion calculated is a result of permeability variations in the flow permeability fields. In such a case, we find that using the semi-
field. infinite solution to match the numerical simulations gives erroneous
One of the major problems in deriving physical dispersion caused results (e.g., time-dependent dispersivities for homogeneous sys-
by a random-permeability field is numerical dispersion. An obvi- tems).8 Therefore, to account for a finite system, we cannot use
ous technique to reduce numerical dispersion is to reduce the size Eq. 10 to infer Nrel.
of a gridblock in a given system and thcn to extrapolate the results An analytic solution for systems of finite length was presented
to zero block size. For thc same-permeability field, dispersivities by Brenncr. 12 However, this solution converges very slowly. To
were calculated at different gridblock sizes. The values decrease overcomc thesc drawbacks, the I D CD equation with the condi-
in magnitude with decreasing gridblock size or increasing number tions ofEqs. 7 through 9 was solved with Laplace transforms. The
of gridpoints (nGP)' These dispersivities were plotted vs. IInGP' solution in Laplace space was inverted with Stehfest's algorithm. 13
and the straight line was extrapolated to zero, IInGr->O (or By comparing analytic and numerical results on homogeneous dis-
nGP~oo). This extrapolated value is the true megascopic disper- placements, we found that the numerical inversion is accurate to
sion without a numerical component. For very heterogeneous sys- within 3 %. This solution was used to match the simulation results
tems (VDP~0.6), the numerical dispersion is negligible compared from the finite-difference model from which we derived values for
with megascopic dispersion even when nGP is rather small. For ITME'
increasingly smaller blocks, the time step size must be reduced to Many of our runs will prove to have a time-varying dispersivity.
maintain stability; hence, the procedure also corrects for timestep Arya 8 has shown that both Eq. 10 and the Laplace-transform in-
size. All results presented here have been corrected in this fashion. verted solution are valid if the dispersivity is replaced by the time-
A detailed description of these simulator runs can be found averaged value
elsewhere. 8

Analytical Solution to CD Equation I =~=~rtD~dt=~rtD_l_dt,

__
J ~ J ~ .............. (11)
The ID version of Eq. 5 has an analytic solution for a semi-infinite N pc L tD 0 L tD 0 N Pe
system 10:
where the overbar denotes a time-averaged quantity. Of course,
if IT is constant, then a=Di and NPe=N pc .
The finite different model calculates concentrations as a func-
tion of x, y, and t. These values are space-averaged in the y direc-
tion (transverse to flow) to obtain a concentration profile in the x
................................... (10) direction. In a typical run, concentration profiles are obtained at

SPE Reservoir Engineering, February 1988 141

 Fig. 3-lsometric view of permeability, V DP = O.S.

every 0.1 tD' These profiles are then least-squares m,!!ched with where
the analytic solution described to obtain a value for IIN Pe at e!l:ch n = number of samples,
time. The final result from one simulation run consists of IIN Pe x(i),x(i +f) = sample values separated by lag spacing e, and
values for successive values of tD' x = estimated mean (assumed stationary) of the
population.
Permeability Assignment
Many schemes are available to assign random-permeability p(R) has the maximum value of 1.0, which is the correlation of each
fields. 14.15 The method used in this study was first suggested by datum with itself. For f> 0, the autocorrelation p(C) generally
HeUer. 16 In this scheme, a set nHP of random locations, "Heller decreases with increasing P, ideally approaching zero at a finite sepa-
points," in the 2D flow domain is chosen. For each Heller point, ration.
a permeability is chosen randomly from a log-normal distribution. The shape of the autocorrelation curve depends on the spatial rela-
Permeability values for each of nGBC~nHP) gridblocks are then tionship of the data. 18 Rapid decay of p(f) to zero within the first
assigned by interpolating the permeability values of nearby Heller lag signifies spatial independence, indicating variation on a scale
points. The closer the center of a gridblock lies to a given HeUer smaller than sampling distance. Gradually decreasing p(R) that ap-
point, the more closely its value reflects the nearby Heller point proaches zero asymptotically (within the computed lags) indicates
value. Using this procedure, we obtain a random-permeability field long-range spatial dependence.
that becomes smooth as the grid is refined with nHP held fixed. p(P) is interpreted quantitatively by means of an integral scale,
which gives an estimate of the mean separation distance over which
Description and Significance of a variable is correlated with itself. 19 We will use these ideas later.
Statistical Parameters Fig. 3 shows an isometric plot of log-permeability for a spatially
correlated permeability field with Lib = 1.
Heterogeneity Measure. One important statistical parameter used
to characterize heterogeneity is the Dykstra-Parsons coefficient I
Results
defined as
Megascopic Dispersion. In this section we show only a sample
of a large number of computer runs. 8 The most important factors
kO.5 -k(f
VDP = , ................................. (12) that influence CiME are heterogeneity, VDP ; autocorrelation; and
k05 aspect ratio, Lib. Diffusion is less important to CiME because the
long length scales cause D to be essentially zero. Diffusion may
where affect Ci MA, however.
For any fixed VDP, D, and Lib, we can obtain a curve showing
ko.5 = median permeability (value with 50% frequency of a relationship between ii ME and t D that has been corrected for nu-
ocq:lrrence) and merical truncation error. Fig. 4 shows a plot of iiME (plotted as
k(f = permeability at 84.1 % of the cumulative sample. IIN pe ) vs. tD for several VDP values with Llb= I. In all runs, the
flow-field structure remains the same, and only the magnitude of
heterogeneity is changed. iiME at any fixed time increases with
VDP varies between zero and one, with a completely uniform higher VDP, indicating an increase in heterogeneity. The variation
system having a value of zero. In typical oilfield cores, VDP usuaUy of iiME with tD for a fixed VDP is more complex. For smaller VDP '
falls between 0.6 and 0.8. 2 ii ME is essentially constant; however, at higher VDP, dispersivity
increases with tD' The increase for tD>0.8 is attributed to outlet
Autocorrelation. The HeUer procedure generates a field that, boundary effects that cannot be included in the simulations when
although still random, is spatially correlated. 17 The spatial depend- D=O. The early effect, however, results becausc a finite distance
ence ofa regionalized variable may be quantified by an autocorre- is required for iiME to stabilize, which will be discussed later.
lation function. 15 The autocorrelation function, p(£), is The growth of iiME with tD also depends on the degree of auto-
correlation in the flow field. Fig. 5 shows the autocorrelation func-
p(£) = C(£)/C(O), .................................. (13) tion of permeability for three runs with nGP = 1,600 gridpoints and
three nHP = 100, 400, and 1,600. For square grids, the integral
scale (the distance where p= lie) is approximately L' niiJf2; i.c.,
where C(f) = estimated autocovariance of samples separated by lag nHP =nGP = 1,600 corresponds to the smallest correlation possible.
spacing fand C(O) is the estimated variance of the sample set. C(f) The horizontal axis in Fig. 5 expresses the lag distance as a multiple
is estimated by of the gridblock spacing.
Fig. 6 shows the variation of iiME with tD for the three cases
I n-f in Fig. 5. With little or no spatial correlation (nHP large), iiME is
C(f)=- L; [x(i)-x][xCiH)-x], ................ (14) essentially constant. If correlation exists, iiME grows with time and
n-f i=1 ultimately stabilizes at much higher values. It will bc shown later
142 SPp Reservoir Engineering, February 1988
 0
0

"HP = 1600 "GP =1600

"HP =400 "GP = 1600
e 0

'"ci
0,----------------------------------------, "HP = 100 "GP = 1600

II)
o .
0

IZ
.,co. ci
VOP -0.8 ZO
S!
<
~

...J
a: .. ""'0
a: ...
~ tq a:
00
~ 0 u
0
Z ~

..,
~ ~
<
...J
..,U ...o 0

'"
co. 0
.., ci
..,'"a:
>
Z
0
0
'"o 0
ci Vnr = 0.5
"OP Q 0-" 0

g1-__~~~==~~~~~;;;V~O~P~=~O~.J~::~~~
N
0
'0. 00 2.00 4.00 6.00 8.00 10.00
"b. 00 0.20 O. 40 O. 80 0.80 1. 00 LAG DISTANCE
OIMENSIONLESS TitlE to
Fig. 5-Autocorrelation for permeability, V op = 0.6 and
Fig. 4-to vs. N p.' as a function of Vop. LIb = 1.

that these results are consistent with theoretical findings. and produces the effect of a local layering of permeabilities. As
Schwartz 14 has previously observed time-dependent dispersivities. the value of Lib is increased, the variation of dispersivity approaches
The system aspect ratio also influences the variation of iiME with the linear vl\riation seen in layered systems with no crossflow, as
tD' In all previous cases, this ratio was one. Fig. 7 shows the calculated by Mercado. 2o Isometric concentration fields for two
growth of ii ME with tD for three cases with Lib = 1, 5, lind 10. Dis- aspect ratios are shown in Figs. 8A and 8B. In the nondimensional
persivity increases with t D for higher aspect ratios and the rate of coordinates used for these plots, the layering effects of high aspect
increase is greater at higher Lib values. We find that this behavior ratios are not obvious.
also agrees with the behavior of the integral scale, but the latter It is also interesting to consider the pressure fields correspond-
cannot be simply related to nHP when Lib 1. 8 * ing to the above simulations. Figs. 9A and 9B show that as Lib
In these calculations, the aspect ratio of the entire permeability increases, the pressure contours become more uniform in the y direc-
field was changed. This means that the scale of variations in the tion and the pressljfe gradients transverse to the flow direction ap-
y direction was changed relative to that in the x direction by a factor proach zero. A zero pressure gradient transverse to the flow
of b/L and the correlations are anisotropic. This procedure leads direction is referred to as vertical equilibrium. The assumption of
to permeability variations that are elongated in the direction of flow vertical equilibrium is valid in uniform media and layered systems

e 0
on
ci ci
=1600 =1600 Lib =10
.
"HP "GP
"HP =40q "GP =1600 Lib = 5
0 ...
0

7" cO
IZ
.
Co
"HP = 100 "GP = 1600
7" ... ci
'za:
Co
Lib =1

...a: ... ....,

=
~ ci
0
=ci....
!i
Z
0

Z
...
~ ~

""
...co. ...
...J ....J
U ...., 0
U
0
a.. '"
... 0
rtl
...., ci
I/)

>
...
a: a:
....,
>
:!i !!: 0
'"
0
ci
0

0 o
0 o
00 . 00 0.20 0.40 O.SO 0.80 1.00 9>.00 0.20 0.40 0.50 Q.80 1.00
DIMENSIONLESS TIME to DIMENSIONlESS TIME to

Fig. 6-Effect of correlation length on N p- . \ V op =0.6 and Fig. 7-Effec;t of aspect ratio (Lib) on N p-.'. Vop =0.6.
LIb =1.

SPE Reservoir Engineering, February 1988 143

 Q
U
Z
o
~ o.
a:
....z
'"~
u

Fig. SA-Isometric view of concentration, V DP = 0.6 and LIb = 1.

Q'
U
Z
o
~
a:
lE
....
u

u
~

Fig. SB-Isometric view of concentration, V DP = 0.6 and LIb = 10.

y-DIRECTION BLOCK NlI1BER y-DIRECTION BLOCK NlI1BER

~. 00 8.00 HI. 00 24.00 32.00 40.00 ~.oo 8.00 18.00 24.00 32.00 40.00
.;;-------~------~------~------~----_4 .;

0 0
0 0
.; .;
..,a:
a:l
..,a:
a:l

z~ 0
0 z~ 0
0

u
0
:..:
..J
a:l
Z
00 :..:
u
0
..J
a:l
Z
-
00

52 0
52 0

...a:.... ...• ...a:u.... ...•

0 0
u

Q QI
I
x

..
0
0

n
..
0
0

o 0
o 0

~~~----------------------------------~ ....;
Fig. 9A-Pressure contours, V DP = 0.6 and LIb = 1. Fig. 9B-Pressure contours, V DP =0.6 and LIb =10.

144 SPE Reservoir Engineering, February 1988

 o o
o o
~:.-------------------------~--~~ ~:~------------------------------~

...
0
...o
0 o

..
0 Q
U U
0 o
Z Z ..
S1 0 90
I-
<
a: !<a:
I-
zU.J
0
!£.... 0
U
Z
0
•
0
U.
Z
00
.
U U

...
0

0
..o
o

o o
o
~oi.-o-o----~~~~~~~~~~~~~~~~o
o. 10 0.20 0.30 0.040
.• o °o4.-o-0----~--~~-r------~------~------0~.·0
J 0.10 0.20 0.30 0.040 ~
DIMENSIONLESS TIME to OlMEN510lilESS TIME to

Fig. 10A-Concentration history at Xo =0.5, VoP =0.6, Fig-. 10B-Concentration history at x 0 = 0.5, V OP = 0.6,
LIb = 10, and D = 0.00. Llb=10, and D=0.05.

with a large aspect ratio. 21.22 The results presented here indicate o
o
that it is also valid in RH media with Lib'$> 1. Such behavior im- ~~------------------------------~
plies rapid relaxation of pressure perturbations transverse to the
bulk flow.

Macroscopic Dispersion. Fig. lOA through lOe show the con-

...o
o
centration history at various fixed points across the medium for three
cases with different amounts of diffusion. When D=O, the curves U
Q

are widely separated, showing that the arrival of any concentra- z o..
tion varies from point to point (in the transverse direction), even o·
-0
at a fixed location in the flow direction. This can also pe interpreted !<
as channeling of the flood front as a result of heterogeneity. In Fig. ~
.... 0
lOB, D=O.05, whicp tends to make the flood front more uniform u.
z.
than before. The cause of this is transverse mixing, which diffu- go
sion enhances. In this case, all history curves are close to each other,
indicating a stable displacement without channeling. If diffusion
is increased further (FIg. 1OC), all curves collapse into a single
curve. In this case, the concentration profiles are smooth, inverse,
.o
o

S-shaped curves, wpich indicates that C'iMA does not vary with tillle.
This is the result predicted by analysis of Taylor's dispersion re-
o
sulting from crossflow in layered media. 22,23 o
°0:4.-0-0--~~-------r-------r-------r------~o.·0
0.10 0.20 0.30 0.40 ..
Analyses of Resul.s DIMENSIONLESS TIME to
Megascopic Dispersion. The diffusion-free results may be analyzed
by methods used in turbulence theory because turbulent flow and Fig. 1DC-Concentration history at x 0 = 0.5, V OP = 0.6,
Llb=10, and D=D.25.
creeping flow in f(H permeable media both take place in stochastic-
velocity fields. For pedagogical reasons, We base the following on
the older theory of Taylor, 24 even though more recent results are
available. 25 The second identification means that we are equating multiple reali-
The distribution around the mean of a large ensemble of particles zations of ID flows to flow in RH media. These yield the following
released at the inlet of a RH permeaple medium will have the fol- relationship petween aME and mean distance traveled:
lQwing variance after they have traveled a mean distance, X.

a; =2C~ rr
o Q
p(~)d~dx, ............................ (15)

where C v = coefficient of variation of the interstitial velocity (stan- From this, the magnitude of aME depends on C v ' but whether or
dard deviation of the velocity distribution divided by the mean) and not dispersivity is scale-independent depends on the nature of p_
p=autocorrelation function of the velocity field. Because the flow A very common model for p is the exponential function
is incompressible, C v and p are independent of time. Eq. 15 is free
of scale factors, except those that might be part of p. The corre- p(~)=e-~/h, ..................................... (18)
spondence between the quantities ip Eq. 15 and those used previ-
ously is
where A= mean cQrrelation length (the integral scale discussed earli-
x=tDL, 2tDINPe=a;IL2 . .......................... (16) er). Eq. 18 is a purely empirical representation, 19 but it does seem

SPE Reservoir Engineering. February 1988 145

 100.0
00
0 0
0

/
•
(> /.
10,0
/
<:>

.'
(>
0 (>
Q
.....
~

1.0 <:>
o <:> 0
y/
....2........Q...,Q..~ ............
0

0.1
0.001 0.010 0.100 1.000 10.000 100.000
U.dp/Do

Fig. 12-Comparison of experimental vs. numerical disper-

sion coefficients.

The latter result is independent of f.. and is the same result as

6 6 6 Mercado's20 for layered media. This indicates that for x ~ f.. the
o
o velocity field is highly correlated, just as it is in a layered system
°o+.-o-o-----o~.-2-0-----0~.-4-0-----0~.-8-0-----0'.-8-0-----4'.00 with no crossflow. The latter observation suggests a classification
DIMENSIONLESS TINE to for heterogeneous permeable media: if f.. ~ L, the medium is diffu-
sive in mixing; if f.. - L, the medium is effectively layered; and if
Fig. 11-Comparison of analytic vs. numerical results. f.. ~ L, the medium is strictly layered,

Alternative Distributions. Eq, 20 suggests that once x becomes

to fit a large variety of stochastic fields (see Fig. 5). Substituting greater than f.., OiME should stabilize, No such stabilization is ap-
Eq. 18 into Eq. 17 gives parent in Fig, 2 despite a five-fold change of a factor of 10, Here
we are associating x in a "generic" RH medium with the horizontal
axis of Fig. 2. The lack of stabilization implies that f.. in the data
is larger than about 1000 m [3,280 ftl, Inasmuch as such extreme
correlation lengths seem unlikely, we investigate an alternative cor-
Eq. 19 was used to generate the solid curves in Fig. II. In each relation model.
case, C v was independently estimated from the standard deviation An alternative scale-free correlation model is that for fractional
and mean of the assigned x-direction transmissibilities. (Because Brownian motion introduced by Mandelbrot and Van Ness, 26 For
transverse pressure gradients are small, this is essentially the same this distribution, the variance of the displacement obeys
as C" based on velocity.) In a similar manner, f.. was derived by
fitting Eq. 17 to the mean correlation function, which was calcu- a~=CDx(4-2Dj), . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (23)
lated from Eq. 13, wherex(i) is the x transmissibility at x-coordinate
i. The satisfactory agreement at moderate to low values of VDP where CD is a constant and Dj is the fractal dimension. In view
suggests that (l) the truncation-error removal procedure is satis- of Eq. 20, we have
factory and (2) the difference between our 2D simulations and the
multiple realizations of a stochastic ID flow assumed by Taylor CXME=(2-Dj )CD X(3-2Dj ). " ....................... (24)
is small.
When VDP=0.6, there is some deviation from Eq. 19 (see Fig. Totally uncorrelated variations produce classic Brownian motion
II). We attribute these to differences between the actual correla- with D j = 1,5 (compare with Eq. 21), The field data in Fig, 2 are
tion structure in the medium and Eq. 18 (note from Fig. 5 that p best fit with a value of Dj = 1.12, suggesting long-range correla-
can actually be negative) and to the absence of multiple realiza- tions. It is a characteristic of fractal distributions that finite subsets
tions in the simulations. 3,8, 15 of a distribution appear correlated at distances on the same order
All simulation results reported thus tar have been in terms of the as the sampled subset. 26,27
time-averaged dispersivity, CXME, or N Pe ' In view of Eqs. 11, 15, Fractal dimensions of 2.5 to 2.78 have been observed in pores. *
and 16, the instantaneous value of megascopic dispersivity behaves From geophysical records of river fluctuations, fractal dimensions
as vary between I. 1 and 2.0. 28 The connection between these obser-
vations and the fractal behavior of dispersivity, however, has not
.............. ... . ... (20) been established .

Macroscopic Dispersivity. In this section, we attempt to estimate

The satisfactory agreement between Eq. 19 and Fig, 11 suggests values for macroscopic, or pore-level, dispersivity, cxMA- To do
some limiting cases for the behavior of CXME' If f.. is much smaller this, we analyze each of the curves in Figs. lOA through lOe in
than the largest mean distance traveled, L (the well spacing or the the manner suggested by Brigham. II This technique should give
length of an average streamline), the dispersivity is constant at results that are largely independent ofthe streamline being sampled
by the particular curve and by the time level. The NPel =CXMAIL
CXME=c3f.. . .... , . . . . . . . . . . . , . . . . . . . . . . . . . . . . . . . . (21) based on individual curves varied a little, but not systematically,
We report the spatially averaged Nre
l value for each run here.
At the other extreme, Eq. 20 becomes These NPel were then corrected for truncation error, in this case
by simply subtracting (4ncp) _1/" This technique works nearly as
CXME=c3x . .. , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (22)
~ Katz,
AJ. and Thompson, AH.: "Fractal Sandstone Pores: Implications for Conductivity
Eq. 22 says that CXME should increase linearly with x, and Pore Formation," submitted to Physical Review Letters.

146 SPE Reservoir Engineering, February 1988

well as the more elaborate procedure described previously and is Do = diffusion coefficient
more suited to analyzing concentration histories. The remaining con- DeC) = diffusion function
tribution to NPe I is an estimate of the macroscopic dispersivity. F( C) = flux function
Fig. 12 shows the results of this procedure applied to a series k = permeability
of runs with increasing amounts of diffusion. All runs had ka = permeability at 84. I % of cumulative sample
ncp=400=nHP with VDP =O.6 and Lib = 10. From the results
ko.s = median permeability
above, these runs would have a nonstabilized OI.ME when diffusion
Kr = longitudinal dispersion coefficient
was small, but OI.MA was stable throughout.
L = system dimension parallel to bulk flow (length)
Also shown in Fig. 12 are coreflooding results reported by perkins
and Johnston,9 which plot a dimensionless dispersion coefficient n = number of samples
on the vertical axis and a dimensionless rate on the horizontal. Both nCB = number of gridblocks
axes are logarithmic. The correspondence between these parame- ncp = number of gridpoints
ters and those used here are nHP = number of Heller points
N Da = Damkohler number
Kr ud NPe = Peclet number
- = (DNPe) -I, _P- = IlDn'/lp, ................... (25)
Do ¢Do p = pressure
S( C) = source function
t = time
where we have used the integral scale A in place of the grain di-
ameter, dp' The agreement between the simulated and experimen- u = superficial phase velocity
tal results is surprisingly good, considering that no fine-tuning was VDP = Dykstra-Parsons coefficient
done to effect the match. In fact, if we shift the experimental curve x = distance in x direction
to the right by about five, the match becomes excellent. x = mean tr1\vel distance
We conclude from this agreement that (1) core-measured disper- llXD = dimensionless mixing or transition zone
sivities are in reasonable agreement with 01. MA' (2) the mean cor- y = distance in y direction
relation length in the experimental data is about five grain diameters, OI.MA = macroscopic dispersivity
and (3) the transition from a diffusion-dominated to a convection- OI.ME = megascopic dispersivity
dominated dispersivity occurs when A=¢Dolu.
A = mean correlation lengtp
~ = viscosity
Conclusions
P = autocorrelation function of velocity field
I. Megascopic dispersivity in RH permeable media is constant (J x = standard deviation of particle positions
for small correlation distances. The magnitude of this dispersivity
¢ = porosity
depends on the correlation distance, heterogeneity, and aspect ratio,
and weakly on diffusion. Subscripts
2. Dispersivity increases with increasing heterogeneity (Dykstra- D = dimensionless
Parsons coefficient). For smaller VDP (~O.5), OI.ME is constant x = x direction
with respect to time. Macroscopic dispersivity is constant at all y = y direction
times.
3. Megascopic dispersivity increases with the system aspect ratio. Superscripts
At high aspect ratio, OI.ME inf:re1\ses with time. This is true even - = time average (except for x)
for small correlation lengths. Macroscopic dispersivity, however,
is invariant with time at all aspect ratios. AcknowledglT!ents
4. For large aspect ratios, steady-state pressure contours are uni-
Atul Arya acknowledges the financial support provided by the At-
form in the direction nOrmal to the bulk flow, signifying vertical
lantic Richfield Foundation. This investigation was supported by
equilibrium. Thus, for large aspect ratios, an RH medium is analo-
the U.S. DOE under Contract DE-ASI9-82BCI0744 and by the
gous to an equivalent layered system.
Center for Enhanced Oil and Gas Recovery Research at the U. of
5. Diffusion does not influence OI.ME significantly because of the
Texas, Austin.
long length scales. At the macroscopic scale, diffusion promotes
transverse mixing, resulting in time-invariant dispersivities. OI.MA
~eferences
is constant for both correlated and uncorrelated media. For large
diffusion, the effects of heterogeneity are less important. I. Dykstra, H. and Parsons, R.L.: "The Prediction of Oil Recovery by
6. These simulation results also agree with the field observations. Water Flood," Secondary Recovery of Oil in the United States, API
(1950) 160-74.
OI.ME increases with measuring distance in tracer tests. This indi-
2. Lambert, M.E.: "A Statistical Study of Reservoir Heterogeneity," MS
cates long correlation distances, high aspect ratios, and large het- thesis, U. of Texas, Austin (1981).
erogeneity. All of these are naturally occurring phenomena in 3. Smith, P.J. and Brown, C.E.: "The Interpretation of Tracer Test
permeable media. Response from Heterogeneous Reservoirs," paper SPE 13262 presented
7. Dispersion in permeable media can be described by methods at the 1984 SPE Annual Technical Conference and Exhibition, Houston,
used in turbulence theory. Simulation results in this study agree Sept. 16-19.
with the theoretical observations of Taylor24 and Mercado. 20 4. Haldorsen, H.H. and Lake, L.W.: "A New Approach to Shale Manage-
8. Heterogeneous permeable media can be classified as a func- ment in Field-Scale Simulation Models," SPEJ (Aug. 1984) 447-57.
5. Greenkorn, R.A. and !<essler, D.P.: "Dispersion in Heterogeneous
tion of the correlation length, A: if A is much less than L, the medium
Nonuniform, Anisotropic Porous Media," Flow Throu,gh Porous Media.,
is diffusive; if A is roughly equal to L, it is effectively layered; and American Chern. Soc., Washington, DC (1970) 59.
if A is much greater than L, it is strictly layered. 6. Lallemand-Barres, A. and Peauclecerf, P. : "Recherche des relations
entre Ie valeur de la dispersivite macroscopique d'un milieu aquifere,
Nomenclature ses autres caracteristiques et les conditions de mesure," Bull., ~RGM
(1978) 2e Serie, Sec. III, No.4.
b = system dimension normal to bulk flow (width) 7. Pickens, J.F. and Grisak, G.E.: "Scale-Dependent Dispersion in a Strati-
C = concentration fied Granular Aquifer," Water Resour. Res. (1981) 17, No.4,
C v = coefficient of variation of velocity field 1191-1211.
CCO = autocovariance of sample f lags apart 8. Arya, A.: "Dispersion and Reservoir Heterogeneity," PhD dissertation,
U. of Texas, Austin (1986).
dp = particle diameter
9. Perkins, T.K. and Johnston, O.C.: "A Review of Diffusion and Dis-
D = dimensionless diffusion coefficient persion in Porous Media," SPEJ (March 1963) 70-80; Trans., AIME,
Df = fractal dimension 228.

SPE Reservoir Engineering, February 1988 147

 10. Naiki, M.: "Numerical Simulation of Polymer Flooding Including the 21. Zapata, V.I. and Lake, L.W.: "A Theoretical Analysis of Viscous
Effects of Salinity," PhD dissertation, U. of Texas, Austin (1979). Crossflow," paper SPE 10111 presented at the 1981 SPE Annual Tech-
11. Brigham, W .E.: "Mixing Equations in Short Laboratory Cores," SPEJ nical Conference and Exhibition, San Antonio, Oct. 4-7.
(Feb. 1974) 91-99; Trans., AIME, 257. 22. Lake, L.W. and Hirasaki, GJ.: "Taylor's Dispersion in Stratified
12. Brenner, H.: "The Diffusion Model of Longitudinal Mixing in Beds Porous Media," SPEJ (Aug. 1981) 459-68.
of Finite Length, Numerical Values," Chem. Eng. Sci. (1962) 17, 23. Gelhar, Lynn W. et al.: "Stochastic Analysis of Macrodispersion in
229-43. a Stratified Aquifer," Water Resour. Res. (Dec. 1979) 15, No.6,
13. Stehfest, H.: "Numerical Inversion of Laplace Transforms," Commu- 1387-97.
nications of the ACM (Jan. 1970) 13, No.1, 47-49. 24. Taylor, G.I.: "Diffusion by Continuous Movements," Proc., London
14. Schwartz, F.W.: "Macroscopic Dispersion in Porous Media: The Con- Mathematical Soc. (1921) 20, 196-212.
trolling Factors," Water Resour. Res. (1977) 13, No.4, 743-52. 25. Dagan, G.: "Solute Transport in Heterogeneous Porous Formations,"
15. Mantoglou, A. and Wilson, J.L.: "The Turning Bands Method for Simu- J. Fluid Mech. (1984) 145,151-77.
lation of Random Fields Using Line Generation by a Spectral Method," 26. Mandelbrot, B.B. and Van Ness, J.W.: "Fractional Brownian Motions,
Water Resour. Res. (1982) 18, No.5, 1379-94. Fractional Noises and Applications," SIAM Review (1968) 10, 422.
16. Heller, J.P.: "Observations of Mixing and Diffusion in Porous Media," 27. Mandelbrot, B.B.: The Fractal Geometry of Nature, W.H. Freeman
Proc., Second Inti. Symposium on Fundamentals of Transport Phenom- and Co., San Francisco (1982).
ena in Porous Media, IntI. Assn. of Hydraulic Research, Guelph, Canada 28. Mandelbrot, B.B. and Wallis, J.R.: "Some Long-Run Properties of
(1971) 1-26. Geophysical Records," Water Resour. Res. (1969b) 5, 321-40.
17. Journel, A.G. and Huijbregts, CJ.: Mining Geostatistics, Academic
Press Inc., New York City (1978). SI Metric Conversion Factor
18. Byers, E. and Stephens, D.B.: "Statistical and Stochastic Analyses of
Hydraulic Conductivity and Particle-Size in a Fluvial Sand," Soil Sci. ft x 3.048* E-Ol m
Soc. Am. J. (1983) 47, No.6, 1072-81.
19. Bakr, A.A. et at.: "Stochastic Analysis of Spatial Variability in Sub- * Conversion factor is exact. SPERE
surface Flows. 1. Comparison of One and Three-Dimensional Flows,"
Water Resour. Res. (1978) 14, No.2, 263-71.
Original SPE manuscript received for review Sept. 22, 1985. Paper accepted for publica-
20. Mercado, A.: "The Spreading Pattern ofInjected Water in a Permea-
tion Dec. 15, 1986. Revised manuscript received April 20, 1987. Paper (SPE 14364) first
. bility Stratified Aquifer,' ~ Symposium of Haifa: Artificial Recharge and presented at the 1985 SPE Annual Technical Conference and Exhibition held in Las Vegas,
Management of Aquifers, IASH-AISH Publ. (1967) 72, 23-36. Sept. 22-24.

148 SPE Reservoir Engineering, February 1988