Signature of Non-Fickian Solute Transport in Complex Heterogeneous Porous Media

We simulate transport of a solute through three-dimensional images of different rock samples, with resolutions of a few microns, representing geological media of increasing pore-scale complexity: a sandpack, a Berea sandstone, and a Portland limestone. We predict the propagators (concentration as a function of distance) measured on similar cores in nuclear magnetic resonance experiments and the dispersion coefficient as a function of Péclet number and time. The behavior is explained using continuous time random walks with a truncated power-law distribution of travel times: transport is qualitatively different for the complex limestone compared to the sandstone or sandpack, with long tailing, an almost immobile peak concentration, and a very slow approach to asymptotic dispersion.

Introduction.-Despite the range of significant applications of solute transport, including improved oil recovery, the long-term fate of nuclear waste repositories, and secure storage of carbon dioxide [1], even the qualitative behavior of most rocks is uncertain.Vast carbonate sedimentary basins contain more than half the world's current oil reserves [2] yet experimental data on dispersion and transport in carbonates is scant.Although driven by the interplay of two simple physical processes-advection and diffusion-the resultant macroscopic behavior is surprisingly rich.Experiments on packed beds, sand columns, sandstone and carbonate cores [3][4][5][6][7][8][9][10], and field studies [11,12] have shown that transport is characterized by early breakthrough of the solute, a long tailing of the concentration at late times and-in many cases-a virtually immobile peak concentration.Transport cannot be described by a traditional advection-dispersion formulation; instead other theories, including multirate mass transfer models and continuous time random walks (CTRW), have been employed to describe the evolution of a dissolved species over time [13,14].Network modeling has been used to provide an explanation for the power-law dependence of dispersion coefficient as a function of Pe ´clet number (Pe ¼ u av L=D m , where u av is the average flow speed, L is the characteristic length, and D m is the molecular diffusion coefficient), reconciling experiment, pore-scale modeling, and CTRW theory for Berea sandstone [15,16].However, for more complex porous media the relationship between pore structure, velocity field, and transport remains unknown, particularly for heterogeneous carbonates.We compute transport through micro-CT images of the pore space with a streamline-based random-walk simulator and demonstrate the qualitatively different signature of transport through the major porous rock types encountered in the subsurface: sandpacks, sandstones, and carbonates.
Model.-We simulate incompressible steady viscous flow directly on pore-space images obtained from micro-CT scanning by solving the Stokes equation r 2 v ¼ rP subject to no-slip velocity boundary conditions.v is the velocity, P is pressure, and is the fluid viscosity.We use a standard predictor-corrector method [17] in conjunction with the finite difference technique using grid blocks obtained from a binarized pore-space image.We apply a streamline-based algorithm with a new semianalytic formulation near solid boundaries.We satisfy r Á v ¼ 0 with a parabolic variation of velocity perpendicular to a solid boundary and bilinear dependence parallel to the solid within each grid block [18].An ensemble of particles is moved by advection along streamlines; we calculate the time needed for a particle inside a pore-space voxel to exit, Át.The time step is dt.If Át > dt the particle remains within the voxel, whereas for Át < dt the particle enters a new block.A random-walk method is used to capture molecular diffusion: a particle instantaneously jumps over a mean-free path ¼ ffiffiffiffiffiffiffiffiffiffiffi ffi 6D m t p in a random direction.We apply a reflection boundary condition on the surface of the solid voxels.The flow domain is cubic.We impose a uniform pressure drop (1.172 kPa for a typical velocity of 1:3 mm=s) across two opposing faces and use periodic boundary conditions for the other faces.If a particle exits the inlet or outlet, it is randomly reassigned to the opposite face-flux-weighted during the advective step and area weighted for the diffusive step [16].The dispersion coefficient D L is computed by calculating the variance of the distance traveled by particles in the main direction of flow , where 2 is the variance of the particle displacement, x i , 2 ¼ h½x i ðtÞ À hx i ðtÞi 2 i.The time step for the simulation is 10 À4 s.The average motion of particle at each time step is less than one voxel.The base-case value of D m was the free molecular self-diffusion coefficient of water, 2:2 Â 10 À9 m 2 s À1 .We transport 50 000 particles that are initially uniformly distributed throughout the pore space.We study transport in three representative geological porous media: a sandpack comprised of a close packing of irregular quartz grains of different size, with a porosity of 37.3% and permeability of 3:9 Â 10 À11 m 2 ; a Berea sandstone that is widely used as a benchmark material representing a large class of rocks in relatively homogeneous oil reservoirs and deep aquifers [19] with a porosity of 19.62% and permeability 1:4 Â 10 À12 m 2 ; and a Portland limestone, with a heterogeneous microstructure typical of many reservoir carbonates [3] with a porosity of 9.17% and permeability 3:0 Â 10 À13 m 2 .The porosity and permeability values were computed on the images.We used micro-CT images of the sandpack, Berea sandstone, and Portland limestone at a resolution (voxel size) of 10 m, 5:345 m, and 9 m, respectively, that have been binarized into solid and void and which contain 300 3 grid blocks (voxels) in total [19] representing a cube of side length 1.6-3.0mm.
In CTRW c ðtÞdt is defined as the ensemble averaged probability of a particle just arrived at a site moving to an adjacent site in a time between t and t þ dt.We compute c ðtÞ on a grid-block basis by finding the time it takes for each particle to traverse each void cell.This approach has been used before using a pore-scale network model where the transitions are from pore to pore [16] and for grid blocks in field-scale simulation [20].We compute the distribution of transit times that particles spend in grid blocks as a function of dimensionless time b ¼ t=t b1 , where t b1 ¼ Áx=u av is the mean advective transit time whose characteristic length is the grid-block spacing Áx, indicated by the subscript b.The average velocity u av is defined as the total flow across the model divided by the cross-sectional area.The travel times are sampled for 1 s or b up to 145.
Results.-In Fig. 1 we show the velocity field for Portland carbonate that illustrates that flow is concentrated in a few channels with much of the pore space largely stagnant.We also plot c ð b Þ for a range of Pe for Portland carbonate and Pe ¼ 1 for sandpack and Berea sandstone.We define the characteristic length L in Pe based on a cubic packing of regular spheres.For this idealized system, the grain diameter is V=S, where V is the volume of the porous medium (pore plus grain) and S is the area of the pore-grain interface.We use the same definition for our images, since the volume and the pore and grain area are readily computed, while grain size is unknown, and obtain the values 160, 130, and 290 m for the sandpack, sandstone, and carbonate, respectively: the images span around 9-19 characteristic lengths.
For the late-time behavior of c ð b Þ, b > 1, we observe an approximately power-law dependence of travel times c $ Àð1þÞ b with ¼ 1:8 for the sandpack and sandstone and ¼ 0:7 for the carbonate.A lower value for indicates greater heterogeneity of the porous medium, resulting in a wider range of transit times leading to more complex non-Fickian transport behavior.The power law is truncated at late times, representing the time to diffuse across grid block (at low to intermediate Pe) or an advective cutoff for large Pe, representing the slowest flow speeds.The power-law dependence has been noted previously for a network model of Berea sandstone [16]: we see a power-law dependence < 1 for the carbonate, indicating qualitatively different transport, as discussed later.
We analyze preasymptotic transport by comparing experimental results for displacement probabilities (propagators) with pore-scale numerical simulation.We show here only the results for the carbonate-we also obtain good predictions for the sandstone and sandpack.In Fig. 2 we compare our results to two independent measurements on Portland stone using NMR [3,4].It is evident that the peak concentration is virtually immobile with a highly dispersed fast-moving tail: particles reside for a long time in slow or no-flow regions, close to the solid, in dead-end pores or in narrow pore spaces, eventually diffusing out and moving rapidly through the better connected, wider regions.Our predictions are in excellent agreement with both sets of experiments.The immobile peak concentration is a  signature of transport with < 1 [14] that is not observed in sandstones and sandpacks [3,4].
Next we compute how the longitudinal dispersion coefficient D L of the injected plume evolves with time for Portland carbonate.Figure 3 shows results for Pe ¼ 2, Pe ¼ 20, Pe ¼ 200, and Pe ¼ 2000.
From CTRW, if c $ t Àð1þÞ then 2 $ t 2 , for < 1; for 2 > > 1, 2 $ t 3À , while for > 2 the behavior is Fickian with 2 $ t [21].There is a late-time cutoff to the power law at time t 2 .For low to moderate Pe this is controlled by diffusion across the characteristic length (a typical pore size) and t 2 ¼ L 2 =2D m [16].For t > t 2 we see a crossover to asymptotic (Fickian) behavior.We now define dimensionless parameters using the typical length L as scaling factor: ¼ t=t 1 where t This behavior is evident in Fig. 3 with a power-law exponent in the preasymptotic regime 2 À 1 ¼ 0:4 consistent with ¼ 0:7 (Fig. 1), even for Pe ¼ 2 and Pe ¼ 200 where a clear power-law regime is not evident.When t ) t 2 the asymptotic dispersion coefficient is defined and the traditional Fickian description is valid.
In Fig. 4 we plot the values for asymptotic dispersion coefficient as a function of Pe for the carbonate, sandstone, and sandpack studied.Plotted also are the experimental data for unconsolidated beadpacks and sandpacks in the literature [16,22]; experimental data for carbonates that cover the whole range of Pe for a well-characterized sample are not available.Good agreement is obtained between the experiments and our model for sandpack and Berea sandstone.In the restricted diffusion regime (Pe < 0:1) D L =D m has lower value for the carbonate than for the sandstone (that itself has a lower value than the sandpack) due to higher tortuosity in the carbonate.In the power-law regime for the carbonate (5 < Pe < 100) we observe the scaling D L =D m $ Pe 1:4 , while in the power-law regime for the sandpack and sandstone (10 < Pe < 600) the scaling is D L =D m $ Pe 1:2 .This corresponds to the values of ¼ 0:7  [16,22].In the power-law regime for the carbonate we observe the scaling D L =D m $ Pe 1:4 , while for the sandpack and sandstone this scaling is D L =D m $ Pe 1:2 .This corresponds to the values of ¼ 0:7 for carbonate and ¼ 1:8 for sandpack and sandstone (as seen in Fig. 1) and agrees with the predictions from CTRW theory.
for the carbonate from Fig. 1 and ¼ 1:8 for the sandpack and sandstone from Fig. 1 (see also Refs.[16,18,23]).Moreover, these results agree with the predictions from CTRW theory presented above.We note that the magnitude of dispersion coefficient in the power-law regime increases as the medium heterogeneity becomes more complex-from sandpack to sandstone to carbonate, as seen in Fig. 4. For the highest values of Pe we see a transition to a linear scaling of D L with Pe, consistent with an advective cutoff to the transition time probability c ðtÞ [16].Discussion and conclusions.-Weuse a direct streamline-based simulation method on pore-space images to study transport in three representative geological media of increasing complexity.Transport in carbonates has a rich non-Fickian behavior, indicated by NMR experiments, that is generically different from sandstones, with a very slowmoving peak concentration, a highly dispersed tail, an extremely gradual approach to an asymptotic regime, and a nonlinear scaling of dispersion coefficient with Pe, consistent with a CTRW analysis of the grid-block travel times.
In carbonates the peak plume position will be retarded relative to the mean flow field with a very wide spread.This may suppress convective mixing of solute and lead to long residence times for tracers and pollutants.For instance, during carbon dioxide, CO 2 , storage, some of the CO 2 dissolves in brine and this denser CO 2 -laden brine sinks.This motion allows the development of a gravitational instability, further mixing, and eventual secure storage.Similarly, mixing during miscible gas injection controls local oil recovery in hydrocarbon fields.However, recent analytical and numerical analyses of these problems [24,25] presume Fickian dispersion.We suggest that this is almost certainly qualitatively incorrect for carbonates and other highly heterogeneous porous media leading to overestimates of the rate of phase exchange in CO 2 storage and enhanced oil recovery.

