Digital Rock Physics: computation of hydrodynamic dispersion

. 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 three-dimensional image of a porous sample obtained with X-ray microtomography, we use the method of volume averaging to assess the longitudinal dispersion. Our numerical implementation is open-source and relies on a modern scienti ﬁ c platform that allows for large computational domains and High-Performance 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.


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, C A , reads, where / is the medium porosity (dimensionless), U is Darcy's velocity (in m/s) and D * is the so-called 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), s, 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, Empirical laws including Archies's equation attempt to relate the tortuosity to the porosity using power-law.
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 s À1 . For high Péclet numbers (Pe ) 1), species are mainly transported by advection and hydrodynamic dispersion is dominant, D * = D A s À1 D hydro .
In the latter case, the dispersion tensor can be orders of magnitude larger than the molecular diffusion. Although linear dispersion modelsi.e. the dispersion tensor varies linearly with the Péclet numberare 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, 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][9][10][11] can also be used to determine the dispersion tensor [12][13][14]. This emerging technology results from the combination of highresolution 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][17][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][23][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 equationsthe closure problemwhose 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/dispersion EvaluationFoam) 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.

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).

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 d 50 = 0.3 mm (Fig. 1b). Then mean density of sandpack is q s = 2650 kg m À3 .

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 lm/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 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 Open-FOAM automatic gridder that maximizes the number of hexahedral cells. For large grids, this procedure can be performed in parallel using High-Performance Computing.

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], where v is the velocity field in the pore-space,p is the deviation to the mean pressure field, l is the fluid viscosity, and ÁP L e 0 is a body source term that describes the mean pressure gradient within the computational domain.
are used at the solid walls. Cyclic conditions applied at the inlet and outlet edges of the computational domain. It means thatp and v are the same at the inlet and outlet patches. The mean flow rate is adjusted by changing Á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 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 equationobtained 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 .

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, where the operator hÁi 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, n Á rB ¼ Àn at the fluid=solid boundary; where v is the velocity profile at the pore-scale computing in Section 2.3,ṽ = v À hvi 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, 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 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). dispersionEvaluationFoamour OpenFOAM Ò implementation of Carbonell and Whitaker's closure problemsolves 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, In dispersionEvaluationFoam, the advectiondiffusion 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. dispersionEval-uationFoam allows a fast and direct approach to compute the dispersion tensor.

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.

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 modelsi.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 where hv x i f is the velocity averaged over a cross-section, and the longitudinal dispersion is exactly D Ã ¼ D A 1 þ Pe 2 48 À Á known as the Taylor-Aris dispersion where the Péclet number is Pe ¼ hvx i f R 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 (e r , e h , 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, The integration of equations (11)- (13) gives, Finally, using equations (6), the dispersion tensor is 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 lm. 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.
First, we compute the velocity profile solving Stokes equations according to the methodology described in Section 2.3. Simulation parameters are l = 10 À6 m 2 /s, Á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 disper-sionEvaluationFoam 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. This verification case highlights the capabilities of dispersionEvaluationFoam to handle 3D unstructured grids. It gives full confidence to use dispersion EvaluationFoam to compute dispersion tensor on more complex pore geometries including microstructure provided by 3D images.

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 lm 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 fluidor 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 steady-state Navier-Stokes equations. Simulation parameters are l = 1.5 Â 10 À5 m 2 /s, ÁP L = 10 3 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: where hvi is the volume averaged velocity profile. We obtain K = 3.1 Â 10 À11 m 2 which is typical for sands [3]. 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 power-law: 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, s % 1/0.52 % 1.9. For advection-dominated 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. Fig. 3. Radial profile of the field B obtained using both equation (14) and dispersionEvaluationFoam is the case of an infinite cylinder.

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. dispersionE-valuationFoam is available publicly on our GitHub depository (https://github.com/csoulain/dispersionEvalu-ationFoam). Due to the versatility of OpenFOAM, our model handles 3D simulations, massively parallel calculations, and unstructured grids. We verify the predictive aspect of dispersionEval-uationFoam 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 open-source, and we expect that our developments will be an important contribution to the flow and transport in porous media community.   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.