Regular Article
Digital Rock Physics: computation of hydrodynamic dispersion
^{1}
Earth Sciences Institute of Orléans, CNRSUniversité d’OrléansBRGM, 45100 Orléans, France
^{2}
Laboratoire GéHCO, Campus Grandmont, Université de Tours, 37200 Tours, France
^{*} Corresponding author: cyprien.soulaine@cnrsorleans.fr
Received:
10
February
2021
Accepted:
26
May
2021
Hydrodynamic dispersion is a crucial mechanism for modelling contaminant transport in subsurface engineering and water resources management whose determination remains challenging. We use Digital Rock Physics (DRP) to evaluate the longitudinal dispersion of a sandpack. From a threedimensional image of a porous sample obtained with Xray microtomography, we use the method of volume averaging to assess the longitudinal dispersion. Our numerical implementation is opensource and relies on a modern scientific platform that allows for large computational domains and HighPerformance Computing. We verify the robustness of our model using cases for which reference solutions exist and we show that the longitudinal dispersion of a sandpack scales as a power law of the Péclet number. The assessment methodology is generic and applies to any kind of rock samples.
© C. Soulaine et al., published by IFP Energies nouvelles, 2021
This 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 longstanding 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, C_{A}, reads,
$$\varphi \frac{\mathrm{\partial}{C}_{A}}{\mathrm{\partial}t}+\nabla \bullet \left(\mathit{U}{C}_{A}\right)=\nabla .\left(\varphi {\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 socalled dispersion tensor (in m^{2}/s). Generally speaking, the spreading of a solute is not governed only by the molecular diffusion D_{A} (in m^{2}/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,
$${\mathrm{D}}^{\mathrm{*}}={D}_{A}{\mathit{\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 powerlaw. The determination of hydrodynamic dispersion (dimensionless), D_{hydro}, 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^{*} = D_{A}τ^{−1}. For high Péclet numbers (Pe ≫ 1), species are mainly transported by advection and hydrodynamic dispersion is dominant, D^{*} = D_{A}τ^{−1}D_{hydro}. 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}^{\mathrm{*}}={D}_{A}\left(1+\frac{\mathrm{P}{\mathrm{e}}^{2}}{48}\right)$. Empirical measurements suggest that the dispersion can be modelled using a powerlaw 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 highresolution imaging and efficient algorithms for solving the governing equations of physical processes at the porescale [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 nondestructive manner. Current capabilities involving HighPerformance Computing enable porescale 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 OgataBanks analytical solution of a 1D advection–dispersion equation [27] to obtain the dispersion value [28, 29]. An alternative gridbased approach relying on the volume averaging theory solves directly a set of partial differential equations – the closure problem – whose volumeaveraged solution gives the value of the dispersion tensor for a given flow condition [30, 31]. Unlike other approaches, the latter solves a steadystate 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 opensource platform OpenFOAM^{®}. The resulting program called dispersionEvaluationFoam is available online (https://github.com/csoulain/dispersionEvaluationFoam) and benefits from all OpenFOAM features including HighPerformance 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 Xray 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 opensource 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 grainsize distribution measured with a laser microparticle size analyzer ranges from 0.1 to 0.8 mm with a mean diameter of d_{50} = 0.3 mm (Fig. 1b). Then mean density of sandpack is ρ_{s} = 2650 kg m^{−3}.
Fig. 1 (a) Scanning electron microscope image of the sandpack. (b) Grain size distribution measured with microparticle size analyzer. (c) Xray microtomography of the sandpack reconstructed with the software VG Studio. 
2.2 Image acquisition, segmentation, and gridding
Threedimensional images of the sandpack are obtained using Xray microtomography. It is a noninvasive imaging technique that allows for constructing a 3D image of an object using a set of twodimensional radiographs of the Xray 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 XRay. 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 linearabsorption 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 porespace using gridbased 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 HighPerformance Computing.
2.3 Computation of the velocity profile
The velocity profile within the porespace is obtained by solving Stokes equation for cyclic domains. The mass and momentum equations read [34],
$$\nabla \bullet \mathit{v}=0,$$(3)
$$0=\nabla \stackrel{\u0303}{p}+\mu {\nabla}^{2}\mathit{v}+\frac{\Delta P}{L}{\mathit{e}}_{0},$$(4)where v is the velocity field in the porespace, $\stackrel{\u0303}{p}$ is the deviation to the mean pressure field, μ is the fluid viscosity, and $\frac{\Delta P}{L}{\mathit{e}}_{0}$ is a body source term that describes the mean pressure gradient within the computational domain.
Noslip conditions,
$$\mathit{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 $\stackrel{\u0303}{p}$ and v are the same at the inlet and outlet patches. The mean flow rate is adjusted by changing $\frac{\Delta P}{L}$. The mean flow orientation is governed by the unit vector e_{0}.
We use OpenFOAM to solve equations (3)–(4). The flow equations for cyclic domains are discretized using a finitevolume method and solved using the SIMPLE (SemiImplicit 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 underrelaxation 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 MultiGrid (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,
$${\mathrm{D}}^{\mathrm{*}}={D}_{A}\left(I+\langle \nabla \mathit{B}{\rangle}^{f}\right)\langle \stackrel{\u0303}{\mathit{v}}\mathit{B}{\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,
$$\stackrel{\u0303}{\mathit{v}}+\mathit{v}\bullet \nabla \mathit{B}=\nabla \bullet \left({D}_{A}\nabla \mathit{B}\right)\mathrm{in}\mathrm{}\mathrm{the}\mathrm{}\mathrm{void}\mathrm{}\mathrm{space},$$(7)
$$\mathit{n}\bullet \nabla \mathit{B}=\mathit{n}\mathrm{}\mathrm{at}\mathrm{}\mathrm{the}\mathrm{}\mathrm{fluid}/\mathrm{solid}\mathrm{}\mathrm{boundary},$$(8)where v is the velocity profile at the porescale computing in Section 2.3, $\stackrel{\u0303}{\mathit{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,
$$\langle \mathit{B}{\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 D_{A}.
The Bfield 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 B_{0} that has a fixed reference value on a reference cell (e.g. B_{0} = 0 for the cell labelled 0). The closure variable is then obtained using,
$$\mathit{B}={\mathit{B}}_{0}\langle {\mathit{B}}_{0}{\rangle}^{f}.$$(10)
In dispersionEvaluationFoam, the advection–diffusion equation (7), solving for the closure variable is discretized using the finitevolume method. The interfacial condition, equation (8), is implemented in a dedicated boundary condition. To facilitate the convergence of the calculation, we use underrelaxation. 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 baseelement of porenetwork models – i.e. a simplification of the pore microstructure as a network of cylinders [36]. For an infinite cylinder of radius R in the xaxis, the velocity profile is parabolic, $\mathit{v}=2\langle {v}_{x}{\rangle}^{f}\left(1{\left(\frac{r}{R}\right)}^{2}\right){\mathit{e}}_{x}$, where 〈v_{x}〉^{f} is the velocity averaged over a crosssection, and the longitudinal dispersion is exactly ${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 $\mathrm{Pe}=\frac{\langle {v}_{x}{\rangle}^{f}R}{{D}_{A}}$ [5, 6]. We verify that our OpenFOAM code for evaluating dispersion tensor predicts Taylor–Aris law.
Considering that a crosssection 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 (e_{r}, e_{θ}, e_{x}), the closure variable simplifies to B = B_{x}(r)e_{x} (because of invariance by rotation and translation along the tube axis) and equations (7)–(9) write,
$${\stackrel{\u0303}{v}}_{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)
$${\frac{\mathrm{\partial}{B}_{x}}{\mathrm{\partial}r}}_{r=R}=0,$$(12)
$${\int}_{0}^{R}\mathrm{}r{B}_{x}\left(r\right)\mathrm{d}r=0.$$(13)
The integration of equations (11)–(13) gives,
$${B}_{x}\left(r\right)=\frac{\langle {v}_{x}\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
$$\frac{{D}_{\mathrm{xx}}^{\mathrm{*}}}{{D}_{A}}=1\frac{2}{{R}^{2}{D}_{A}}{\int}_{0}^{R}\mathrm{}r{\stackrel{\u0303}{v}}_{x}{B}_{x}\left(r\right)\mathrm{d}r=1+\frac{\langle {v}_{x}{\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 threedimensional 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.
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} m^{2}/s, $\frac{\Delta P}{L}$ = 10^{6} Pa/m and e_{0} = e_{x}. 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, D_{A}, ranging from 10^{−7} to 10^{−13} m^{2}/s which corresponds to Péclet numbers between 10^{−2} and 10^{4}. In Figure 3, we plot the profile of the closure variable B_{x} 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.
Fig. 3 Radial profile of the field B obtained using both equation (14) and dispersionEvaluationFoam is the case of an infinite cylinder. 
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 μm^{3}/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 × 10^{6} 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 steadystate Navier–Stokes equations. Simulation parameters are μ = 1.5 × 10^{−5} m^{2}/s, $\frac{\Delta P}{L}$ = 10^{3} Pa/m. The pressure field, velocity vectors and magnitude plotted in Figure 5 illustrate the tortuous and nonuniform distribution of the flow pathways in the sandpack porosity. The permeability is obtained by averaging the velocity profile and using Darcy’s law:
$$K=\mu \langle \mathit{v}\rangle {\left(\frac{\Delta P}{L}\right)}^{1}\bullet {\mathit{e}}_{0},$$(16)where 〈v〉 is the volume averaged velocity profile. We obtain K = 3.1 × 10^{−11} m^{2} which is typical for sands [3].
Fig. 5 Pressure and velocity profiles obtained numerically in a sandpack. The sandpack 1.2 mm × 1.2 mm × 2.16 mm threedimensional images have been obtained using microtomography. The flow simulation is carried out with OpenFOAM. The velocity vectors and magnitude illustrate the tortuous and nonuniform distribution of the flow pathways in the porespace. 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 10^{4}. This is achieved by solving the mathematical problem introduced in Section 2.4 for different values of the diffusivity D_{A}. The resulting dispersion values are plotted in Figure 6. The values obtained numerically vary with the Péclet number as a powerlaw:
$$\frac{{D}_{\mathrm{xx}}^{\mathrm{*}}}{{D}_{A}}=0.52\left(1+35\mathrm{P}{\mathrm{e}}^{1.25}\right).$$(17)
Fig. 6 Longitudinal dispersion in a sandpack obtained numerically using dispersionEvaluationFoam. The asymptotic hydrodynamic dispersion varies with the Péclet number according to a powerlaw. 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 advectiondominated transport, Pe > 1, and the hydrodynamic dispersion scales with Pe^{1.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 porescale simulations. Threedimensional 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 opensource platform. Our numerical development, called dispersionEvaluationFoam, is based on the scientific library OpenFOAM for solving partial differential equations using the finitevolume 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 Pe^{1.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 opensource, 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 ANR10LABX10001 and the grant FraMatI ANR19CE050002. It has also received financial support from the CNRS through the MITI interdisciplinary programs and the Region CentreVal de Loire (Contribution 201500103985). The authors benefitted from the use of the cluster at the Centre de Calcul Scientifique en région CentreVal de Loire.
References
 Dentz M., Carrera J. (2007) Mixing and spreading in stratified flow, Phys. Fluids 19, 017107. ISSN 10706631. https://doi.org/10.1063/1.2427089. [CrossRef] [Google Scholar]
 Berkowitz B., Cortis A., Dentz M., Scher H. (2006) Modeling nonfickian transport in geological formations as a continuous time random walk, Rev. Geophys. 440, 2, 1–49. ISSN 19449208. 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, batchreaction, onedimensional 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) Porescale modeling and continuous time random walk analysis of dispersion in porous media, Water Resour. Res. 42, W01202, 1–5. ISSN 00431397. 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 nonfickian 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 Xray computed tomography and random walk simulation, Water Resour. Res. 38, 8–1–8–12. ISSN 00431397. https://doi.org/10.1029/2001wr000937. [Google Scholar]
 Kang P.K., Anna P., Nunes J.P., Bijeljic B., Blunt M.J., Juanes R. (2014) Porescale intermittent velocity structure underpinning anomalous transport through 3d 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) Porescale imaging and modelling, Adv. Water Resour. 51, 197–216. [Google Scholar]
 Noiriel C., Soulaine C. (2021) Porescale imaging and modelling of reactive flow in evolving porous media: Tracking the dynamics of the fluidrock interface, Trans. Porous Media. https://doi.org/10.1007/s11242021016132. [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 00431397. https://doi.org/10.1029/2011wr011404. [Google Scholar]
 Soulaine C., Roman S., Kovscek A., Tchelepi H.A. (2017) Mineral dissolution and wormholing from a porescale 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 porescale 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 macroscale 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 twophase flow on microct images of porous media and upscaling of porescale forces, Adv. Water Resour. 740, 116–126. ISSN 03091708. [Google Scholar]
 Soulaine C., Gjetvaj F., Garing C., Roman S., Russian A., Gouze P., Tchelepi H. (2016) The impact of subresolution porosity of Xray 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., RomeroGomez P.D.J., Oostrom M., Wietsma T.W., Serkowski J.A., Zachara J.M. (2015) Porescale and multiscale numerical simulation of flow and transport in a laboratoryscale column, Water Resour. Res. 510, 2, 1023–1035. [Google Scholar]
 Bijeljic B., Raeini A., Mostaghimi P., Blunt M.J. (2013) Predictions of nonfickian solute transport in different classes of porous media using direct simulation on porescale 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 01693913. https://doi.org/10.1007/s112420160693z. [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]
 OrtegaRamírez M.P., Oxarango L. (2021) Effect of Xray μct resolution on the computation of permeability and dispersion coefficient for granular soils, Trans. Porous Media 137, 307–326. ISSN 01693913. https://doi.org/10.1007/s11242021015577. [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 03091708. 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 00431397. 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 sinusoidalwalled tube: Effects of inertial and unsteady flows, Adv. Water Resour. 62, 215–226. ISSN 03091708. 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 twophase 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
Fig. 1 (a) Scanning electron microscope image of the sandpack. (b) Grain size distribution measured with microparticle size analyzer. (c) Xray microtomography of the sandpack reconstructed with the software VG Studio. 

In the text 
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 
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 
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 
Fig. 5 Pressure and velocity profiles obtained numerically in a sandpack. The sandpack 1.2 mm × 1.2 mm × 2.16 mm threedimensional images have been obtained using microtomography. The flow simulation is carried out with OpenFOAM. The velocity vectors and magnitude illustrate the tortuous and nonuniform distribution of the flow pathways in the porespace. 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 
Fig. 6 Longitudinal dispersion in a sandpack obtained numerically using dispersionEvaluationFoam. The asymptotic hydrodynamic dispersion varies with the Péclet number according to a powerlaw. 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 