Pore scale modelling of calcite cement dissolution in a reservoir sandstone matrix

Long-term evolution of permeability and tortuosity due to porosity changes evoked by reactivity of aqueous solutions is of paramount importance for predicting water-rock interaction. This challenge is best tackled by introducing reactive transport modelling on the pore-scale, where the modeling domain is a high-resolution tomographic image of the porous media. We suggest to use a voxel based Navier-Stokes-Brinkman solver in a finite volume formulation coupled to the thermodynamic equilibrium code PhreeqC. High-performance parallelized computations using this coupled numerical reactive transport solver are performed directly on the voxel grid of the segmented micro-CT scans. Retreat of the calcite cement in a sandstone matrix due to dissolution reactions can be directly visualized by digital rock physics experiments.


Introduction
In this contribution, we suggest coupling the highly robust GeoDict code for fluid flow and colloid transport (www.math2market.com) with the capable geochemical equilibrium code PhreeqC (wwwbrr.cr.usgs.gov) to solve reactive transport problems directly on the voxelized scale geometry of porous media as represented by segmented 3D micro-computer tomographic (-CT) scans. A challenge is yet the extreme computational demand by which this approach is strongly limited in space and time. The effectiveness of this coupling in terms of high performance computing and realization in a massively parallelized manner is therefore essential to offer promising potential. Our workflow is aimed at exploring the dynamically embedded applicability of solution transport simulation on basis of a Navier-Stokes-Brinkman (NSB) algorithm. The common features of the NSB solvers are (i) direct simulation on 3D computer tomography images, (ii) their usage of the finite volume approach based on an adaptive grid structure and (iii) their favourable discretization of the symmetric boundary conditions. Similar to the lattice Boltzmann equation modeling approach [1], computations using the NSB solvers are performed directly on the voxel grid of the segmented -CT scans. This solver has been favourably evaluated against other solvers in a variety of benchmark tests [e.g., 2,3]. Clearly, accuracy of the numerical results are highly dependent on the accuracy of the segmentation, which necessitates using advanced methods [4,5]. This case study is aimed at exploring the capability of our new approach in predicting and visualizing long-term (decades) dynamics of calcite cement dissolution in a sandstone reservoir matrix on high spatial and temporal resolution.

Reservoir sandstone and fluid sample
The microstructure and geochemical system for the digital rock physics experiment was chosen to represent reactive flow in a typical geological reservoir. The calcite dissolution simulations were performed in a -CT generated pore geometry of a reservoir sandstone sample. The sample plug was derived from the Rotliegend geological formation in Altmark region of Germany. The -CT scans were measured at the new P05 imaging beamline of PETRA III synchrotron facility (DESY, Hamburg, Germany). With a circumference of 2.3 km, PETRA III is the biggest and most brilliant storage ring light source in the world. The spatial resolution of the reconstructed geometry was 1.28 m with an unprecedented brilliant phase contrast. For the simulation runs, a region-of-interest (ROI) of 400 3 voxels was chosen yielding in a cuboid sample volume of just about 0.5 mm edge length. Phase segmentation using the materials property simulation code GeoDict revealed beside the dominating quartz content an initial open porosity fraction of 0.033 and a calcite cement volume fraction of 0.168. Composition of the highly saline formation fluid with an initial pH 6.3 is compiled in Table 1. For the simulation runs, the fluid composition was reduced to the main components of interest and equilibrated with calcite at atmospheric pressure but elevated temperature of 60 o C by employing PhreeqC script modules. The inflow fluid was adjusted at a lower pH 5.5 by adding respective amount of hydrochloric acid to induce calcite dissolution reactions.

Geodict code
The FlowDict module of the GeoDict code is an iterative finite volume NSB solver of the Stokes equations. The solver computes the permeability tensor, as well as pressure and fluid velocity distributions, on 3D -CT images. It uses an adaptive instead of a regular grid to reduce the number of grid cells significantly. The basis of this adaptive gridding is called the Left-Identity-Right (LIR) tree that is used for spatial partitioning of the 3D images at reduced computational demand [6]. The pores are represented as differently sized rectangular cuboid voxel cells. In addition, solid regions do not occupy computational memory, which can thus reduce the required RAM size. The pore space is coarsened in areas where the pressure and velocity do not vary much, while keeping the original resolution near the solid surfaces and in regions where pressure and velocity vary rapidly. The GeoDict solver is thus very memory efficient and relatively fast in computations of highly porous microstructures. Clearly, the speed at which convergence is achieved depends also on the complexity and inhomogeneity of the pore structure. Nonetheless, the runtime per iteration is relatively short due to the small number of cells.
Limitations of that code include challenges in modeling slip-boundary conditions, albeit which are not used in this case study. The code allows also for colloid transport using the AddiDict module. New particles may be created at the inflow plane at each time step implementing a probability algorithm. Upon breakthrough at the opposite outflow plane, particles are removed from further simulation runs. In AddiDict, the size of particles can be chosen for calculating physical interactions with the structure while applying the Stokes-Einstein equation. However, an alternative approach is chosen setting the size to a diameter→0, which simplifies the applied equation by neglecting the physical interactions requiring only the diffusion coefficient as input. Particle movement is simulated by Brownian motion at time intervals much smaller (typically in the order of 5 ms) than the fluid transport time steps to ensure a limitation of advective motion to several voxels, thus particles do not get stuck and flow continuously. Number of particles thus moved through a modeling domain can reach dozens of millions.

