PHYSICAL REVIEW E

VOLUME 58, NUMBER 5

NOVEMBER 1998

Spatial correlations and dispersion for uid transport through packed glass beads studied by pulsed eld-gradient NMR
S. Stapf,1,* K. J. Packer,1 R. G. Graham,1 J.-F. Thovert,2 and P. M. Adler3
Department of Chemistry, University of Nottingham, Nottingham NG7 2RD, United Kingdom 2 ` Laboratoire des Phenomenes de Transport dans les Melanges (LPTM), 86960 Futuroscope Cedex, France 3 Institut de Physique du Globe de Paris (IPGP), Tour 24, 4 Place Jussieu, 75252 Paris Cedex 05, France Received 31 March 1998 The two-dimensional displacement joint probability density P (X,Z) for water owing through a bed of glass beads has been measured by means of pulsed eld-gradient nuclear magnetic resonance. The simultaneous particle displacements X and Z perpendicular and parallel to the pressure gradient, respectively, at a given encoding time , are obtained from an experiment employing orthogonal magnetic eld gradients. The resulting probability density distribution is compared to numerical simulations of ow through an equivalent system of randomly deposited monodisperse spheres. The dependence of the centered second moments in X and Z on ow time is discussed for the experimental and simulated data. A crossover from a time scale dominated by Brownian motion toward a behavior determined by the convective ow and velocity uctuations is observed. The mutual dependence between displacements perpendicular and parallel to the ow direction is revealed in the evolution of a correlation coefcient X 2 ,Z . This coefcient is found to increase for short times and to decrease for larger displacements, with a maximum at an average displacement corresponding to the bead radius. As a means of displaying the cause of these correlations, a correlation probability density C (X,Z) P (X,Z) P (X) P (Z) is suggested, where P (X) and P (Z) are the marginals of P (X,Z). A plot of this matrix renders zero in the absence of correlations, but produces a characteristic pattern of positive and negative regions when displacements in X are correlated with those in Z. The time evolution of this pattern is discussed and compared to the shape of model propagators obtained from an analytical function and a numerical simulation for a simplied capillary array, respectively. S1063-651X 98 05611-6 PACS number s : 47.55.Mh, 81.05.Rm, 47.11. j, 83.70.Gp
1

I. INTRODUCTION

Transport and dispersion of uid phases and solutes within porous solid structures is of importance in a wide range of areas such as oil reservoir appraisal and management, aquifer behavior, distillation and ltration processes, heterogeneous catalyst bed design and performance, pollutant dispersal and recovery in the environment, etc. The eld has an extensive literature, and a wide range of experiments have been performed, many of which have been summarized in three major references 13 , with historically important compilations and data being also found in Refs. 4,5 . Few theoretical results are available to describe these uid transport processes, except for the case of Poiseuille ow 6,7 and for dilute suspensions 8 . Numerical results for long time behavior have been obtained for two- 9,10 and threedimensional 11 structures. Some consideration has been given to the early time dispersion behavior of solutes, the so-called non-Fickian or nonlocal regime, both experimentally 12,13 and theoretically 14,15 . In the latter approach, the convection-diffusion equation was solved symbolically in terms of a Green function P, which gives the probability of nding a tracer at position r at time t, given that it was

*Author to whom correspondence should be addressed. FAX:

known to be at r0 at t 0 . This Green function was used to express the average mass ux as a function of the initial concentration distribution, with the resulting expression being rearranged and a time-dependent dispersion tensor introduced, which generalizes the classical local dispersion. So far, this has been worked out only for uid transport through a dilute suspension of spheres, although it was claimed that it could be extended to larger solid phase volume fractions. Experimental data from closely packed beds of spheres 13 have shown only modest and qualitative agreement with this theory. Of the many experimental approaches to the characterization of uid transport in porous solids, nuclear magnetic resonance NMR has a number of signicant advantages. Principal among these is the fact that NMR studies the uid directly and is able to investigate optically opaque systems, which constitute the majority of those of interest. A particularly powerful NMR method is that based on the use of magnetic eld-gradient pulses to determine the statistics of nuclear spin and, hence, molecular displacements 16 . Such pulsed gradient spin echo PGSE NMR methods are based on the propagator formalism 17 which directly gives P (R), the probability distribution of displacements R r( ) r(0) in time . This method has been used experimentally to characterize the diffusive-convective transport through a number of model porous solids 1825 . The availability of increased computational power has made possible simulations with sufcient spatial and temporal resolu6206 1998 The American Physical Society

PRE 58

SPATIAL CORRELATIONS AND DISPERSION FOR . . .

6207

FIG. 1. Pulse sequence used for measuring two-dimensional average propagators. rf hard pulses are given by black rectangles, and the slice selective pulse is indicated by its sinc shape. Gradient pulses in dark gray are encoding gradients, and crusher gradients to remove residual phase coherences are drawn in light gray. The encoding gradients are applied simultaneously in orthogonal directions.

tion to compare directly with experiment 2325 , making use of algorithms which generate realizations of statistical porous systems 11,26 . All measurements reported to date have determined P (R) for R either parallel or transverse to the pressure gradient driving the net ow. However, to our knowledge, no attempt has yet been made to correlate displacements parallel (Z) and perpendicular (X) to the pressure gradient quantitatively with each other for ow in a porous system, either by experimental or computational methods. Preliminary results on ow through a porous sandstone were presented in Ref. 27 . In this paper we determine, by PGSE NMR, the two-dimensional joint probability density P (X,Z) for water ow in a bed of randomly deposited glass beads. The time evolution of this propagator, as well as that of the moments and correlation coefcients connected to X and Z, are described and compared with extensive computer simulations.
II. THEORY

proportional to its position at time 0 and to the area of the gradient,

,ri 0

ri 0

B0

gri 0 ,

where g and are the strength and duration of the applied gradient, respectively, is the gyromagnetic ratio, and B 0 denotes the static magnetic eld. (ri ) is the Larmor frequency at the position ri . After an evolution time , the refocusing gradient results in a negative phase shift which leaves the resultant shift
i

g ri

ri 0

gRi

where Ri ( ) ri ( ) ri (0) indicates the displacement for particle i during the observation time . The total signal amplitude is obtained by summation over all spins, equivalent to the integral S q P R exp i2 qR dR, 3

To obtain the two-dimensional propagator, we have modied the alternating pulsed eld-gradient stimulated echo APGSTE version 28 of the PGSTE sequence 29 . In addition to the original sequence, a second set of gradient pulses has been added; gradients in both orthogonal directions are switched simultaneously see Fig. 1 . The splitting of the defocussing and refocussing gradients by insertion of a rf pulse of ip angle minimizes the effect of molecular displacements through the internal magnetic eld gradients, which arise from the differences of magnetic susceptibility existing between the porous solid matrix and the saturating uid 30 . As described in Ref. 28 , offsetting the phase of the second /2 pulse by 90 allows discrimination between positive and negative displacements. Each nuclear spin i experiences a phase shift i that is

where q (2 ) 1 g, and the average propagator P (R) P(r0 ) P(r, ;r0 )dr0 . P(r0 ) is the probability density for starting positions, while P(r, ;r0 ) is the conditional probability for displacements from r0 to r in time , equivalent to the Green function mentioned above 14,15 . The average propagator P (R) can therefore be obtained directly by Fourier transformation of S (q) with respect to q. Under the circumstances of this study, q consists of two orthogonal components q z and q x ; thus q kq z iq x , 4

where k and i are unit vectors along z and x, respectively. z is parallel to the pressure gradient driving the uid ow, and

6208

STAPF, PACKER, GRAHAM, THOVERT, AND ADLER

PRE 58

x is perpendicular to it. The two-dimensional average propagator P (X,Z), for displacements X and Z, is then calculated by successive fast-Fourier transform in both directions. It must be mentioned that Eq. 3 is only applicable if the duration of the gradient pulse is negligible compared to the experimental time ( ). The boundary conditions imposed experimentally see below ensure that there is no net radial ow and the system is therefore axially symmetric. The propagator for displacements in the x-direction, P (X), is therefore identical for all directions perpendicular to the ow axis.
III. EXPERIMENT

The NMR measurements were carried out using a GE CSI spectrometer operating for proton resonance at 85 MHz, the eld being provided by an Oxford Instruments 85/310 horizontal bore magnet equipped with room temperature shims and S-150 Accustar actively shielded gradient coils providing gradients of up to 0.2 T m 1. Phase cycling of the rf pulses 28 was used to minimize the effects of background gradients and dc offsets. A sample of glass beads of 600 50- m diameter Jencons Scientic Ltd., Leighton Buzzard, UK, Cat. No. H102/1/126 was prepared by lling a glass tube of 27.2-mm inner diameter and 120-mm length with water and adding the wet glass beads and water suspension slowly, allowing the beads to sediment. Air bubbles were removed by stirring during the lling process. The tube was tted at both ends with sintered glass disks to ensure an even distribution of streamlines over the whole cross-section. The sample was connected, via a 2-m narrow-bore pipe, to a precision pump Pharmacia P50 which was operated at constant volume ow rates of 4.8, 14.0, and 42.0 ml/min, respectively. The relaxation time of the owing water was reduced to T 1 600 ms by adding copper II sulfate, allowing a pulse sequence repetition time of 3.3 s. Experiments were performed using the APGSTE pulse sequence given in Fig. 1. The signal was acquired using a slice-selective soft pulse with an effective axial slice thickness of 60 mm located symmetrically at the center of the sample. This was used in order to avoid edge effects which could arise from the inow and outow of water at the ends of the sample. Two-dimensional data sets were obtained by stepwise variation of the strength of the pulsed gradients simultaneously applied parallel and perpendicular to the ow direction. Data were acquired using n X 2 and n Z 2 evenly spaced positive and two negative gradient values, respectively, where n X and n Z denote the total dimensionality in X and Z directions, and were chosen as powers of 2 typically n X 16, n Z 32 or 64 . Measurements at the two negative gradient values served to determine the zero-q phase shift which enabled the reconstruction of the full data set of dimension (2n X 4) (2n Z 4). This matrix was then Fourier transformed numerically after zero lling to twice the original size, and subsequently phase corrected. The marginals P (X) and P (Z) of the resulting propagator P (X,Z), dened as P X P X,Z dZ, P Z P X,Z dX, 5

FIG. 2. Example of a reconstructed packed bed of monodisperse spheres.

were compared to one-dimensional propagators obtained under identical conditions but with a higher resolution of 64 or 128 points and with symmetrical gradient steps covering the range q max ,...,qmax . The propagators matched satisfactorily. All calculations were performed using interactive data language IDL 31 .
IV. SIMULATIONS

To generate an adequate representation of the real porous medium, random sphere packing has been simulated by successive deposition of grains in a gravitational eld. The Nth grain is introduced at a random location above the bed of N 1 grains already deposited, and is allowed to fall until it reaches a local minimum of its potential energy. A more detailed description of the deposition process is found in Refs. 26 and 24 . A random packing of monodisperse spheres was generated by this algorithm, incorporating periodic boundary conditions along the two horizontal axes see Fig. 2 . The sample is made of N 3 elementary cubes of side C 1 a 60 m, corresponding to 10 of the bead diameter, with N C 64. The porosity of the void space was chosen to be 0.44. The bed permeability K was calculated by solving the Stokes equations, and was found to be K 7140 Darcy. In the next step, the velocity eld was generated by solving the Stokes equations p 2 v, v 0, 6

where v, p, and are the velocity, pressure, and viscosity of the uid, respectively, and v 0 on the surface of the wetted solid. The symmetric permeability tensor K only depends on

PRE 58

SPATIAL CORRELATIONS AND DISPERSION FOR . . .

6209

the geometry of the system, and describes the relation between the macroscopic pressure gradient p and the seepage velocity : v v 1 K p. 7

The numerical method used to solve these equations is outlined in Ref. 32 . It assumes low Reynolds numbers which are guaranteed by the experimental conditions employed in this study. The determination of the average propagator is performed by inserting a large number of particles uniformly distributed within the pore space. For each elementary time step, the particles position is calculated by adding convective and random diffusive displacements, where the geometrical restrictions of the solid matrix are taken into account. The relative weight of these contributions to the displacement is expressed by the Peclet number Pe
v *L , D

where v * is the interstitial velocity and L is a characteristic length, taken as being equal to the sphere diameter. The random component is adjusted for a given Peclet number to be as large as possible in order to speed up the statistical convergence, provided that the total elementary jump length is kept smaller than a/2 24 . The precision of these calculations was carefully studied in Ref. 11 ; it was concluded that the calculations were reliable for Peclet numbers smaller than 1000. Following this procedure, typically 2.5 105 particles were distributed randomly in the pore space of the lattice, and were allowed to undergo ow and Brownian motion. The self-diffusion coefcient was chosen as 2.1 10 9 m2 s 1, the sphere diameter as 6 10 4 m, and the interstitial velocity as 2.85 10 4 , 8.33 10 4 , and 2.5 10 3 m s 1, respectively. The latter values correspond to the ow rates used during the experiments, i.e., 4.8, 14.0, and 42.0 ml/min, respectively. Additional simulations were run with smaller numbers of particles for time scales well above and below the experimentally accessible range, and for which the self-diffusion coefcient was also varied over several orders of magnitude. This simulation technique for obtaining the displacement distributions of the particles is, in a sense, strictly equivalent to the formalism using the Green function which was developed by Koch and Brady 14,15 , since it amounts to the solution of the convection-diffusion equation. For this reason, it was not found necessary to rederive it within the nonlocal formalism of Refs. 14,15 .
V. MODELING

laries of length l and uniform radius, isotropically distributed in orientation. 500 particles are assigned to the start of each capillary at t 0. The velocity vectors are oriented in the direction of the capillary, and have magnitudes taken at random from the usual uniform distribution appropriate to laminar ow within a circular pipe. To represent the axial nature of the system, we make v max V0 cos , where is the angle between the capillary and the z axis. Coordinates of each particle are then calculated at various times . The important extension of the model from that in Ref. 18 is that particles which reach the end of their starting capillaries within transfer to the start of another capillary in which they continue their motion. The orientation of this second capillary lies on a cone of half-angle emanating from the end of the starting capillary. is constant for each calculation, but is varied in order to give insight into the way in which locally discrete changes in ow direction inuence the propagators P (X,Z). The nal results from all 5 106 trials for each are collected in a 64 64 64 three-dimensional histogram, and stored on disc for further processing. Although the use of such a model may be criticized as being based on a physically unrealistic representation of a real pore space, such models and more sophisticated versions of them are widely used for investigating uid transport in porous solids 33 .
VI. RESULTS AND DISCUSSION A. Experiments and simulations on spherical beads

To investigate the inuence of characteristic system dimensions, both ow rate and encoding time were varied in the experiments over a range of average displacements covering two orders of magnitude from much less than the bead size of 600 m to about 2 mm. The experiments were restricted by the required condition / 1; the smallest ratio used was larger than 8. On the other hand, encoding times exceeding the longitudinal relaxation time T 1 by far were not considered feasible due to the loss of signal intensity. This range, however, could be expanded in the simulations where similar restrictions do not apply. In addition, the ow rates had to be chosen to result in Peclet and Reynolds numbers Pe and Re to comply with the conditions of the simulations. For the three ow rates used in the experiment 4.8, 14.0, and 42.0 ml/min, respectively , one obtains Pe 88, 245, and 720, where Eq. 8 has been used with L 600 m. Some simulations were run with higher Peclet numbers but only for short times, where it has been found that the method used is still valid 24 . The Reynolds number is dened as Re
v *L

In order to gain a more direct physical insight into the relationship between the features observed in the experimentally determined joint propagators P (X,Z) and the underlying pore space characteristics, we have extended the simple model for ow in a porous solid which we introduced earlier 18 . The model used here consists of 104 cylindrical capil-

with being the kinematic viscosity of the liquid. Re indicates the ratio of inertial and viscous forces, and should not be much larger than unity for optimal simulation results; in this case, it took the three values 0.15, 0.5, and 1.4, respectively. For such low Reynolds numbers, the deviations of the real velocity eld from the solution of the Stokes equations are expected to be very small. Two-dimensional propagators joint probability densities for displacements X and Z) for water owing in the column

6210

STAPF, PACKER, GRAHAM, THOVERT, AND ADLER

PRE 58

FIG. 3. Two-dimensional average propagator P (X,Z) for experimental data at a ow rate of 42.0 ml/min. All propagators are normalized to P (X,Z)dX dZ 1. Contour lines are drawn from approximately 0.05 of the peak intensity in linear spacing; numbers indicate probability densities in 104 m 2. Evolution times are as indicated.

of packed glass beads are shown in Fig. 3 for the 42.0-ml/ min ow rate, and for encoding times between 64 and 761 ms. The contour lines indicate regions of equal probability density. Numbers at each line are given in 104 m 2 and the propagator is normalized so that P (X,Z)dX dZ 1. The propagator is symmetric in X, which is a consequence of the experimental constraint that no net ow occurs perpendicular to the pressure gradient. In the Z direction, the fraction of particles experiencing negative displacements is very small, while a pronounced peak near zero displacement is found for short times. It can be seen immediately that the lines of equal probability density spread in both directions with increasing time. At the shortest time 64 ms, when the mean displacement Z in the ow direction is 160 m, considerably less than the bead size, and when the root-mean-square displacement due to diffusion, Z 2 ( ) , is only 16 m, the total spread 1 at the outermost line corresponding to about 20 of the peak probability density is almost the same in X and Z. Thus, for an average displacement less than the bead size, the particles spread to a similar degree in X and Z. The same is observed from experiments with smaller ow rates and hence even smaller average displacements. For longer times and for larger ow rates, however, the shape of the propagator becomes more elongated, and large particle displacements occur preferentially along Z. A closer look at the shape of individual contour lines reveals that the ratio of their main axes, Z/X, is changing, this ratio being larger for contours of higher probability density. This feature suggests that near the higher levels of probability density, which are associated

with small X displacements, particles possess, on average, much larger displacements in the axial direction than perpendicular to it. In contrast, for those particles which have traveled large distances, the difference between X and Z is considerably less pronounced. This deviation from a symmetric shape of the two-dimensional propagator already indicates a connection between X and Z displacements. Another feature apparent in Fig. 3 is the position of the peak for each ow time . While it remains near zero for shorter times, it becomes shifted toward larger displacements only for times of 340 ms and larger. A similar behavior has been found previously with one-dimensional measurements 18,24 : A narrow Gaussian peak around Z 0 develops a shoulder for increasing times and nally disappears for long times; the shoulder, on the other hand, gives rise to another peak that eventually becomes the center of a Gaussian distribution, and is given by the average displacement Z v * . Observing the one-dimensional propagator P (Z) alone leads to a simplied interpretation of two components, one quasistatic and the other moving with v v * . This can be seen by comparing the two-dimensional propagator with its marginals see Fig. 4 : in this example for 42-ml/min ow rate and 231 ms, the peak for P (X,Z) is still at Z 200 m, while its projection P (Z) P (X,Z)dX shows a maximum probability density near 600 m. As the marginal represents an average over all X displacements, part of the information about the real properties of the displacement probabilities in two or three dimensions is lost, which can lead to misinterpretations.

PRE 58

SPATIAL CORRELATIONS AND DISPERSION FOR . . .

6211

FIG. 4. Two-dimensional average propagator P (X,Z) for experimental data at a ow rate of 42.0 ml/min, 231 ms. The marginals P (X) and P (Z) are drawn along the X and Z axis, respectively.

In Fig. 5, simulated propagators at equivalent times are shown for comparison. A similarity is observed for the general shape of the contour lines and their axis ratios. However, a major difference can be found with respect to the peak of maximum probability density. It is much more persistent in the simulation than in the experiment, being prominent even at the longest evolution time of 750 ms. In the simulations, it eventually disappears for 5 s. Experiments and simula-

tions generally coincide much better for small average displacements, as long as the marginal P (Z) can be described by a peak near Z 0 and an asymmetrically decaying tail at larger displacements 18,24,25 . The difference between experiment and simulations has been discussed previously, as possibly arising as a consequence of the surface relaxivity that tends to remove a fraction of the molecules near the walls in the NMR experiment 25 . Due to the nonslip condition near the wall, this affects mainly spins with small velocities, and might lead to a faster decay of the peak near Z 0. However, allowing a loss of particles at the wall did not inuence the result in earlier simulations 24 . The inuence of large Peclet and Reynolds numbers was also mentioned in Ref. 24 . As in this previous investigation, edge effects can be ruled out as the inner diameter of the tube is equal to 45 bead diameters. It must be pointed out that the deviations only affect the peak at small displacements, and therefore a fraction of particles that show small net displacements apart from their Brownian motion; this fraction disappears much later in the simulation. However, the behavior for short times is well represented by the simulations, as well as the probability densities for large displacements. The numerical simulations coincide reasonably well with the experimental results once these limitations are considered.
B. Time and displacement dependence of moments of P X,Z

In Sec. VI A, we described the evolution of the propagator P (X,Z) as a function of time in a qualitative way by

FIG. 5. Two-dimensional average propagator P (X,Z) for simulated data 250 000 particles at a ow rate of 42.0 ml/min. All propagators are normalized to P (X,Z)dX dZ 1. Contour lines are drawn from approximately 0.05 of the peak intensity in linear spacing. Evolution times are as indicated.

6212

STAPF, PACKER, GRAHAM, THOVERT, AND ADLER

PRE 58

FIG. 6. Average displacement in the ow direction, Z , as a function of time for experimental open symbols and simulated data solid symbols . Solid lines indicate a linear t to the experimental data. Numbers denote ow rates in ml/min.

discussing its shape and features such as global and local maxima of probability density. For a quantitative description we will rst focus our attention on the development of moments of P (X,Z) with increasing encoding time and their dependence on ow rate and Peclet number. The general denition of moments of nth order is given as Xn P X,Z X n dX dZ, Zn P X,Z Z n dX dZ. 10 It can be expected that the evolution of moments is entirely determined by the spatial structure of the porous system and the consequent properties of the velocity eld. The moments are only indirectly affected by the time variable inasmuch as it separates regimes where diffusion and ow, respectively, are dominant for the resulting mean-squared displacements or second moments. For the rst moment, however, a simple relationship is found: X 0, Q A 11

leads to nonzero rst moments X . This effect, however, can be controlled, and is avoided when centered moments are analyzed. In Fig. 6, the average displacements Z obtained from experimental data and from the simulations are compared for all three ow rates. Both sets coincide satisfactorily. The slope of the experimental data yields an interstitial velocity of (3.1 0.1) 10 4 m/s for the ow rate 4.8 ml/min, (8.6 0.2) 10 4 m/s for 14.0 ml/min, and (2.53 0.05) 10 3 m/s for 42.0 ml/min, respectively. Given the dimensions of the system, this leads to an effective porosity of the sample of ( 45 2)%, in good agreement with the value assumed for the simulated random sphere packing. From the above result, it can be concluded that the rst moment Z in fact scales linearly with time ; in all further discussions, the time variable can therefore be replaced by the average displacement in the ow direction. The second moments of displacements parallel and perpendicular to the pressure gradient, Z 2 and X 2 , contain contributions from both Brownian motion and convection. The second moment from self-diffusion alone is isotropic and given by the relation X2 Z2 2D . 13

For an innite isotropic medium, D( ) is a timeindependent constant and equal to the self-diffusion coefcient D 0 . Diffusion within a restricted geometry shows a more complicated pattern. While for displacements much smaller than the wall separation, DD 0 is found, displacements much larger than both the average pore size and the correlation length of the pore space lead to a constant selfdiffusion coefcient reduced by a certain factor which depends on the porosity and the tortuosity of the system 35 . In the intermediate range, D( ) can be expressed 36 by the short-time expansion D D0 1 4 9 S V D0 O , 14

v*

12

where Q is the ow rate and A the cross-section of the sample. The average interstitial velocity v * is also known as the Dupuit-Forcheimer velocity 34 , and is a well-dened quantity as the setup of the experiment guarantees a constant ow rate Q and thus a constant v * . The average displacement in the ow direction, Z , must then be proportional to time as the contribution due to self-diffusion remains zero for all times. As no net ow occurs in X, and Brownian motion leads to an isotropic spreading in all directions, X must be zero. The assumption of an isotropic medium seems justied due to the large number of beads ( 1.7 105 ) although the local porespace favors anisotropic spreading. The smaller number of unit cells in the simulation, containing only 150 beads, can give rise to a certain asymmetry which

where S and V are the surface area and the volume of the porous system. However, for the glass bead system investigated in this study, a signicant decrease of D( ) can only be expected for times in the order of 1 s, where the inuence of ow is already dominating. Taking ow into account, in the limit of innite times, the propagator is expected to become a Gaussian centered at X 0 in the X direction and a Gaussian shifted by the average displacement Z 0 v * in the Z direction. It has been shown that the shape of the propagator becomes roughly Gaussian for average displacements much larger than the bead size in systems similar to the one investigated in this work 19,24 . However, as we are interested in the intermediate regime where this limit is not yet reached, we investigate the time dependence of the second moments of X and Z. The spreading of the propagator in X and Z is best monitored by looking at the centered second moments (X X ) 2 and (Z Z ) 2 . While the former is identical to 2 X as X 0, the latter describes the spreading of the probability density distribution around the center at Z v * . A plot of the centered second moments is shown in

PRE 58

SPATIAL CORRELATIONS AND DISPERSION FOR . . .

6213

FIG. 7. Dependence of the centered second moments (X X ) 2 and (Z Z ) 2 on the average displacement Z . a Centered second moments; experimental results are denoted by larger, open symbols. b Derivative dlog(X X ) 2 ,(Z Z ) 2 /d(log Z ) from simulated data indicating the exponent in the expression (X X ) 2 ,(Z Z ) 2 Z .

Fig. 7 a . Moments in both directions are found to follow a sigmoidal pattern with a larger slope for intermediate displacements and a shallower dependence on Z for both small and large average displacements. Assuming a powerlaw relationship illustrated for Z by

15

see Fig. 7 b reveals that 1 a plot of the exponent describes the region where self-diffusion dominates the spreading of the propagator function at short times, while

6214

STAPF, PACKER, GRAHAM, THOVERT, AND ADLER

PRE 58

1 is again reached asymptotically for very large displacements exceeding 10 mm, reecting the limiting pseudorandom walk character of the displacement distributions expected at long times. In between, a maximum in is found at a displacement which depends on the Peclet number, and is shifted toward larger average displacements for smaller Pe. is generally larger for (Z Z ) 2 than for (X X ) 2 . For the highest ow rate used in this study, 2 is almost reached, which can be expected to be the maximum value possible; it is considerably less for the smaller ow rates. Note the minimum 1 in (X X ) 2 occurring at Z 1 mm. This minimum might be related to the particular shape of the two-point correlation function of the matrix which possesses a negative region anticorrelation between 0.9 and 2.0 bead radii for a random packing of monosized spherical particles 26 . Indeed, Z at the position of the minimum corresponds to an average transverse displacement (X X ) 2 of about 350 m or 1.2 bead radii see Fig. 7 a . While this is not a sufcient model for the system studied here, it is of interest to compare the behavior of (Z Z ) 2 with that predicted for Taylor dispersion for ow between parallel plates 37 , i.e., Z Z
2

2D 0

v *2a 2 2 D0

v *2a 4
4

D2 0

1 e

2D

2 0 /a

16 The limiting cases are found to be 2D 0 Z Z

v *2

a 2 /D 0 a 2 /D 0 . 17

FIG. 8. Ratio of dispersion coefcients to the self-diffusion coefcient D 0 for experimental and simulated data. a Parallel to the ow direction, D * ( Z ). b Perpendicular to the ow direction, D * ( Z ).

2D 0

v *2a 2 2 D0

O 1 ,

1 to 2 is observed for short disA crossover from placements, equivalent to short times where a 2 /D 0 whereas a decrease of back toward 1 is predicted for the condition c a 2 /D 0 . The observed behavior is qualitatively similar to the classical Taylor dispersion, showing a crossover from 1 to 2 at a time scale when the ow effect on the mean-squared displacement becomes dominant over Brownian motion. At larger Z , when the spreading of the propagator almost behaves as in the case of pure directed ow, hence 2, this similarity is lost as the spreading is now dominated by the pseudorandom walk caused by ow around the obstructions of the porous solid. A comparison of the spreading of the particle probability density function due to ow and as a consequence of Brownian motion is achieved by determining the ratio of the components of the dispersion tensor D* and the self-diffusion coefcient D 0 . The dispersion tensor D* can be regarded as a function of time and therefore average displacement: D* D * It is dened by D* d R d R
2

where (R R ) 2 is the centered second moment in three dimensions. In an isotropic porous medium with an interstitial velocity v* parallel to the z axis, D* can be decomposed into parallel and transverse components 11 D* 0 0 0 D* 0 0 0 . D*

D*

20

D* Z .

18

19

The dispersion tensor remains time dependent as long as 1. In Fig. 8, the components of the dispersion tensor perpendicular and parallel to the ow direction are compared for the three different ow rates used in this study. D * is found to become constant for average displacements Z 300 400 m, a value that seems to be roughly independent of the Peclet number. From this point on, the displacements in X are spreading in a way that is characteristic for a Gaussian propagator. In their treatment of dispersion, Koch and Brady 8 suggested that D * was determined entirely by mechanical contributions and estimated that it would reach its limiting value in a time c a/( v * 1/2), where is the solids volume fraction. This time is equivalent to Z a/ 1/2 which for 0.55 gives Z 800 m for the grain diameter of a 600 m. This is within a factor of 2 of the values represented in Fig. 8 b . D * approaches its asymptotic limit much more slowly than D * , and does not become constant within even the

PRE 58

SPATIAL CORRELATIONS AND DISPERSION FOR . . .

6215

range of the simulations. For the highest ow rate of 42 ml/min, D * is still changing for average displacements Z exceeding 100 bead radii. The limiting values of D * can, however, be estimated reasonably well; they fall in the same range as results presented in the literature see Ref. 20 for a compilation . The characteristic time at which the longitudinal dispersion coefcient becomes roughly constant can be estimated from the simulated data to about 20 s; this is in good agreement with relations given in the literature, where the quantity D 0 /L 2 is of order unity 38 . Koch and Brady 8 derived expressions for the time required for the nonmechanical contributions to longitudinal dispersion to reach their asymptotic limits. The mechanism of relevance to this study is the so-called boundary layer dispersion arising from the nonslip condition of the velocity eld at the walls of the solid matrix. The time for this process to reach its asymptotic limit is of order a Pe1/3/ v * which, when converted into values of Z appropriate to the three velocities used here, give Z in the range (3 5) 103 m. This can be seen to underestimate the values suggested by the simulation results in Fig. 8 a . The latter are supported by the experimentally derived data in that, at the limit of the experiments, corresponding to Z 103 m, there is no sign of approach to asymptotic behavior. In Ref. 20 , in which the dispersion coefcient was determined by PGSE NMR from the low-q data, a much faster approach toward D * ( ) const was found for a similar system of spherical beads. An identical analysis of the measurements presented in this work led to the same values of D * ( ) as those obtained from the full propagator, as described above, and we are unable to offer any explanation for this difference. The dependence of the dispersion coefcients on Peclet number has been written 39 in the dimensionless powerlaw form D * /D 0 Pe . 21

FIG. 9. The residence time distribution c (Z ) as a function of for ow through a 600- m glass bead pack at 14 ml/min. The denitions of Z and are dimensionless as given in Ref. 14 : 2 1 2 Z Zk 1/2 and Z k 1/2, where k 1/2 9 a is the screening length. The corrected Peclet number according to Ref. 14 is Pe v * k 1/2/D 150.

For Pe much larger than unity, 2 is found if diffusion is the predominant mixing effect between different velocities, as is the case for dispersion in laminar ow in a capillary tube. If velocity variations within the network are the main cause for dispersion, 1 is observed. For a system containing regions where the uid velocity is large and others where it is small, an intermediate behavior is expected 11 . In Ref. 26 , ow through systems containing either spherical or ellipsoidal particles with a wide range of aspect ratios was simulated numerically. An average relationship of D * /D 0 0.26 Pe1.29 and D * /D 0 0.27 Pe0.72 was found for Pe 10, which coincides well with the values in this study. In particular, the proportionality D * /D * Pe1/2 40 is found to a good approximation see Figs. 8 a and 8 b .
C. Residence time distributions RTDs from PGSE NMR

elsewhere 1625 . As with c (Z), which allow the calculation of tracer RTDs, c Z ( ), so the P (Z) give access to the equivalent NMR determined RTDs. The formal equivalence requires the assumption that the distributed labeling of molecular positions by NMR is the same as the innitely thin, innite area transverse labeling assumed in tracer measurements. All that is required is that, as the NMR method measures displacements, the distribution of initial labeling is representative of all locations in the uid phase. Based on their nonlocal formulation of the transport problem in porous solids, Koch and Brady calculated c (Z) and c Z ( ) for a Peclet number of 100 and a solids fraction of 0.5 which are similar to those used in this study. It is clear from previously published data, e.g., Refs. 1825 that the calculated c (Z) shown in Fig. 2 of Ref. 14 show little agreement with those observed experimentally. We have used our data to calculate RTDs for the 600- m glass sphere bed using the same denitions of variables as used by Koch and Brady 14 . Figure 9 shows these NMR-derived RTDs for the ow rate corresponding to Pe 245 (Pe 150 according to the screeninglength corrected denition of Ref. 14 and 0.55. It can be seen by comparison with Fig. 3 in Ref. 14 that there is little similarity, even qualitatively, with the predictions of the theory.
D. Correlation between displacements in X and Z

The displacement distributions, such as P (Z), derived from PGSE NMR, are formally equivalent to the tracer distributions which are used to characterize uid transport through porous solids 14,15 . In particular, c (Z) displayed in Fig. 2 of Ref. 14 correspond exactly to the P (Z) determined by the NMR experiments described in this paper and

The principal aim of this paper is to show that our NMR experiment gives direct information on correlations between axial and transverse components of the uid transport. The centered second moments are generally found to be larger along the pressure gradient than perpendicular to it. In fact, they show a mutual dependence for long times when the limit of Gaussian propagators is reached. However, as will be pointed out later, the moments (X X ) 2 and (Z Z ) 2 are not correlated for in the strict mathematical sense, despite this proportionality. The centered second moments represent averaged quantities for the whole system, and do not sufciently describe the actual connection between displacements in both directions.

6216

STAPF, PACKER, GRAHAM, THOVERT, AND ADLER

PRE 58

From the simulation it is clear that for displacements much larger than the diffusion range, X 2 follows a powerlaw dependence of Z. The slope in the log-log plot is found to be 0.74 0.02 both in the experimental and the simulated data for Z 500 m. While for small Z, isotropic selfdiffusion is dominating and X and Z are uncorrelated, there remains a positive correlation between both quantities throughout the full range of Z. The quantication of correlations for the entire range of displacements is necessarily connected with the reduction of specic properties onto a single characteristic gure. The proper mathematical denition of the correlation coefcient which relates two parameters A and B, A,B , is given as follows: cov A,B
A,B

Var A

Var B

22

where covariance and variance are dened as cov A,B Var A AB A2 A B , A 2, 23 24

and where A , for example, is given by A

P A,B A dA dB.

25

Both experimental and simulation results are subject to noise that makes the calculation of higher moments increasingly inaccurate. One is therefore interested in the simplest possible correlation relations, such as
X ,Z , X ,Z 2 , X 2 ,Z , X 2 ,Z 2 .

In particular, one would like to know if, for a given encoding time, a spin-bearing particle that has traveled a larger-thanaverage distance in the ow direction is also expected to have a larger-than-average displacement in the perpendicular direction. In other words: are large displacements in Z correlated with large displacements in X or not? A way to gain insight into this question is to compute the second moments X 2 as a function of Z. One thereby chooses a subset of particles that have reached, by different pathways, a displacement between Z and Z dZ after a time interval . The inverse approach of plotting (Z Z ) 2 vs X is also justied, but less illustrative. Figure 10 a shows a superposition of the experimental data for 42 ml/min, as compared to the simulations in Fig. 10 b . In Fig. 10 a , it can be seen that X 2 which is identical to the centered second moment is only weakly dependent on Z for small Z as this region is dominated by selfdiffusion which is isotropic. The slope then rises monotonically toward larger Z with all values of achieving an approximately common slope in the log-log plot. The fact that the curves for different times are nearer to each other in the simulated data as compared to the experimental ones is a consequence of the sharper and more persistent peak around zero already mentioned. It has the effect that static and owing particles are more strongly separated than in the experimental case. Thus the diffusiondominated region where X 2 scales with time and therefore with Z ) is restricted to a smaller X and Z.

It was found that the general features of these quantities and their evolution with time and ow rate remain the same while only the absolute values are changed slightly. Figure 11 a compares the correlation coefcients X 2 ,Z for the experimentally obtained propagators. Correlations are generally larger for high ow rates, i.e., larger Peclet numbers. Furthermore, an increase of the correlation coefcient with Z is found for 4.8 ml/min, a decrease for 42.0 ml/min, and an intermediate behavior showing a maximum for 14.0 ml/min. Despite an estimated error of up to 10% in the absolute values, the maximum seems to appear for Z between 100 and 300 m for 14.0 ml/min. The interpretation becomes clearer by considering the X 2 ,Z derived from the simulated propagators. These are shown in Fig. 11 b , and cover a wider range of displacements than the experimental data. The dependence on ow rate is obvious. A maximum is clearly seen at displacements near to the bead radius. Additional simulations for 4.8 ml/ min have been run with self-diffusion coefcients of D 0 /100 and D 0 /10000 in order to reveal the effects of ow alone, and these are included in Fig. 11 b . From the behavior of the correlation coefcient, the following interpretation can be drawn. At very short times, particle displacements are dominated by Brownian motion on scales much smaller than the bead size. This motion is random and isotropic, and generates no correlation between dis-

PRE 58

SPATIAL CORRELATIONS AND DISPERSION FOR . . .

6217

FIG. 11. a Correlation coefcients X 2 ,Z for experimentally obtained two-dimensional propagators. Errors bars are estimated to 10% from the scattering of the computed values. b Correlation coefcients X 2 ,Z for two-dimensional propagators obtained by simulation. Open symbols indicate simulations at reduced selfdiffusion coefcients of D 0 /100 and D 0 /10 000, respectively.

placements in X and Z. For larger times, the motion is increasingly inuenced by the velocity eld in the pore space, with its preferred direction along the pressure gradient. The velocity eld represents the way uid particles probe the spatial structure of the pore space under the constraints of the boundary conditions. The local velocity distributions in a pack of glass spheres were discussed in Ref. 41 ; regions of large and small local velocities could be identied which indicate the position of the particles with respect to the main owlines. It is clear that larger Peclet numbers allow fewer particles to change between these pools in a given time, resulting in larger correlation coefcients at any given displacement. For very small self-diffusion coefcients, motions become dominated by the ow eld, as is clear from Fig. 11 b . At a characteristic length scale, uid molecules will have to change their ow direction or exchange between the pools of large and small velocity; the persistence length of a streamline must be connected to the typical structural size of the system, and for very large displacements it becomes increasingly unlikely for a particle to ow in its initial direction. One would therefore expect a maximum of the correlation coefcient that is intimately related to the characteristic length scale, which is given by the bead size in this system. Note that in the complete absence of Brownian motion, the correlation coefcient will start off as a constant, the value of

which is solely determined by the spatial structure of the velocity eld. It will start to decrease when the average displacement Z exceeds the persistence length, the observed maxima in Figs. 11 a and 11 b being a consequence of the combined effect of Brownian motion and ow. For very large displacements, a complete loss of correlation is expected as the dispersion in both X and Z becomes determined by the pseudodiffusive character of the displacements. Note that although the latter condition is almost met for the largest displacements in the simulation, which is indicated by a Gaussian-like shape of the marginal propagators, the correlation coefcient has not yet fallen close to zero. The peak in the correlation coefcient is remarkably well pronounced in the experimental data, given the fact that a polydispersity of the real glass beads exists. However, while the absolute values of in the experimental and in the simulated data sets are comparable for the largest ow rate, the dependence on the Peclet number is more pronounced in the experiment. Given that D, v * , and the bead diameter are identical, and that edge effects are not expected to play any role for the displacements observed, deviations of the simulated from the real velocity eld must be assumed which affect small ow rates most. This feature might be connected to an observation of the peak near zero displacement that is more persistent in the simulation than in the experimental results: a larger fraction of particles possessing both small X and small Z has to inuence the absolute value of the correlation coefcients. The observation of a positive correlation between X and Z suggests the concept of a preferential ow direction, or, more precisely, a direction which is more probable to nd than others. For example, the probability density of particles remaining near X 0 decreases for larger displacements Z. In order to visualize the geometrical properties of this correlation, we write the joint probability density as the sum of the product of its marginals and a function that incorporates the correlation between the measured quantities: P X,Z P X P Z C X,Z . 26

The correlation matrix C (X,Z) contains only zero elements if displacements in X and Z are mutually independent. This case is realized for Brownian motion in an isotropic system where P X,Z 1 2
X Z

exp

X X 0 2/ 2
2 Z

2 X

exp

Z Z 0 2/ 2

However, if X and Z are correlated, a plot of P(X,Z) P(X) P(Z) renders the correlation matrix directly. This operation has been performed for the propagators in Fig. 4, and is presented in Fig. 12 for the experimental data. Simulated data result in equivalent plots. Compared to the actual propagators in Fig. 4, the difference propagators show a much larger amount of structure. One observes a fourfold symmetry: a positive peak for small X and Z develops into a narrow negative region for large Z and small X as well as for large X and small Z. For intermediate values of X and Z, the sign of C (X,Z) is again positive, and its local maximum

6218

STAPF, PACKER, GRAHAM, THOVERT, AND ADLER

PRE 58

FIG. 12. Difference propagator C (X,Z) P (X,Z) P (X) P (Z) for experimental data at a ow rate of 42.0 ml/min, axis scales in m. Contour lines are drawn from approximately 0.05 of the peak intensity in linear spacing; numbers indicate probability densities in 104 m 2 using the normalization of Fig. 4. Positive values are indicated by solid lines, and negative values by dashed lines. Dotted lines represent C (X,Z) 0. Data are cut at points where P (X,Z) is typically below 0.02 of its peak value to avoid contributions by noise; this determines the outer bounds of the C (X,Z) 0 line.

describes a line that is curved in the ow direction. These general features are independent of the parameters in the experiment or the simulation and have been observed before for ow in porous sandstone rocks 27 . It is also found for short-time simulations for average displacements down to a few micrometers. However, for a diffusion coefcient reduced by a factor of 104 in the simulations, the features become much sharper, revealing the initial state of the velocity eld before spreading due to self-diffusion or changes in local ow directions lead to a broadening in the correlation matrix picture. The absolute magnitude in the correlation plot, on the other hand, decreases with longer encoding times see numbers at contour lines in Fig. 12 due to the general spreading of the probability density function. This is analogous to the spreading of the propagator P (X,Z) itself. The correlation matrix can be understood as a function that identies regions of correlations larger or smaller than the reference state P (X) P (Z) which formally denes the independence of the variables X and Z. However, as pointed out by Wadsworth and Bryan 42 , while the use of a function such as C (X,Z) gives a clear answer to whether the variables X and Z are correlated or not, the mathematical representation of the relationship, once demonstrated, is an art in itself. That said, it is clear that, because of normalization of the propagators in general, the integral over C (X,Z) must be zero, which necessitates a particular sym-

metry and a sign change at certain coordinates. A simplied view of C (X,Z) might lead to an interpretation of the ridges of positive values as the lines of preferential displacements. One would then conclude that uid molecules would preferably ow along a certain angle max relative to the X axis, while, for 0 and 90, displacements are less likely to be found. One possible denition of max employs determining the maximum of the radial rst moment, R( ) , with respect to the origin at X,Z 0. While max is not found to be signicantly inuenced by the ow rate for the case investigated in this study, it is certainly a function of the ratio of the average displacements in z and x directions. We are currently considering in more detail the signicance of C (X,Z) and the role of the reference state P (X) P (Z) in its denition and physical interpretation.

E. Comparison with model systems

Another approach towards a better understanding of the meaning of correlations between displacements can be obtained by performing independent simulations and assuming model propagators. One case which can be solved analytically assumes a Gaussian propagator along X for each point in Z, the variance of which scales with a power of Z:

PRE 58

SPATIAL CORRELATIONS AND DISPERSION FOR . . .

6219

P X,Z

4
1 /2

1 2

Z2/

X2/

Z2

27

This propagator fullls the normalization condition P(X,Z)dX dZ 1 when dened for positive X and Z only. The general moments are represented by

n/2

m/2

/2 n

X nZ m

n 1 2 1 2

m 2

1 2 , 28

the correlation coefcients are given by

n 1 2
X n ,Z m

m 2 2
2

m n n 1 2 2
2

2 29 m
2

1 1 n 2 n

m 2

where ( 1)/2. For 0, all correlation coefcients are equal to zero, as expected. Note that the correlations do not depend on the parameters and . The dependence of are shown in Fig. 13. X,Z and X 2 ,Z on Numerical computations of this model propagator have been performed employing a MatLab routine; the resulting moments and correlations were found to match with the theoretical values within less than 1%, depending on the resolution of the calculation. The ndings in Sec. VI D, shown in Fig. 10, suggest 2 0.74. The propagator and the correlation matrix for this parameter are shown in Fig. 14. The similarity to the bead system investigated in this study is striking at rst sight. The main difference between model and experiment lies in the short-displacement behavior. The marginal P(Z) of the distribution of Eq. 27 has the form

Z e Z . Although similarities between the experimentally obtained propagator P(Z) and this distribution might exist, the persistent peak of high probability density near the origin cannot be represented. Also, the expected Gaussian shape of P(Z) for long times is not realized by the model. Nevertheless, the general properties of the correlation matrix coincide reasonably well with both the experimental and simulated data presented above. Tentative analyses of ow through other porous media lead to similar characteristics, but with an exponent 2 considerably different from the one obtained for packed glass beads 43 . A model propagator in this or a modied form can therefore serve for the discussion and comparison of general features of the spreading process in X and Z, while certain details such as the small-displacement and long-time behavior contain additional information about the characteristics of the specic system. In an attempt to gain some further insight into the inuence of local pore geometry on the form of the joint probability density, we have analyzed the probability densities generated by the simple two-bond capillary model, described earlier, for a series of values of the angle in the range 0 80 with v max( 0) 2l. Clearly the dependency of the probability density on will be conned to regions where the displacement exceeds the length l of the initial capillary; as the probability density is small in those regions, we expect the linear correlation coefcient and rst moment to be relatively insensitive to variations in . We have therefore calculated the ratio (X X ) 2 /(Z Z ) 2 as a function of . We nd for this model that

FIG. 13. Correlation coefcients X,Z and X 2 ,Z for the model propagator described in the text Eq. 27 as a function of the parameter .

X Z

X Z

2 2

0.287 0.178 1 cos

30

6220

STAPF, PACKER, GRAHAM, THOVERT, AND ADLER

PRE 58

FIG. 14. a Two-dimensional propagator P(X,Z) for the model described in the text Eq. 27 ; 2 0.74. Contour lines are drawn from 0.05 of the peak intensity in linear spacing. b Correlation matrix C(X,Z) P(X,Z) P(X) P(Z) for the plot in a . Contour lines are drawn from 0.005 of the peak intensity in nonlinear spacing. Positive values are indicated by solid lines, and negative values by dashed lines. Dotted lines represent C(X,Z) 0.

To compare this result to our experimental and simulated data, it is necessary to relate the bond length l to the characteristic length scale in our sample. We take l L/3; this choice, although inevitably to some extent arbitrary, is physically reasonable and gives values for Z and (Z Z ) 2 in the ranges 5865 m and 4000 6100 m2, respectively, consistent with the data presented in Fig. 7 a . The second-moment ratios for our experimental and beadpack simulation data are almost independent of Z in the range of Z under consideration. They do however depend upon Peclet number, ranging from 0.5 for Pe 88 to an apparent limit of 0.4 for Pe 250, when spreading is dominated by convective transport. The limiting value is the most appropriate for comparison with our model, which neglects effects due to diffusion; Eq. 30 then yields a value of 70 for . That value is not unreasonablefor example, in a hexagonal close-packed structure uid elements experience a change in direction of 72 between entering and leaving a tetrahedral pore. However, in view of the simplicity of our model and the fact that random packing must lead to a distribution of pore geometries, we must regard this result as merely illustrative of the sort of pore scale structural information which could be extracted from the form of the joint density at small displacements.

determined experimentally for water ow in a system of packed glass beads. Numerical simulations have been carried out for a wider range of average displacements than was accessible by experiment. The simulations, obtained by successively generating a random packing of spherical beads, solving the Stokes equation, and monitoring of the displacements of particles randomly distributed in the pore space, have been shown to agree reasonably well with the experimental results. The major part of the deviations can be explained by a different spreading behavior of a fraction of particles remaining near their starting positions. The evolution of second moments of the displacement distributions as a function of ow rate and encoding time shows the crossover from a regime dominated by Brownian motion to an asymptotical approach toward propagators of Gaussian shape. The limiting dispersion for long times is reached much earlier in X compared to Z. The dependence of the dispersion coefcient D* on the Peclet number suggests an intermediate behavior where the mixing process is dominated by velocity variations but contains a nonnegligible contribution from self-diffusion. The mean-squared displacement perpendicular to the ow direction follows a power-law dependence on Z for large Z. From the full joint probability propagator representation, a correlation between X 2 and Z is found that decreases for average displacements in the ow direction Z , larger than a characteristic length equivalent to the bead radius, while decreasing again for much smaller Z due to the dominating effect of isotropic self-diffusion. A plot of the correlation matrix C (X,Z) P (X,Z) P (X) P (Z) possibly indicates preferential directions for displacements by larger-than-average correlation values for a specic region of displacements. The shape of this region changes with increasing ow time. Correlation matrices of similar shape can be generated by appropriate model propagator functions. Numerical simulations have shown that the peculiarities of the matrix shape reect characteristic structural details such as size and orientation of microcapillaries. Experiments and theoretical analyses are currently being performed to describe a wider range of porous systems by the properties of uid ow through their pore space 43 . Twodimensional propagators can be employed to improve the understanding of how uids spread in a porous network and how this spreading is governed by parameters such as self-diffusion coefcients and ow rates, on the one hand, and porosity, permeability, and tortuosity on the other hand.

ACKNOWLEDGMENTS

VII. CONCLUSIONS

The joint probability distribution for displacements parallel and perpendicular to the net ow direction have been

One of the authors S.S. is indebted to the Deutsche Forschungsgemeinschaft DFG , Grant No. Sta 511/1-1 and the Training and Mobility of Researchers TMR Programme of the European Union Grant No. ERBFMBICT961836 , for nancial support.

PRE 58

SPATIAL CORRELATIONS AND DISPERSION FOR . . .

6221

1 J. J. Fried and M. A. Combarnous, Adv. Hydrosci. 7, 169 1971 . 2 A. L. Dullien, Porous Media: Fluid Transport and Pore Structure Academic, New York, 1979 . 3 J. H. Cushman, Dynamics of Fluids in Hierarchical Porous Media Academic, New York, 1990 . 4 H. D. Pfannkuch, Rev. Inst. Fr. Petrole 18, 215 1963 . 5 D. J. Gunn and C. Pryce, Trans. Inst. Chem. Eng. 47, 341 1969 . 6 R. Aris, Proc. R. Soc. London, Ser. A 235, 67 1956 . 7 G. I. Taylor, Proc. R. Soc. London, Ser. A 219, 186 1953 . 8 D. L. Koch and J. F. Brady, J. Fluid Mech. 154, 399 1985 . 9 D. A. Edwards, M. Shapiro, H. Brenner, and M. Shapira, Transp. Porous Media 6, 337 1991 . 10 H. P. Amaral Souto and C. Moyne, Phys. Fluids 9, 2253 1997 . ` 11 J. Salles, J.-F. Thovert, R. Delannay, L. Prevors, J.-L. Auriault, and P. M. Adler, Phys. Fluids A 5, 2348 1993 . 12 E. Charlaix, J.-P. Hulin, and T. J. Plona, Phys. Fluids 30, 1690 1987 . 13 A. Ding and D. Candela, Phys. Rev. E 54, 656 1996 . 14 D. L. Koch and J. F. Brady, J. Fluid Mech. 180, 387 1987 . 15 D. L. Koch and J. F. Brady, Chem. Eng. Sci. 42, 1377 1987 . 16 P. T. Callaghan, Principles of Nuclear Magnetic Resonance Microscopy Clarendon, Oxford, 1991 . 17 J. Karger and W. Heink, J. Magn. Reson. 51, 1 1983 . 18 K. J. Packer and J. J. Tessier, Mol. Phys. 87, 267 1996 . 19 M. H. G. Amin, S. J. Gibbs, R. J. Chorley, K. S. Richards, T. A. Carpenter, and L. D. Hall, Proc. R. Soc. London, Ser. A 453, 489 1997 . 20 J. D. Seymour and P. T. Callaghan, AIChE. J. 43, 2096 1997 . 21 H. Van As and D. van Dusschoten, Geoderma 80, 389 1997 . 22 J. J. Tessier and K. J. Packer, Phys. Fluids 10, 75 1998 . 23 L. Lebon, L. Oger, J. Leblond, J. P. Hulin, N. S. Martys, and L. M. Schwartz, Phys. Fluids 8, 293 1996 . 24 J. J. Tessier, K. J. Packer, J.-F. Thovert, and P. M. Adler, AIChE. J. 43, 1653 1997 .