FIG. 1 (
FIG. 1 (color).The probability c ð b Þ of traveling between two neighboring voxels for Portland carbonate, for Pe ¼ 2, 200, 2000, and 1 (symbols).The probabilities show a power-law trend c $ Àð1þÞ b ; a line with slope corresponding to ¼ 0:7 is represented.The dimensionless time is b ¼ t=t b1 , where t b1 is the mean travel time through a grid block Áx=u av .Shown in the left inset is the dimensionless velocity field computed in the Portland carbonate image; in the right inset the comparison of c ð b Þ for Portland carbonate (black symbols) with c ð b Þ of the sandpack (light blue symbols) and Berea sandstone (red symbols) for Pe ¼ 1 indicates quantitatively different generic behavior characterized by the two parallel solid lines (dark blue and orange color) with slope corresponding to ¼ 1:8.

PRL 107 ,
204502 (2011) P H Y S I C A L R E V I E W L E T T E FIG. 2. Probability of displacement in a micro-CT image of Portland carbonate PðÞ as a function of displacement (solid line), compared with the propagators obtained from NMR experiments from Scheven et al. [3] and Mitchell et al. [4] for the same time t ¼ 1 s and flow rate.The coordinates are rescaled by the nominal mean displacement hi 0 ¼ u av t. u av ¼ 1:3 mm=s and Pe ¼ 171 in our simulation.FIG. 3. The preasymptotic behavior of D L for Portland carbonate computed from the pore-scale model as a function of dimensionless time, ¼ t=t 1 where t 1 ¼ L=u av , is plotted for different values of Pe The solid lines with the slope 0.4 indicate the predicted scaling using Fig. 1 and CTRW: D L =D m $ Pe Â 2À1 , where the slope is 2 À 1 ¼ 0:4 for ¼ 0:7.