ogst
Home All issues Volume 76 (2021) Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles, 76 (2021) 51 Full HTML
Open Access
Issue
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 76, 2021
Article Number 51
Number of page(s) 8
DOI https://doi.org/10.2516/ogst/2021032
Published online 30 June 2021

© C. Soulaine et al., published by IFP Energies nouvelles, 2021

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

The accurate description of hydrodynamic dispersion according to the flow conditions is one of the long-standing challenges in hydrogeology and subsurface engineering [1, 2]. Briefly, in porous formations, the transport of chemical species is modelled using advection–dispersion equations instead of advection–diffusion equations [3]. In its simplest form, the equation governing the evolution of concentration, CA, reads,

ϕ∂CA∂t+∇∙(UCA)=∇.(ϕD*∙∇CA),$$ \phi \frac{\mathrm{\partial }{C}_A}{\mathrm{\partial }t}+\nabla \bullet \left({U}{C}_A\right)=\nabla.\left(\phi {\mathrm{D}}^{\mathrm{*}}\bullet \nabla {C}_A\right), $$(1)where ϕ is the medium porosity (dimensionless), U is Darcy’s velocity (in m/s) and D* is the so-called dispersion tensor (in m2/s). Generally speaking, the spreading of a solute is not governed only by the molecular diffusion DA (in m2/s) but also by the microstructure and the local velocity field. On the one hand, the tortuous feature of the porous structure characterizes by the medium tortuosity (dimensionless), τ, tends to slow down the spreading. On the other hand, hydrodynamic dispersion stretches a solute band in the flow direction during its transport. The dispersion tensor is a combination of these two effects and writes,

D*=DAτ-1(I+Dhydro).$$ {\mathrm{D}}^{\mathrm{*}}={D}_A{{\tau }}^{-1}\left(\mathrm{I}+{\mathrm{D}}_{\mathrm{hydro}}\right). $$(2)

Empirical laws including Archies’s equation attempt to relate the tortuosity to the porosity using power-law. The determination of hydrodynamic dispersion (dimensionless), Dhydro, is not trivial because it depends on the heterogeneity of the local velocity profile and its magnitude with respect to molecular diffusion. The competition between transport by advection and transport by diffusion is commonly characterized using the Péclet number, Pe, defined as the ratio between characteristic diffusion and advection times. For low Péclet numbers (Pe ≪ 1), the transport is dominated by diffusion and hydrodynamic dispersion is negligible, D* = DAτ−1. For high Péclet numbers (Pe ≫ 1), species are mainly transported by advection and hydrodynamic dispersion is dominant, D* = DAτ−1Dhydro. In the latter case, the dispersion tensor can be orders of magnitude larger than the molecular diffusion. Although linear dispersion models – i.e. the dispersion tensor varies linearly with the Péclet number – are often used for describing contaminant transport [4], they do not rely on strong theoretical background. Indeed, Taylor [5]and Aris [6] have demonstrated that for straight tubes, D*=DA(1+Pe248)$ {D}^{\mathrm{*}}={D}_A\left(1+\frac{\mathrm{P}{\mathrm{e}}^2}{48}\right)$. Empirical measurements suggest that the dispersion can be modelled using a power-law of the Péclet number with an exponent ranging from 0.5 to 2 [7].

Digital Rock Physics (DRP) – one of the main technological breakthroughs in geosciences that appears in the last decades for characterizing rock properties [8–11] – can also be used to determine the dispersion tensor [12–14]. This emerging technology results from the combination of high-resolution imaging and efficient algorithms for solving the governing equations of physical processes at the pore-scale [15]. DRP allows the assessment of effective properties including porosity, geometric and reactive surface area [16–18], absolute permeability [19], Forchheimer correction [20], and relative permeabilities and capillary pressure curves [21] in a non-destructive manner. Current capabilities involving High-Performance Computing enable pore-scale simulations on samples that reach the size of a Representative Elementary Volume (REV) of the porous rock [19, 22–24]. The evaluation of the dispersion tensor is usually performed in two steps. First, the local velocity field within the porosity is obtained by solving Navier–Stokes equations. Then, the transport is solved and upscaled to obtain the values of the dispersion tensor.

