Predictions of non-Fickian solute transport in different classes of porous media using direct simulation on pore-scale images

We present predictions of transport through micro-CT images of porous media that include the analysis of correlation structure, velocity, and the dynamics of the evolving plume. We simulate solute transport through millimeter-sized three-dimensional images of a beadpack, a sandstone, and a carbonate, representing porous media with an increasing degree of pore-scale complexity. The Navier-Stokes equations are solved to compute the flow field and a streamline simulation approach is used to move particles by advection, while the random walk method is employed to represent diffusion. We show how the computed propagators (concentration as a function of displacement) for the beadpack, sandstone, and carbonate depend on the width of the velocity distribution. A narrow velocity distribution in the beadpack leads to the least anomalous behavior, where the propagators rapidly become Gaussian in shape; the wider velocity distribution in the sandstone gives rise to a small immobile concentration peak, and a large secondary mobile peak moving at approximately the average flow speed; in the carbonate with the widest velocity distribution, the stagnant concentration peak is persistent, with a slower emergence of a smaller secondary mobile peak, characteristic of highly anomalous behavior. This defines different types of transport in the three media and quantifies the effect of pore structure on transport. The propagators obtained by the model are in excellent agreement with those measured on similar cores in nuclear magnetic resonance experiments by Scheven, Verganelakis, Harris, Johns, and Gladden, Phys. Fluids 17, 117107 (2005).