Transport and geochemical model
At each time step of 1 s chosen, the flow field is computed using symmetric boundaries in all directions by employing the LIR solver. Computations are based on results from previous time steps implying runtime improvements. A pressure gradient of 50 Pa was applied resulting in an initial mean fluid velocity at the permeable pore space of 24.8 m/s. With the applied diffusion coefficient of 1⋅10 -9 m 2 /s, the grid Péclet number is indicating predominantly a diffusion controlled transport. However, along the main flow paths and in particular at pore throats, advection may become the dominating transport regime.
For the geochemical process simulation using our newly developed ReacDict module, we suggest to introduce a simple trick exploiting the capabilities of the Lagrangian colloid transport capability of the AddiDict module. This is that the inflow reactive fluid is not just replacing the initial fluid residing in each pore voxel. Instead, the inflow fluid is transported in form of virtual blobs. They are given a virtual mass corresponding to the fluid density, which may be the same or different of the pore fluid density. At the geochemical calculation intervals of 10 ms, the inflowing fluid blobs of composition given in Table 1 were equilibrated with the pore fluid thus decreasing the respective chemical compound composition in the pore voxels. In particular, local undersaturation with respect to the calcite phase may occur leading to a kinetically constrained local solid phase dissolution. The kinetic reaction calculations were based on the default dataset in the PhreeqC database. In the acidic pH range used, only the first rate term k1•a(H + )•(1-10 1-(2/3•SI) ), with k1 = 7.3•10 -6 and the calcite saturation index SI, is of relevance. The specific surface area was parameterized based on the locally available calcite mass and a global factor of 20 for calcite surface roughness. The resulting dissolved calcite mass was converted to voxel fraction considering a calcite density of 2.715 g/cm 3 . In addition, it was multiplied by an alteration factor of 2,592,000, thus one second of kinetic reaction time step corresponds to 30 days of real reaction time. This parameter was adjusted in optimizing between significant runtime improvement and the accuracy of the dissolution model at the sub-voxel (i.e., solute blob) scale. Per geochemical calculation interval, the kinetic reaction calculations thus result in calcite dissolution of less than 0.14 vol.-% in each reacting voxel containing calcite volume faction. Upon dissolution, solid voxels may become porous and finally form new pore space depending from the mass of dissolved material over time. The overall workflow is depicted in Fig. 1.

Results and Discussion
Due to the acidic inflow solution, calcite dissolution was initiated at sub-voxel scale leading successively to opening and extending of flow paths, which could be followed in time and space at high resolution. With a simulation run time of about 45 min a real kinetic reaction time of 222 years could thus be simulated and visualized in form of a time lapse movie. Notably, 13,350 CPU hours were necessary for such a simulation run with up to 1.6•10 7 solution blobs moving in the 400 3 voxel space, while performing 5.9•10 9 single PhreeqC calls resulting in 1.2•10 8 individual calcite dissolution events at the sub-voxel scale. The thus induced continuous porosity increase effects the permeability development as shown in Fig. 2, but the relationship was intriguingly non-linear. Dissolution reactions in diffusioncontrolled sub-ROI regions of low local Péclet numbers had a negligible effect on the overall permeability evolution of the sandstone microstructure. While the porosity changes are directly linked to the reaction rate, the permeability change depends also on the local flow field. Jumps in the curve can be related to opening of flow paths during dissolution. The development of the Damköhler number with reaction time was calculated from the ratio of the calcite reaction rate to the velocity in flow direction (Fig. 2). The non-linear curve starts with a plateau reflecting rather equivalent increases of the reaction rate and the velocity. After 400 s simulation time (treaction = 33 a), the Damköhler curve start to decrease monotonously up to about three orders of magnitude at 2,700 s flow time (treaction = 222 a), when the accessible calcite is nearly completely dissolved. Related to the reactive surface area determined from the 3D CT scan and multiplied by a factor of 20 to account for surface roughness, the average calcite reaction rate is determined at 1.5±0.6•10 -5 mol/m 2 s. The rate curve at Fig. 2 shows a rather small scatter induced by local transport effects. However, the reaction rate is higher by a factor of 15±6 than an experimentally obtained reaction rate for micritic calcite at 0.98±0.04•10 -6 mol/m 2 s [7], which is obviously due to the surface roughness factor of 20 taken into account. The discrepancy in reaction rates due to deviating concepts in reactive surface area assessment is known to be a critical point in experimental result interpretation [7] and in comparison with numerical results [8]. Nonetheless, the presented pore scale modeling workflow offers promising potential and capabilities for the resimulation of experimental data or even the prediction of experimental results. Time-resolved hydromechanical transport parameters and reaction rates can thus be obtained directly from datasets of dynamic -CT studies. We acknowledge DESY (Hamburg, Germany), a member of the Helmholtz Association HGF, for the provision of experimental facilities. Parts of this research were carried out at PETRA III and we would like to thank the beamline scientists in charge Fabian Wilde and Felix Beckmann for assistance in acquiring the -CT scans of the rock sample at the P05 imaging beamline.