25 R. S. Maier, D. M. Kroll, Y. E. Kutsovsky, H. T. Davis, and R. S. Bernard, Phys. Fluids 10, 60 1998 . 26 D. Coelho, J.-F. Thovert, and P. M. Adler, Phys. Rev. E 55, 1959 1997 . 27 K. J. Packer, S. Stapf, J. J. Tessier, and R. A. Damion, Magn. Reson. Imaging to be published . 28 E. J. Fordham, S. G. Gibbs, and L. D. Hall, Magn. Reson. Imaging 12, 279 1994 . 29 E. O. Stejskal and J. E. Tanner, J. Chem. Phys. 42, 288 1965 . 30 A. J. Lucas, S. G. Gibbs, E. W. Jones, M. Peyron, A. D. Derbyshire, and L. D. Hall, J. Magn. Reson., Ser. A 104, 273 1993 . 31 IDL Interactive Data Language, Research Systems Inc., Boulder 1995 . 32 R. Lemaitre and P. M. Adler, Transp. Porous Media 5, 325 1990 . 33 S. R. McDougall and K. S. Sorbie, Petroleum Geosci. 3, 161 1997 . 34 J. Bear, Dynamics of Fluids in Porous Media Academic, New York, 1975 . 35 R. Kimmich, S. Stapf, A. I. Maklakov, V. D. Skirda, and E. V. Khozina, Magn. Reson. Imaging 14, 793 1996 . 36 L. L. Latour, P. P. Mitra, R. L. Kleinberg, and C. H. Sotak, J. Magn. Reson., Ser. A 101, 342 1993 . 37 A. Compte, R. Metzler, and J. Camacho, Phys. Rev. E 56, 1445 1997 . 38 N.-W. Han, J. Bhakta, and R. G. Carbonell, AIChE. J. 31, 277 1985 . 39 P. M. Adler, Porous Media. Geometry and Transport Butterworth-Heinemann, Stoneham, 1992 . 40 D. R. F. Harleman and R. R. Rumer, J. Fluid Mech. 16, 385 1963 . 41 Y. E. Kutsovsky, L. E. Scriven, H. T. Davis, and B. E. Hammer, Phys. Fluids 8, 863 1996 . 42 G. P. Wadsworth and J. G. Bryan, Introduction to Probability and Random Variables McGraw-Hill, New York, 1960 . 43 S. Stapf and K. J. Packer unpublished .