I. INTRODUCTION
Solute transport in porous media is of importance in a broad range of scientific fields and engineering applications, notably for contaminant migration in subsurface hydrology [1,2], geological storage of carbon dioxide [3], packed bed reactors and chromatography in chemical engineering [4,5], and tracer studies and miscible displacement in enhanced oil recovery [1,6].A considerable body of research has been devoted to describing transport that results in a nonlinear growth of the variance of displacement with time.These processes cannot be described by the solution of the governing transport equations in a homogeneous medium with Fick's law used to describe the dispersive flux; hence this type of transport is said to be non-Fickian [2].The interplay of simple physical processes-advection and diffusion-results in a rich macroscopic transport behavior that is a consequence of the wide range of local flow speeds experienced by the moving particles, combined with local mixing (see recent reviews [2,7] and references therein).The non-Fickian nature of dispersive processes in heterogeneous porous media has been demonstrated experimentally from pore to field scales [8][9][10][11][12][13][14][15][16].However, the predictive understanding of the relationship between pore structure, the velocity field, and transport is still limited.Our goal is to investigate and explain the origin of non-Fickian transport behavior as a function of pore-scale heterogeneity by simulating flow and solute transport directly on micro-CT images of pore space in simple homogeneous beadpacks and more complex geological porous media, namely sandstone and carbonate rocks.
In the asymptotic limit, assuming that the porous medium is homogeneous, when the velocity field is fully sampled by solute, macroscopic transport parameters, such as the dispersion coefficient D, are constant and can be used in an averaged advection-dispersion equation.However, until the velocity field is fully sampled, transport is non-Fickian and D possesses temporal or spatial variation.This variation is observed in varying plume shapes and the corresponding description of the probability density function (PDF) of either displacement or transit times of the solute particles.These PDFs have been studied experimentally by nuclear magnetic resonance (NMR) measurements on unconsolidated beadpacks [12,[17][18][19][20][21] and samples of geological rock [12,[22][23][24][25][26][27].Notable progress has been made in matching the experimental NMR propagator data with numerical models that either simulate transport directly on the pore space of the beadpacks [19,20,28], sandstones [22,29], and a carbonate [30], or use pore networks extracted from micro-CT images in a sandstone [31] and dolomite [27].
A landmark comparative experimental NMR study of propagators in a beadpack, a sandstone, and a carbonate rock was presented by Scheven et al. [12].This study has demonstrated that the propagators measured in the beadpack rapidly reach a symmetric Gaussian-like shape about the mean displacement, consistent with Fickian transport; however, for Bentheimer sandstone and Portland carbonate two asymmetric peaks for stagnant and mobile fluid are observed at different observation times.In Bentheimer sandstone, at early times there is a pronounced peak representing the stagnant fluid regions that gradually disappears over time with the emergence of a dispersed plume of mobile fluid approximately centered about the mean displacement.A much larger stagnant peak is observed for Portland carbonate that persists for longer times-consequently the formation of a highly dispersed mobile plume is delayed.
No predictive model has been able to explain the relationship between pore structure, velocity field, and transport.In this study we examine the non-Fickian behavior of solute transport simulated directly on micro-CT images of the pore space.We simulate solute transport through three-dimensional images of a beadpack (a disordered close packing of spheres), Bentheimer sandstone, and Portland carbonate, representing porous media with an increasing degree of pore-scale complexity.The Navier-Stokes equations are solved to compute the flow field.A streamline method is used to represent advection, while a random walk approach models diffusive transport.
In earlier work, we have used this approach to study transport in a sandpack, sandstone, and carbonate to compute the distribution of transit time across image voxels and to predict the asymptotic dispersion coefficient [29,30].We showed good agreement between NMR-measured and predicted propagators at one time for Bentheimer sandstone and Portland carbonate [29,30].Herein, we extend this work to provide a more detailed and accurate comparison with experiment at all the times measured for the three media listed above.Moreover, we explain the results in terms of the distribution of velocities in the pore space.We present a methodology for making transport predictions on micro-CT images that includes analyses of variograms of porosity and velocity, the velocity fields and velocity distributions, and the dynamics of the evolving plume.

II. IMAGES, MATHEMATICAL MODEL, AND FLOW FIELDS
Our transport simulations require a description of the pore space, a method for calculating the flow field, and a particle tracking method for moving solute by advection and diffusion.

A. Pore-space images
We study transport in the three porous media mentioned in the Introduction, with increasing pore-scale complexity.First is an experimental beadpack image with 300 3 voxels comprised of a disordered close packing of spherical grains of the same size [32,33] with diameter d g = 100 μm; the porosity is 35.93% and the permeability is 5.63 × 10 −12 m 2 .Second is a Bentheimer sandstone 300 3 micro-CT image with a porosity of 21.51% and permeability 3.53 × 10 −12 m 2 .Third is a Portland limestone 320 3 x-ray synchrotron image with a porosity of 8.62% and permeability 2.47 × 10 −13 m 2 .The porosity and permeability values were computed on the images.The porosity is calculated as the ratio of the number of pore voxels, N pvox divided by the total number of voxels, N vox .The voxels that have no connection through the pore space to either the inlet or outlet are excluded from the flow calculations and further analysis.
The voxel sizes are 2, 3, and 9 μm for the beadpack, Bentheimer sandstone, and Portland carbonate, respectively.The images have been binarized into solid and void, which for the above-mentioned image sizes, represent cubes with side lengths of 0.6, 0.9, and 2.88 mm for the beadpack, Bentheimer sandstone, and Portland carbonate, respectively.Figures 1(a)-1(c) show two-dimensional cross sections of the segmented image of the beadpack and the gray scale images for Bentheimer sandstone and Portland carbonate.
The dry scan image for Bentheimer sandstone was acquired on a cylindrical core of 5 mm diameter and length 25 mm with an Xradia Versa micro-CT scanner (provided by iRock Technologies).The dry scan image for Portland carbonate was acquired on a cylindrical core of the same size with a synchrotron beamline (SYREMP beamline at the Elettra synchrotron in Trieste, Italy) at a resolution of 9 μm, corresponding to two different detector pixel sizes of 3.85 and 4.5 μm; the CCD camera binned the results giving the final resolution of twice the detector pixel size.Reconstruction and image analysis was performed by in-house software, resulting in images of around 600 3 pixels from which a central cubic section was taken for our simulations.

B. Mathematical model for flow
For calculating the flow field we use a standard finite volume method implemented in OPENFOAM [34].Incompressible steady viscous flow is simulated directly through the porespace images by solving the volume conservation equation ( 1) and the Navier-Stokes equations (2): where u is the velocity vector, μ is the viscosity of water (μ = 0.001 Pa s), ρ is the density of water (ρ = 1000 kg/m 3 ), and p is pressure.The pressure and velocity are solved iteratively based on the pressure implicit with splitting of operators (PISO) algorithm of Issa [35]; further details can be found in [36].We impose no flow, u = 0, on solid boundaries.
The simulations are run at a Reynolds number Re = ρu av L μ 1 assuming steady state, ∂ u ∂t = 0.The average flow speed is calculated as u av = q/ε, where q = Q/(L y L z ) is the Darcy velocity, Q (m 3 /s) is the total volumetric flux calculated as Q = ∫ u x dA x , where A x (m 2 ) is the cross-sectional area of void voxels perpendicular to the direction of flow x, and u x is the face velocity that is normal to A x ; L x ,L y ,L z are the image lengths in each direction, and ε is porosity.Each voxel in the image is converted to a grid block in the finite volume mesh.Since we simulate slow flow, the second term on the left in Eq. ( 2) is small compared to the second term on the right (viscous) term.
The flow domain is cubic.We use constant pressure boundary conditions for pressure at the left and the right faces of the images (the pressure drop is P ).For the other faces of the images and for the solid walls, no-slip boundary conditions are used.By solving the Navier-Stokes equations we obtain the velocities and pressures for each voxel, and calculate absolute permeability k (m 2 ) from Darcy's law: The steps in the flow field computation are presented in Figs.We can define a characteristic length L (related to a typical grain size) that is given by πV /S, where V is the volume of the porous medium (pore plus grain) and S is the area of the pore-grain interface [37].The area S is measured directly on the image from the number of voxel faces separating void from grain.This method is employed for consolidated media where it is not possible to extract an unambiguous mean grain size.Using this method we find characteristic lengths of 139.9 and 327.0 μm for the sandstone and carbonate, respectively.The grain diameter 100 μm is used as the characteristic length for the beadpack.
To study the correlation structure, in Figs.3(a)-3(c) we plot variograms for porosity, γ p , and velocity in the direction of flow, γ u x , for the images of beadpack, Bentheimer sandstone, and Portland carbonate.The functions are calculated as I (x i ) is the indicator function for porosity [I (x i ) = 1 for pore voxels and I (x i ) = 0 for grain voxels], u x (x i ) are velocities in the direction of flow across faces oriented normal to the x direction, and N is the number of voxels.Plotted are γ p and γ u x values normalized to the theoretical values at infinite range (uncorrelated limit) The variograms for both porosity and velocity for all three samples indicate that the system becomes uncorrelated beyond the characteristic length, which is much smaller than the total system size.This suggests that the images are sufficiently large to obtain representative results to compare with experiments on larger core samples.

C. Velocity distributions
Figures 2(g)-2(i) show the very different nature of the velocity fields: in the beadpack flow they are well connected, evenly distributed throughout the sample and characterized by less tortuous channels; in Bentheimer sandstone similar features are observed but with a higher degree of tortuosity; however, in the poorly connected Portland carbonate, flow is concentrated in a few channels with considerable stagnant regions of the pore space.
Furthermore, the velocity distributions obtained on the images reveal considerable differences for the three porous media studied.In Fig. 4 the probability density functions of the ratio of the magnitude of u (at the voxel centers) divided by the average flow speed u av are presented for the beadpack, Bentheimer sandstone, and Portland carbonate as semilog and log-log plots.The analytical probability density function of |u|/u av for a single circular cylindrical tube is also shown to represent the homogeneous limit.The PDF is a histogram of the velocity distribution sampled uniformly in 256 bins of log(|u|/u av ).The PDFs exhibit different characteristics in terms of the spread between low and high velocities, and the magnitude of the peak centered approximately on |u|/u av = 1.
Figure 4 demonstrates that there are over eight orders of magnitude of variation in flow speed, and, unlike for the single tube, a significant fraction of the pore space for beadpack, Bentheimer sandstone, and Portland carbonate experiences very low velocities.The beadpack has a PDF that-for low velocities-is similar to a single tube-again emphasizing the Variograms showing the normalized functions for porosity, γ p /γ p,∞ , and velocity in the direction of flow, γ ux /γ ux,∞ , for the images of beadpack (a), Bentheimer sandstone (b), and Portland carbonate (c).Indicated by the vertical line is the characteristic length (representing a mean grain size) for all three samples.We find that the porosity and velocity appear to be largely uncorrelated beyond this characteristic length, which is much less than the total system size.
homogeneous nature of the system.The spread of velocity in Portland carbonate is considerably wider than in Bentheimer sandstone, which in turn has a distribution that is much wider than the beadpack.In Portland carbonate many velocities are one thousandth of the average or lower.Fewer voxels are effectively stagnant for Bentheimer sandstone, and even fewer for the beadpack.At the fast extreme of the distributions, locally, the velocity can be more than 100 times the average.There is the greatest spread of these fast speeds for Portland, less for Bentheimer, and least for the beadpack.Note also that the peak of the distributions is close to the average, pore velocity (|u|/u av = 1); this peak is largest for the beadpack.We will use these characteristics to interpret the shapes of transport propagators in Sec.III and to explain different generic transport behavior.We will also demonstrate that we are able to make quantitative predictions of experimental results.

D. Transport model
For moving solute by advection, we use a particle tracking method that employs a semianalytic description of the velocity field within a grid block for all combinations of solid boundaries [37].For a known velocity u within a voxel we move solute by a displacement udt in each time step.A random walk method is used to describe molecular diffusion: a particle instantaneously jumps over a mean-free path ξ = √ 6D m t in a random direction.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 diffusion coefficient is 2.2 × 10 −9 m 2 s −1 which is the free self-diffusion coefficient of water [38].Particles are initially placed in uniformly spaced voxels.Within each voxel the particle is placed at random.We place 1 000 000 particles in the pore space.This boundary condition represents the conditions in NMR experiments, where the transport of water through water is measured: The solute is distributed uniformly throughout the medium.We apply a reflection boundary condition for particles that hit the surface of the solid voxels.If a particle exits the inlet or outlet face of the image, it is randomly reassigned to the opposite face-flux weighted during the advective step and area weighted for the diffusive step [39].Reflecting boundary conditions are used for the other image faces.We track particles and plot concentration profiles as a function of particle displacement (propagators).

III. RESULTS AND DISCUSSION
We study the different generic types of non-Fickian transport in porous media by analyzing displacement probabilities (propagators) on the images of the beadpack, Bentheimer sandstone, and Portland carbonate.We use the velocity distributions from Fig. 4 to explain the behavior.Furthermore, we predict accurately the propagators measured in NMR experiments by Scheven et al. [12].

A. Propagators and the effect of pore structure
The evolution of the propagators relative to the expected mean displacement in the main flow direction are presented in Fig. 5 for the beadpack [Fig.5(a)] that has a relatively narrow distribution of velocities [as shown in Fig. 4(b)], for Bentheimer sandstone [Fig.5(b)] with a wider spread of velocities, and Portland carbonate [Fig.5(c)] with the widest spread.The probabilities of displacement are plotted for different times t = 0.106, 0.2, 0.45, 1, and 2 s, as a function of displacement at the same flow conditions as reported in the NMR experiments by Scheven et al. [12]-that is u av = 0.91, 1.03, and 1.3 mm/s for the beadpack, Bentheimer sandstone, and Portland carbonate, respectively.
Pore structure complexity determines quantitatively different generic transport behavior that is evident in Fig. 5 as very different shapes of propagators in the beadpack, Bentheimer sandstone, and Portland carbonate.At early times (t = 0.106 s), in the beadpack, the majority of the solute is moving, although the propagator is asymmetric.This means that the distribution is not Gaussian as would result from a solution of the transport equation in one dimension employing Fick's law.However, as time progresses, solute particles start to sample the entire velocity field, and for t = 0.45 s onwards the propagator rapidly becomes Gaussian about the mean displacement, thus representing Fickian behavior.A very different response is observed at early times for Bentheimer sandstone.A considerable peak in concentration around zero is present-this is stagnant solute that moves by diffusion (note the negative displacement values due to random motion by diffusion).The immobile peak gradually narrows due to mass transfer exchange between immobile and mobile fluid regions by diffusion.At longer times (t = 2 s) there is a formation of a dominant secondary mobile fluid peak in concentration that is centered approximately on 1.For Portland carbonate, a large immobile fluid concentration peak is observed that is persistent with time.The mobile fluid peak starts to form only at late times (t > 1 s) and its magnitude is much lower than that in the beadpack and Bentheimer sandstone.This implies a delayed exchange between flowing and stagnant regions of the pore space.The difference in transport through the three media can be explained by the velocity distributions in Fig. 4. The large stagnant peak in Portland carbonate is a consequence-as mentioned in the previous section-of a large number of very low velocities in the pore space.
The overall dispersion of the concentration is also related to the spread of the velocity distribution, with Portland carbonate, again, showing the most dispersed profiles.We can use the characteristic length L = 100 μm for the beadpack and the estimated values of 139.9 and 327.0 μm for the sandstone and carbonate, respectively, to estimate typical times to diffuse across one pore (if we assume that this is the characteristic length).This time is of the order τ diff = L 2 D m giving 4.54, 8.89, and 48.6 s for the beadpack, sandstone, and carbonate, respectively.For reference the Péclet numbers (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) are 41.4,65.5, and 187.3, respectively, for the beadpack, sandstone, and carbonate.For comparison, the time to transit a characteristic length by advection is equal to diffusion if the velocity is 2.42 × 10 −2 , 1.53 × 10 −2 , and 5.34 × 10 −3 of the average speed for the three media.Regions with a much lower velocity than this are effectively stagnant, since diffusive transport is much faster over the characteristic length.
For the beadpack, the time to diffuse a characteristic length is longer than the time taken to reach asymptotic, Fickian, behavior, indicating a lack of spatial correlation: To sample the velocity distribution, a solute particle only needs to diffuse around one pore length at most, as presented in Fig. 3(a).For the sandstone and carbonate, however, the stagnant peak persists for longer, giving more anomalous behavior.
If we have a system that is macroscopically homogeneousthat is, displays statistically the same pore structure over much greater lengths than those studied-then, eventually, an asymptotic limit is reached in all cases.The dispersion coefficient is related to the spread of velocities and is greatest for the carbonate, less for the sandstone, and lowest for sand or beadpacks [30].However, real systems display distinct geological structures at several length scales and so this asymptotic limit may never be reached: Before the solute has time to sample the local flow field, it is transported into a region with a different pore structure with its own distinctive distribution of velocities.This explains how, at the field scale (100-1000 s m) and for transport times of days or months, anomalous behavior can still be observed ( [13], and references therein; [15,40,41]).The structure at all scales affects the transport, and a firm basis at the pore scale is necessary to account properly for large-scale behavior [42].

B. Comparison with experiment
We compare the computed propagators with those measured in NMR experiments by Scheven et al. [12] in Fig. 6.The rock types, flow speeds, diffusion coefficient, and displacement times are the same in both the experiments and simulations: There are no adjustable or tunable parameters in our model.
There is a very good agreement between the experiments and simulation.As discussed above we see a rapid approach to Gaussian behavior in the beadpack [Fig.6  For the beadpack and Bentheimer sandstone, there appears to be a slight shift in the peak between experiment and the model.This could be due to systematic undersampling of slow-flow regions in the experiments that led to an apparent higher average speed [12].The agreement is particularly good for the most heterogeneous sample, Portland carbonate.The implication here is that capturing the wide spread of velocities is the key to capturing transport correctly.These results suggest that given a good image of the pore space it is possible to make a priori predictions of transport.

IV. CONCLUSIONS
We have described and explained different generic transport behavior in three porous media with increasing pore-scale complexity.In porous media with a relatively homogeneous pore space, giving a narrow spread in local velocity distribution, the exemplar of which is represented by a disordered close packing of equally sized spheres, Fickian transport is rapidly attained through molecular diffusion at the pore scale.On the other hand, with a wider spread in local velocity distribution such as in Bentheimer sandstone, a stagnant peak is seen at early times, with a pronounced mobile peak only emerging later.A greater fraction of the pore space is effectively stagnant and to sample the velocity distribution fully requires diffusive exchange between mobile and immobile regions.Finally and most importantly, in porous media with a very widespread distribution of local velocities, the exemplar of which is Portland carbonate, there is a persistent immobile concentration peak at longer times, with a smaller secondary mobile peak, leading to a highly anomalous behavior.The propagators obtained by the model are in excellent agreement with experimental measurements by Scheven et al. [12].The work demonstrates that pore-scale modeling can provide reliable predictions of core-scale (millimeter to centimeter scale) transport.

FIG. 1 .
FIG. 1. Two-dimensional cross sections of the segmented image of the beadpack (a) and the three-dimensional gray scale images for Bentheimer sandstone (b) and Portland carbonate (c).

FIG. 2 .
FIG. 2. (Color) Image geometries of the beadpack (a), Bentheimer sandstone (b), and Portland carbonate (c) shown with the pore volume represented by gray color; normalized pressure fields with a unit pressure difference across the model for beadpack (d), Bentheimer sandstone (e), and Portland carbonate (f); normalized flow fields, where the ratios of the magnitude of u at the voxel centers divided by the average flow speed u av are shown using cones (too small to be seen individually) that are colored using a logarithmic scale spanning from 5 to 500 for the beadpack (g), Bentheimer sandstone (h), and Portland carbonate (i).In the images green and red colors indicate high values, while blue color indicates low values.

FIG. 4 .
FIG. 4. Probability density function of the velocity distributions for the beadpack, Bentheimer sandstone, and Portland carbonate presented on (a) semilogarithmic axes, and (b) doubly logarithmic axes.The solid line is the distribution for a single cylindrical tube, representing the homogeneous limit.

FIG. 5 .
FIG. 5. (Color) Probability of molecular displacement P (ζ ) in the images of (a) beadpack, (b) Bentheimer sandstone, and (c) Portland carbonate as a function of displacement ζ for the set of times t = 0.106, 0.2, 0.45, 1, and 2 s.The coordinates are rescaled by the expected nominal mean displacement ζ 0 = u av t in the direction of flow.The average velocities are u av = 0.91, 1.03, and 1.3 mm/sfor the beadpack, Bentheimer sandstone, and Portland carbonate, respectively.These correspond to the flow speeds in the experiments[12] with which we compare our results in Sec.III B.
(a)], an immobile stagnant peak that gradually disappears and a large mobile peak developing in Bentheimer sandstone [Fig.6(b)], and

FIG. 6 .
FIG. 6. Probability of molecular displacements P (ζ ) (solid lines) as a function of displacement ζ for the set of times t = 0.106, 0.2, 0.45, 1, and 2 s compared to the propagator obtained in the NMR experiments (dashed lines) by Scheven et al. [12] at the same observation times: (a) the beadpack, (b) Bentheimer sandstone, and (c) Portland carbonate.The coordinates are rescaled by the nominal mean displacement ζ 0 = u av t, u av = 0.91 mm/s for the beadpack, u av = 1.03mm/s for Bentheimer sandstone, and u av = 1.3 mm/s for Portland carbonate, as in the experiments.