Simulation of dissolution in porous media in three dimensions with lattice Boltzmann, finite-volume, and surface-rescaling methods

We present a pore-scale dissolution model for the simulation of reactive transport in complex porous media such as those encountered in carbon-storage injection processes. We couple a lattice Boltzmann model for flow calculation with a finite-volume method for solving chemical transport equations, and allow the computational grid to change as mineral surfaces are dissolved according to first-order reaction kinetics. We appraise this scheme for use with high Péclet number flows in three-dimensional geometries and show how the popular first-order convection scheme is affected by severe numerical diffusion when grid Péclet numbers exceed unity, and confirm that this can be overcome relatively easily by using a second-order method in conjunction with a flux-limiter function. We then propose a surface rescaling method which uses parabolic elements to counteract errors in surface area exposed by the Cartesian grid and avoid the use of more complex embedded surface methods when surface reaction kinetics are incorporated. Finally, we compute dissolution in an image of a real porous limestone rock sample injected with HCl for different Péclet numbers and obtain dissolution patterns in concordance with theory and experimental observation. A low injection flow rate was shown to lead to erosion of the pore space concentrated at the face of the rock, whereas a high flow rate leads to wormhole formation.

Figure 1

A discrete surface (blue is solid; white is fluid) on which a smooth surface is interpolated by parabolic sections (dashed orange lines). Scaling factors ${\psi }^{\mathit{i}}\le 1$ associated with the nodes’ exposed faces are computed from projections of the lattice area onto the parabolic curves. The method is applied to all boundary faces; three examples are shown. The finely dashed black line shows a case where marching cubes would decouple the surface from the solid volumes.

Figure 2

Adjoining cells in one dimension indexed relative to cell $i$, with concentration values $C_i$. The central cell is shown with velocities at the left and right faces of $u^{-}$ and $u^{+}$, and scalar fluxes $f_{i-1/2}$ and $f_{i+1/2}$.

Figure 3

The value of the superbee limiter $\mathrm{\Psi }\left(r_i\right)$ as a function of the curvature sensor $r_i$. The limiter takes on values of 0 for $r_i\le 0$ and 2 when $r_i\ge 2$.

Figure 4

(a) Geometry used for validation and numerical diffusion tests. The diameter of the sphere is half the dimension of the channel; (b) the velocity field: orange (light) zones are fast flow and blue (dark) zones are slow; and (c) steady-state concentration field (scale normalized to inlet concentration) at $\mathrm{Pe}=73$, injected from rear.

Figure 5

Concentration of a tracer plotted along characteristic lines of the geometry. Solid lines: second-order flux-limiter transport model; dashed lines: comsol calculation. (N.B. for grayscale versions of this figure: The order of the curves from top to bottom corresponds to that of the legend read from left to right and downwards).

Figure 6

Concentration fields computed in the geometry of dimension $L=200$ and shown in the central slice for different hydrodynamic Péclet numbers. The convection schemes used were (a) first order and (b) second order with superbee limiter. The scale is normalized to the inlet concentration.

Figure 7

Concentration profiles across the central part of the outlet of the geometry for (a) $\mathrm{Pe}=73$, (b) $\mathrm{Pe}=730$ and (c) $\mathrm{Pe}=7300$. (N.B. for grayscale versions of this figure: The order of the curves corresponds to that of the legend except in (c), where second-order 100 cubed lies below second-order 200 cubed).

Figure 8

A grid cell at position ${\mathit{r}}_0$ with solid fraction $\mathbf{\sigma }\left({\mathit{r}}_0\right)$. An exposed cell face is shown with the unit normal vector ${\mathit{n}}_0$, and eight perpendicular neighbors at displacements $\mathit{d}{\mathit{r}}_{\mathit{i}}$. A Cartesian coordinate system is defined on the face so that $\widehat{\mathit{z}}={\mathit{n}}_0$.

Figure 9

Construction of a parabolic surface and the projection of its normal vector $\left(\delta x,\delta y,\delta z\right)$ onto the grid surface face vector ${\mathit{n}}_0$.

Figure 10

Textured surfaces: (a) a discrete plane in two dimensions, (b) a discrete plane in three dimensions, and (c) a smooth cylinder in which boundary nodes are partially filled according to their overlap with the circular radius.

Figure 11

Ketton carbonate geometry and flow field. The image resolution is 15.8 µm and the grid size is ${200}^3$ lattice units.

Figure 12

The concentration of $H^{+}$ in the dissolution of Ketton carbonate at $\mathrm{Pe}=0.8$, injected from left rear at different times: (a) $t=2.3\mathrm{h}$, (b) $t=150\mathrm{h}$.

Figure 13

The concentration of $H^{+}$ in the dissolution of Ketton carbonate at $\mathrm{Pe}=80$ at different times: (a) $t=19min$, (b) $t=2.7\mathrm{h}$, (c) $t=5.4\mathrm{h}$. HCl was injected from left rear. The wormhole is apparent in (c) as the large flow path at the center of the domain.

Figure 14

Cutaways of the Ketton limestone geometry after dissolution at (a) $\mathrm{Pe}=0.8$ at $t=56.4\mathrm{h}$ and (b) $\mathrm{Pe}=80$ at $t=5.4\mathrm{h}$ showing a wormhole through the sample. Injection was from left front. The porosity is 38.5% in both cases.

Figure 15

The porosity of the Ketton limestone sample against time undergoing dissolution at two different Péclet numbers. The dashed region of the $\mathrm{Pe}=80$ curve indicates where the mechanical movement of the rock grains likely needs to be taken into account.

Figure 16

The permeability of the Ketton carbonate sample as the porosity changes during dissolution at $\mathrm{Pe}=80$ (top line) and $\mathrm{Pe}=0.8$.