Various approaches have been used to compute the solute transport dispersion. Basically, they can be classified in two categories: the Lagrangian and Eulerian approaches. In the former, random walkers are used to simulate particle transport along the streamlines of the velocity field [12, 14, 24, 25]. It includes the Continuous Time Random Walk (CTRW) for which the asymptotic dispersion coefficient is computed from the rate of change of the variance of the particle positions [1, 2, 26]. In Eulerian approaches, the solute transport is described by an advection–diffusion equation solved on an Eulerian grid and the resulting breakthrough curve is fitted against the Ogata-Banks analytical solution of a 1D advection–dispersion equation [27] to obtain the dispersion value [28, 29]. An alternative grid-based approach relying on the volume averaging theory solves directly a set of partial differential equations – the closure problem – whose volume-averaged solution gives the value of the dispersion tensor for a given flow condition [30, 31]. Unlike other approaches, the latter solves a steady-state problem and directly provides the asymptotic values of the dispersion tensor. Comparison between the methods has highlighted that the dispersion coefficients estimated using volume averaging and particle tracking agree closely with each other [32].

In this paper, we propose a DRP technique for calculating the dispersion tensor using 3D images of the microstructure. The approach is based on the work of Carbonell and Whitaker [30] who propose a general formulation for the dispersion tensor using the method of volume averaging and the concept of closure problems [33]. We implemented the closure problem for computing the dispersion tensor within the open-source platform OpenFOAM®. The resulting program called dispersionEvaluationFoam is available online (https://github.com/csoulain/dispersionEvaluationFoam) and benefits from all OpenFOAM features including High-Performance Computing and unstructured grids to process 3D images of rock samples [19].

The paper is organized as follows. In Section 2, we describe the acquisition of the 3D images of a sandpack, the computation of the velocity profile, and the numerical method to determine the dispersion tensor. In Section 3, we first verify the accuracy of our numerical framework by simulating cases for which analytical solution exists, and we then compute the dispersion tensor of sandpack imaged using X-ray microtomography. Finally, we close with a summary and conclusions.

2 Material and methods

In this section, we present the methodology to compute the dispersion tensor including the acquisition and gridding of microtomography images, the computation of the velocity profile, and the evaluation of the dispersion tensor. Our approach relies heavily on open-source packages including the Computational Fluid Dynamics platform OpenFOAM (https://www.openfoam.org).

2.1 Sandpack sample

The sandpack sample is made of grains of different shape and size (Fig. 1a). The grain-size distribution measured with a laser microparticle size analyzer ranges from 0.1 to 0.8 mm with a mean diameter of d50 = 0.3 mm (Fig. 1b). Then mean density of sandpack is ρs = 2650 kg m−3.

thumbnail Fig. 1

(a) Scanning electron microscope image of the sandpack. (b) Grain size distribution measured with microparticle size analyzer. (c) X-ray microtomography of the sandpack reconstructed with the software VG Studio.

2.2 Image acquisition, segmentation, and gridding

Three-dimensional images of the sandpack are obtained using X-ray microtomography. It is a non-invasive imaging technique that allows for constructing a 3D image of an object using a set of two-dimensional radiographs of the X-ray attenuation properties of the object materials. The sandpack sample was imaged at the Earth Sciences Institute of Orléans using a Nanotom microtomograph 180NF Phoenix X-Ray. We used an accelerating voltage of 180 kV, a filament current of 170 nA, and an operating voltage of 120 V. The stack of 2D images is processed to construct a volume with the VGStudio MAX 1.2 (Volume Graphics). The 3D volume consists of 1000 × 1000 × 1800 voxels with a resolution of 1.2 μm/voxel. In the raw images, each voxel value corresponds to a measure of linear-absorption coefficients. A nonlinear digital filtering (median filter 3 × 3 × 3) is applied to reduce the noise before a segmentation algorithm is used to binarize the images into a void and solid phases (see Fig. 1c).

In this paper, we solve flow and transport processes within the pore-space using grid-based approaches. The computational grid is obtained from the segmented images in two consecutive steps. First, the surface of the solid grains is extracted into a Stereolithography Interface Format (STL) using ImageJ (https://www.imagej.net) and the 3Dviewer plugin. Then, the mesh is generated based on the STL using snappyHexMesh, the OpenFOAM automatic gridder that maximizes the number of hexahedral cells. For large grids, this procedure can be performed in parallel using High-Performance Computing.

2.3 Computation of the velocity profile

The velocity profile within the pore-space is obtained by solving Stokes equation for cyclic domains. The mass and momentum equations read [34],

∇∙v=0,$$ \nabla \bullet {v}=0, $$(3)

0=-∇p̃+μ∇2v+ΔPLe0,$$ 0=-\nabla \mathop{p}\limits^\tilde+\mu {\nabla }^2{v}+\frac{\Delta P}{L}{{e}}_0, $$(4)where v is the velocity field in the pore-space, p̃$ \mathop{p}\limits^\tilde$ is the deviation to the mean pressure field, μ is the fluid viscosity, and ΔPLe0$ \frac{\Delta P}{L}{{e}}_0$ is a body source term that describes the mean pressure gradient within the computational domain.

No-slip conditions,

v=0,$$ {v}=\mathbf{0}, $$(5)are used at the solid walls. Cyclic conditions applied at the inlet and outlet edges of the computational domain. It means that p̃$ \mathop{p}\limits^\tilde$ and v are the same at the inlet and outlet patches. The mean flow rate is adjusted by changing ΔPL$ \frac{\Delta P}{L}$. The mean flow orientation is governed by the unit vector e0.

We use OpenFOAM to solve equations (3)–(4). The flow equations for cyclic domains are discretized using a finite-volume method and solved using the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) pressure–velocity coupling procedure developed by Patankar [35]. We adapt the existing simpleFoam solver to integrate the mean pressure drop in the momentum equation. The under-relaxation factors are set to 0.2 for pressure and 0.9 for velocity. The pressure equation – obtained by combining the mass and momentum equations (3)–(4) – is solved using a Geometric Algebraic Multi-Grid (GAMG) linear solver. The simulations are considered fully converged when the residuals are below 10−8.

2.4 Evaluation of the dispersion tensor

Using the method of volume averaging [33], Carbonell and Whitaker [30] have demonstrated that the dispersion tensor can be described by,

D*=DA(I+〈∇B〉f)-〈ṽB〉f,$$ {\mathrm{D}}^{\mathrm{*}}={D}_A\left(I+\left\langle \nabla {B}{\right\rangle^f\right)-\left\langle \stackrel{\tilde }{{v}}{B}{\right\rangle^f, $$(6)where the operator 〈·〉f is the volume average of a quantity over the volume occupied by the fluid within the computational domain. The vector B, also called closure variable is solution of the boundary value problem,

ṽ+v∙∇B=∇∙(DA∇B) in the void space,$$ \stackrel{\tilde }{{v}}+{v}\bullet \nabla {B}=\nabla \bullet \left({D}_A\nabla {B}\right)\enspace \mathrm{in}\mathrm{\enspace }\mathrm{the}\mathrm{\enspace }\mathrm{void}\mathrm{\enspace }\mathrm{space}, $$(7)

n∙∇B=-n at the fluid/solid boundary,$$ {n}\bullet \nabla {B}=-{n}\mathrm{\enspace }\mathrm{at}\mathrm{\enspace }\mathrm{the}\mathrm{\enspace }\mathrm{fluid}/\mathrm{solid}\mathrm{\enspace }\mathrm{boundary}, $$(8)where v is the velocity profile at the pore-scale computing in Section 2.3, ṽ$ \stackrel{\tilde }{{v}}$ = v − 〈v〉f is the deviation to the mean velocity, and n is the normal vector at the fluid–solid interface. Equations (7)–(8) are supplemented by cyclic conditions at the edge of the computational domain and by the additional constrain that the average of B is zero,

〈B〉f=0.$$ \left\langle {B}{\right\rangle^f=0. $$(9)

Hence, once the velocity profile in a REV is known, B can be computed numerically solving equations (7)–(9). The dispersion tensor for a given flow condition is obtained using equation (6). The variation of the dispersion tensor along with the Péclet number is determined by solving iteratively the closure problem with different values of DA.

The B-field is solution of equations (7)–(8) up to a constant that is determined by the zero average constrain, equation (9). Although crucial, it is not straightforward to solve the boundary value problem along with equation (9). dispersionEvaluationFoam – our OpenFOAM® implementation of Carbonell and Whitaker’s closure problem – solves a slightly different set of equations to satisfy the zero average constrain. Instead, a common practice consists in solving the partial differential equation along with the associated boundary conditions for a variable B0 that has a fixed reference value on a reference cell (e.g. B0 = 0 for the cell labelled 0). The closure variable is then obtained using,

B=B0-〈B0〉f.$$ {B}={{B}}_0-\left\langle {{B}}_0{\right\rangle^f. $$(10)

In dispersionEvaluationFoam, the advection–diffusion equation (7), solving for the closure variable is discretized using the finite-volume method. The interfacial condition, equation (8), is implemented in a dedicated boundary condition. To facilitate the convergence of the calculation, we use under-relaxation. dispersionEvaluationFoam allows a fast and direct approach to compute the dispersion tensor.

3 Results and discussion

In this section, we use dispersionEvaluationFoam to compute the dispersion tensor of various pore geometries. First, we verify that we recover Taylor–Aris law for the hydrodynamic dispersion within a single straight cylinder. Then, we compute the dispersion tensor of a sandpack using 3D images obtained by microtomography.

3.1 Verification: Taylor–Aris dispersion

The objective of this part is to verify the robustness of dispersionEvaluationFoam on cases for which reference solutions exist. The investigation of flow and transport in straight cylindrical channels is convenient not only because analytical solutions for the momentum and mass balance equations exist but also because they constitute the base-element of pore-network models – i.e. a simplification of the pore microstructure as a network of cylinders [36]. For an infinite cylinder of radius R in the x-axis, the velocity profile is parabolic, v=2〈vx〉f(1-(rR)2)ex$ {v}=2\left\langle {v}_x{\right\rangle^f\left(1-{\left(\frac{r}{R}\right)}^2\right){{e}}_x$, where 〈vx〉f is the velocity averaged over a cross-section, and the longitudinal dispersion is exactly D*=DA(1+Pe248)$ {D}^{\mathrm{*}}={D}_A\left(1+\frac{\mathrm{P}{\mathrm{e}}^2}{48}\right)$ known as the Taylor–Aris dispersion where the Péclet number is Pe=〈vx〉fRDA$ \mathrm{Pe}=\frac\left\langle {v}_x{\right\rangle^fR}{{D}_A}$ [5, 6]. We verify that our OpenFOAM code for evaluating dispersion tensor predicts Taylor–Aris law.

Considering that a cross-section is a Representative Elementary Volume of an infinite tube, the equations introduced in Section 2.4 also lead to Taylor–Aris dispersion. Indeed, in cylindrical coordinates (er, eθ, ex), the closure variable simplifies to B = Bx(r)ex (because of invariance by rotation and translation along the tube axis) and equations (7)–(9) write,

ṽx=DA1r(∂∂r(r∂Bx∂r)),$$ {\mathop{v}\limits^\tilde}_x={D}_A\frac{1}{r}\left(\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r\frac{\mathrm{\partial }{B}_x}{\mathrm{\partial }r}\right)\right), $$(11)

∂Bx∂r|r=R=0,$$ {\left.\frac{\mathrm{\partial }{B}_x}{\mathrm{\partial }r}\right|}_{r=R}=0, $$(12)

∫0RrBx(r)dr=0.$$ {\int }_0^R r{B}_x(r)\mathrm{d}r=0. $$(13)

The integration of equations (11)–(13) gives,

Bx(r)=〈vx〉R24DA((rR)2-12(rR)4-13).$$ {B}_x(r)=\frac\left\langle {v}_x\right\rangle{R}^2}{4{D}_A}\left({\left(\frac{r}{R}\right)}^2-\frac{1}{2}{\left(\frac{r}{R}\right)}^4-\frac{1}{3}\right). $$(14)

Finally, using equations (6), the dispersion tensor is

Dxx*DA=1-2R2DA∫0RrṽxBx(r)dr=1+〈vx〉2R248DA2,$$ \frac{{D}_{{xx}}^{\mathrm{*}}}{{D}_A}=1-\frac{2}{{R}^2{D}_A}{\int }_0^R r{\mathop{v}\limits^\tilde}_x{B}_x(r)\mathrm{d}r=1+\frac\left\langle {v}_x{\right\rangle^2{R}^2}{48{D}_A^2}, $$(15)which corresponds to Taylor–Aris formula.

We verify that dispersionEvaluationFoam converges to equation (14) and to Taylor–Aris dispersion, equation (15). We consider a three-dimensional cylinder, 2 mm long and 260 μm. Cyclic boundary conditions are applied on both ends of the cylinder so that the computational domain is a Representative Elementary Volume of an infinite tube. The cylinder is gridded with 13 600 hexahedral cells (see Fig. 2a). We use a full 3D mesh instead of an axisymmetric grid to illustrate the capabilities of dispersionEvaluationFoam to handle 3D unstructured grids.

thumbnail Fig. 2

The longitudinal dispersion in a straight capillary tube is obtained using dispersionEvaluationFoam: (a) the cylinder is meshed with an unstructured grid and contains cyclic boundaries, (b) the velocity profile is parabolic.

First, we compute the velocity profile solving Stokes equations according to the methodology described in Section 2.3. Simulation parameters are μ = 10−6 m2/s, ΔPL$ \frac{\Delta P}{L}$ = 106 Pa/m and e0 = ex. As expected, we obtain a parabolic profile (see Fig. 2b).

Then, equations (7)–(9) are solved using dispersionEvaluationFoam along with the velocity profile freshly computed. We performed multiple simulations using diffusion coefficients, DA, ranging from 10−7 to 10−13 m2/s which corresponds to Péclet numbers between 10−2 and 104. In Figure 3, we plot the profile of the closure variable Bx along a channel radius. The numerical results are in very good agreement with those obtained analytically with equation (14). In Figure 4, we observe a very good match between the values of the longitudinal dispersion obtained numerically for different Péclet numbers and Taylor–Aris formula.

thumbnail Fig. 3

Radial profile of the field B obtained using both equation (14) and dispersionEvaluationFoam is the case of an infinite cylinder.

thumbnail Fig. 4

Longitudinal dispersion in a straight capillary tube obtained analytically using Taylor–Aris dispersion model and numerically using dispersionEvaluationFoam. There is a perfect agreement between the theory and the simulations.

This verification case highlights the capabilities of dispersionEvaluationFoam to handle 3D unstructured grids. It gives full confidence to use dispersionEvaluationFoam to compute dispersion tensor on more complex pore geometries including microstructure provided by 3D images.

3.2 Longitudinal dispersion of a sandpack

In this section, we assess the asymptotic longitudinal hydrodynamic dispersion of a sandpack using 3D images obtained with microtomography imaging.

The initial set of images contain 1800 × 1000 × 1000 voxels, and the resolution is 1.2 μm3/voxel. The stack of images is segmented into fluid and solid phases (Fig. 1). The porosity of the sandpack sample, obtained by counting the number of voxels containing fluid – or void –, is ϕ = 0.36. We verified that this value of the porosity is stable when the averaging volume was smaller than the sample size. Although this indicates that the 3D images can be considered as a REV with respect to porosity, it does not necessarily mean that sample volume is a REV of the dispersion coefficient because the latter is usually much larger than the REV of porosity. Therefore, the hereby method assesses the dispersion coefficient of the given sample.

The numerical assessment of the sample permeability and hydrodynamic dispersion requires to grid the void space. First, the fluid–mineral surface is extracted (in STL format) from the segmented images using ImageJ (https://www.imagej.net). Then, the surface mesh is used with OpenFOAM automatic gridder, snappyHexMesh, to generate the computational grid. To be able to apply cyclic boundary conditions, we add two manifolds at the inlet and outlet of the domain. The solid surface and the edges of the domain are set as impermeable walls. The final mesh contains 2.7 × 106 cells. The flow and dispersion calculations are conducted in parallel using 8 CPUs.

The velocity profile in the sandpack sample is obtained by solving equations (3)–(4) using simpleFoam, an OpenFOAM solver for steady-state Navier–Stokes equations. Simulation parameters are μ = 1.5 × 10−5 m2/s, ΔPL$ \frac{\Delta P}{L}$ = 103 Pa/m. The pressure field, velocity vectors and magnitude plotted in Figure 5 illustrate the tortuous and non-uniform distribution of the flow pathways in the sandpack porosity. The permeability is obtained by averaging the velocity profile and using Darcy’s law:

K=μ〈v〉(ΔPL)-1∙e0,$$ K=\mu \left\langle {v}\right\rangle{\left(\frac{\Delta P}{L}\right)}^{-1}\bullet {{e}}_0, $$(16)where 〈v〉 is the volume averaged velocity profile. We obtain K = 3.1 × 10−11 m2 which is typical for sands [3].

thumbnail Fig. 5

Pressure and velocity profiles obtained numerically in a sandpack. The sandpack 1.2 mm × 1.2 mm × 2.16 mm three-dimensional images have been obtained using microtomography. The flow simulation is carried out with OpenFOAM. The velocity vectors and magnitude illustrate the tortuous and non-uniform distribution of the flow pathways in the pore-space. The velocity vectors are colored with the pressure field. The profiles are normalized by the maximum and minimum values and range between 0 (blue) and 1 (red).

The longitudinal hydrodynamic dispersion is calculated using dispersionEvaluationFoam for Péclet number ranging from 10−4 to 104. This is achieved by solving the mathematical problem introduced in Section 2.4 for different values of the diffusivity DA. The resulting dispersion values are plotted in Figure 6. The values obtained numerically vary with the Péclet number as a power-law:

Dxx*DA=0.52(1+35Pe1.25).$$ \frac{{D}_{{xx}}^{\mathrm{*}}}{{D}_A}=0.52\left(1+35\mathrm{P}{\mathrm{e}}^{1.25}\right). $$(17)

thumbnail Fig. 6

Longitudinal dispersion in a sandpack obtained numerically using dispersionEvaluationFoam. The asymptotic hydrodynamic dispersion varies with the Péclet number according to a power-law. The exponent is lower than Taylor–Aris dispersion but larger than a mere linear dispersion model. Below Pe ≈ 1, hydrodynamic dispersion is negligible and the value of the plateau corresponds to the sandpack tortuosity.

For Péclet number below 1, the diffusion is the dominant transport mechanism and hydrodynamic dispersion is negligible. Therefore, D* is independent of the Péclet number for Pe < 1. The value of the plateau corresponds to the sandpack tortuosity, τ ≈ 1/0.52 ≈ 1.9. For advection-dominated transport, Pe > 1, and the hydrodynamic dispersion scales with Pe1.25. The exponent 1.25 is in good agreement with the experimental and numerical values obtained in the literature for sandpacks [7, 12, 37]. The exponent is lower than Taylor–Aris dispersion but larger than a mere linear dispersion model that is often used in transport simulators.

4 Conclusion

We computed the longitudinal dispersion of a sandpack using pore-scale simulations. Three-dimensional images of rock samples were obtained using microtomography imaging, segmented into fluid and solid phase, and gridded. The determination of the hydrodynamic dispersion tensor relies on the closure problem proposed by Carbonell and Whitaker [30]. We implemented their approach into a modern and open-source platform. Our numerical development, called dispersionEvaluationFoam, is based on the scientific library OpenFOAM for solving partial differential equations using the finite-volume method. dispersionEvaluationFoam is available publicly on our GitHub depository (https://github.com/csoulain/dispersionEvaluationFoam). Due to the versatility of OpenFOAM, our model handles 3D simulations, massively parallel calculations, and unstructured grids.

We verify the predictive aspect of dispersionEvaluationFoam by assessing the hydrodynamic dispersion in an infinite cylinder. The results are in very good agreement with Taylor–Aris theory. Then, we evaluate the tortuosity and the longitudinal dispersion of a sandpack. We found a dependency in hydrodynamic conditions that scales Pe1.25 in agreement with data published in the literature.

The methodology proposed in this paper is generic and can be applied to any kind of rock. Importantly, the tools and the package we developed to perform our simulations are open-source, and we expect that our developments will be an important contribution to the flow and transport in porous media community.

Acknowledgments

The research leading to these results has received funding from the French Agency for Research (Agence Nationale de la Recherche, ANR) through the labex Voltaire ANR-10-LABX-100-01 and the grant FraMatI ANR-19-CE05-0002. It has also received financial support from the CNRS through the MITI interdisciplinary programs and the Region Centre-Val de Loire (Contribution 201500103985). The authors benefitted from the use of the cluster at the Centre de Calcul Scientifique en région Centre-Val de Loire.

References

  • Dentz M., Carrera J. (2007) Mixing and spreading in stratified flow, Phys. Fluids 19, 017107. ISSN 1070-6631. https://doi.org/10.1063/1.2427089. [CrossRef] [Google Scholar]
  • Berkowitz B., Cortis A., Dentz M., Scher H. (2006) Modeling non-fickian transport in geological formations as a continuous time random walk, Rev. Geophys. 440, 2, 1–49. ISSN 1944-9208. https://doi.org/10.1029/2005RG000178. [Google Scholar]
  • Bear J. (1972) Dynamics of fluids in porous media, Elsevier, New York. [Google Scholar]
  • Parkhurst D.L., Appelo C.A.J. (2013) Description of input and examples for PHREEQC version 3 – a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations, Vol. 6–A43, U.S. Department of the Interior, U.S. Geological Survey Techniques and Methods, Reston, VA. [Google Scholar]
  • Taylor G.I. (1953) Dispersion of soluble matter in solvent flowing slowly through a tube, Proc. Roy. Soc. A. 219, 186–203. [Google Scholar]
  • Aris R. (1956) On the dispersion of a solute in a fluid flowing through a tube, Proc. Roy. Soc. A. 235, 67–77. [Google Scholar]
  • Bijeljic B., Blunt M.J. (2006) Pore-scale modeling and continuous time random walk analysis of dispersion in porous media, Water Resour. Res. 42, W01202, 1–5. ISSN 0043-1397. https://doi.org/10.1029/2005wr004578. [Google Scholar]
  • Spanne P., Thovert J.F., Jacquin C.J., Lindquist W.B., Jones K.W., Adler P.M. (1994) Synchrotron computed microtomography of porous media: Topology and transports, Phys. Rev. Lett. 730, 14, 2001. [Google Scholar]
  • Andrä H., Combaret N., Dvorkin J., Glatt E., Han J., Kabel M., Keehm Y., Krzikalla F., Lee M., Madonna C., Marsh M., Mukerji T., Saenger E.H., Sain R., Saxena N., Ricker S., Wiegmann A., Zhan X. (2013) Digital rock physics benchmarks part I: Imaging and segmentation, Comput. Geosci. 500, 25–32. [Google Scholar]
  • Andrä H., Combaret N., Dvorkin J., Glatt E., Han J., Kabel M., Keehm Y., Krzikalla F., Lee M., Madonna C., Marsh M., Mukerji T., Saenger E.H., Sain R., Saxena N., Ricker S., Wiegmann A., Zhan X. (2013) Digital rock physics benchmarks part II: Computing effective properties, Comput. Geosci. 500, 33–43. Benchmark problems, datasets and methodologies for the computational geosciences. [Google Scholar]
  • Soulaine C., Maes J., Roman S. (2021) Computational microfluidics for geosciences, Front. Water 3, 1–11. https://doi.org/10.3389/frwa.2021.643714. [Google Scholar]
  • Bijeljic B., Mostaghimi P., Blunt M.J. (2011) Signature of non-fickian solute transport in complex heterogeneous porous media, Phys. Rev. Lett. 107, 204502. [Google Scholar]
  • Nakashima Y., Watanabe Y. (2002) Estimate of transport properties of porous media by microfocus X-ray computed tomography and random walk simulation, Water Resour. Res. 38, 8–1–8–12. ISSN 0043-1397. https://doi.org/10.1029/2001wr000937. [Google Scholar]
  • Kang P.K., Anna P., Nunes J.P., Bijeljic B., Blunt M.J., Juanes R. (2014) Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media, Geophys. Res. Lett. 410, 17, 6184–6190. [Google Scholar]
  • Blunt M.J., Bijeljic B., Dong H., Gharbi O., Iglauer S., Mostaghimi P., Paluszny A., Pentland C. (2013) Pore-scale imaging and modelling, Adv. Water Resour. 51, 197–216. [Google Scholar]
  • Noiriel C., Soulaine C. (2021) Pore-scale imaging and modelling of reactive flow in evolving porous media: Tracking the dynamics of the fluid-rock interface, Trans. Porous Media. https://doi.org/10.1007/s11242-021-01613-2. [Google Scholar]
  • Molins S., Trebotich D., Steefel C.I., Shen C. (2012) An investigation of the effect of pore scale flow on average geochemical reaction rates using direct numerical simulation, Water Resour. Res. 48, W03527, 1–11. ISSN 0043-1397. https://doi.org/10.1029/2011wr011404. [Google Scholar]
  • Soulaine C., Roman S., Kovscek A., Tchelepi H.A. (2017) Mineral dissolution and wormholing from a pore-scale perspective, J. Fluid Mech. 827, 457–483. [Google Scholar]
  • Guibert R., Nazarova M., Horgue P., Hamon G., Creux P., Debenest G. (2015) Computational permeability determination from pore-scale imaging: Sample size, mesh and method sensitivities, Trans. Porous Media 1070, 3, 641–656. [Google Scholar]
  • Soulaine C., Quintard M. (2014) On the use of a Darcy – Forchheimer like model for a macro-scale description of turbulence in porous media and its application to structured packings, Int. J. Heat Mass Transf. 740, 88–100. [Google Scholar]
  • Raeini A.Q., Blunt M.J., Bijeljic B. (2014) Direct simulations of two-phase flow on micro-ct images of porous media and upscaling of pore-scale forces, Adv. Water Resour. 740, 116–126. ISSN 0309-1708. [Google Scholar]
  • Soulaine C., Gjetvaj F., Garing C., Roman S., Russian A., Gouze P., Tchelepi H. (2016) The impact of sub-resolution porosity of X-ray microtomography images on the permeability, Trans. Porous Media 1130, 1, 227–243. [Google Scholar]
  • Scheibe T.D., Perkins W.A., Richmond M.C., McKinley M.I., Romero-Gomez P.D.J., Oostrom M., Wietsma T.W., Serkowski J.A., Zachara J.M. (2015) Pore-scale and multiscale numerical simulation of flow and transport in a laboratory-scale column, Water Resour. Res. 510, 2, 1023–1035. [Google Scholar]
  • Bijeljic B., Raeini A., Mostaghimi P., Blunt M.J. (2013) Predictions of non-fickian solute transport in different classes of porous media using direct simulation on pore-scale images, Phys. Rev. E 87, 013011. [Google Scholar]
  • Noetinger B., Roubinet D., Russian A., Le Borgne T., Delay F., Dentz M., de Dreuzy J.-R., Gouze P. (2016) Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale, Trans. Porous Media 115, 345–385. ISSN 0169-3913. https://doi.org/10.1007/s11242-016-0693-z. [Google Scholar]
  • De Anna P., Le Borgne T., Dentz M., Tartakovsky A.M., Bolster D., Davy P. (2013) Flow intermittency, dispersion, and correlated continuous time random walks in porous media, Phys. Rev. Lett. 1100, 18, 184502. [Google Scholar]
  • Ogata A., Banks R.B. (1961) A solution of the differential equation of longitudinal dispersion in porous media. Number 411,A in Geological Survey professional Paper, United States Department of the Interior. US Government Printing Office, Washington, DC. [Google Scholar]
  • Ortega-Ramírez M.P., Oxarango L. (2021) Effect of X-ray μ-ct resolution on the computation of permeability and dispersion coefficient for granular soils, Trans. Porous Media 137, 307–326. ISSN 0169-3913. https://doi.org/10.1007/s11242-021-01557-7. [Google Scholar]
  • Zaretskiy Y., Geiger S., Sorbie K., Förster M. (2010) Efficient flow and transport simulations in reconstructed 3d pore geometries, Adv. Water Resour. 330, 12, 1508–1516. ISSN 0309-1708. https://doi.org/10.1016/j.advwatres.2010.08.008. URL https://www.sciencedirect.com/science/article/pii/S0309170810001521. [Google Scholar]
  • Carbonell R.G., Whitaker S. (1983) Dispersion in pulsed systems – II: Theoretical developments for passive dispersion in porous media, Chem. Eng. Sci. 380, 11, 1795–1802. [Google Scholar]
  • Wood B.D. (2007) Inertial effects in dispersion in porous media, Water Resour. Res. 430, W12S16, 1–16. ISSN 0043-1397. https://doi.org/10.1029/2006wr005790. [Google Scholar]
  • Richmond M.C., Perkins W.A., Scheibe T.D., Lambert A., Wood B.D. (2013) Flow and axial dispersion in a sinusoidal-walled tube: Effects of inertial and unsteady flows, Adv. Water Resour. 62, 215–226. ISSN 0309-1708. https://doi.org/10.1016/j.advwatres.2013.06.014. [Google Scholar]
  • Whitaker S. (1999) The method of volume averaging, theory and applications of transport in porous media, Kluwer Academic, Dorderecht. [Google Scholar]
  • Roman S., Soulaine C., Abu AlSaud M., Kovscek A., Tchelepi H. (2016) Particle velocimetry analysis of immiscible two-phase flow in micromodels, Adv. Water Resour. 95, 199–211. [Google Scholar]
  • Patankar S.V. (1980) Numerical heat transfer and fluid flow, Taylor & Francis, Washington, DC. [Google Scholar]
  • Fatt I. (1956) The network model of porous media, Trans. AIME 207, 01, 144–181. [Google Scholar]
  • Pfannkuch H.O. (1963) Contribution à l’étude des deplacements de fluides miscibles dans un milieu poreux, Rev. Inst. Fr. Pet. 18, 215–270. [Google Scholar]

All Figures

thumbnail Fig. 1

(a) Scanning electron microscope image of the sandpack. (b) Grain size distribution measured with microparticle size analyzer. (c) X-ray microtomography of the sandpack reconstructed with the software VG Studio.
In the text
thumbnail Fig. 2

The longitudinal dispersion in a straight capillary tube is obtained using dispersionEvaluationFoam: (a) the cylinder is meshed with an unstructured grid and contains cyclic boundaries, (b) the velocity profile is parabolic.
In the text
thumbnail Fig. 3

Radial profile of the field B obtained using both equation (14) and dispersionEvaluationFoam is the case of an infinite cylinder.
In the text
thumbnail Fig. 4

Longitudinal dispersion in a straight capillary tube obtained analytically using Taylor–Aris dispersion model and numerically using dispersionEvaluationFoam. There is a perfect agreement between the theory and the simulations.
In the text
thumbnail Fig. 5

Pressure and velocity profiles obtained numerically in a sandpack. The sandpack 1.2 mm × 1.2 mm × 2.16 mm three-dimensional images have been obtained using microtomography. The flow simulation is carried out with OpenFOAM. The velocity vectors and magnitude illustrate the tortuous and non-uniform distribution of the flow pathways in the pore-space. The velocity vectors are colored with the pressure field. The profiles are normalized by the maximum and minimum values and range between 0 (blue) and 1 (red).
In the text
thumbnail Fig. 6

Longitudinal dispersion in a sandpack obtained numerically using dispersionEvaluationFoam. The asymptotic hydrodynamic dispersion varies with the Péclet number according to a power-law. The exponent is lower than Taylor–Aris dispersion but larger than a mere linear dispersion model. Below Pe ≈ 1, hydrodynamic dispersion is negligible and the value of the plateau corresponds to the sandpack tortuosity.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.