Pore‐scale simulation of carbonate dissolution in micro‐CT images

We present a particle‐based method to simulate carbonate dissolution at the pore scale directly on the voxels of three‐dimensional micro‐CT images. The flow field is computed on the images by solving the incompressible Navier‐Stokes equations. Rock‐fluid interaction is modeled using a three‐step approach: solute advection, diffusion, and reaction. Advection is simulated with a semianalytical pore‐scale streamline tracing algorithm, diffusion by random walk is superimposed, while the reaction rate is defined by the flux of particles through the pore‐solid interface. We derive a relationship between the local particle flux and the independently measured batch calcite dissolution rate. We validate our method against a dynamic imaging experiment where a Ketton oolite is imaged during CO2‐saturated brine injection at reservoir conditions. The image‐calculated increases in porosity and permeability are predicted accurately, and the spatial distribution of the dissolution front is correctly replicated. The experiments and simulations are performed at a high flow rate, in the uniform dissolution regime – Pe ≫ 1 and PeDa ≪ 1—thus extending the reaction throughout the sample. Transport is advection dominated, and dissolution is limited to regions with significant inflow of solute. We show that the sample‐averaged reaction rate is 1 order of magnitude lower than that measured in batch reactors. This decrease is the result of restrictions imposed on the flux of solute to the solid surface by the heterogeneous flow field, at the millimeter scale.


Introduction
Applications in geological carbon storage [Wilson and Monea, 2004] and enhanced oil recovery [Pizarro and Branco, 2012] have created a surge in research focused on simulating and predicting geochemical reactions between formation fluids and host rocks triggered by the injection of CO 2 . Following injection, the CO 2 combines with reservoir fluids forming an acidic solution that gradually reacts with alkaline minerals in the rock, particularly in carbonate reservoirs. Petrophysical and mechanical alterations resulting from such reactions may have significant implications for the management and monitoring of these reservoirs. In addition, carbonate dissolution is also of interest when modelling diagenetic processes [Whitaker and Smart, 2011] and matrix acidizing in oil reservoirs [Fredd and Fogler , 1998]. measured with micro-CT imaging, and porosity following limestone dissolution. Studying dissolution in a calcareous sandstone, Lamy-Chappuis et al. [2014] reported evidence that relatively small increases in porosity could precipitate steep changes in permeability. Even in homogeneous rocks the reaction-induced increase in porosity can be non-uniform and different dissolution patterns (e. g. wormholes, face and uniform dissolution) may emerge according to the specifics of the flow and reaction regimes, as shown in Daccord et al. [1993], Fredd and Fogler [1998], Golfier et al. [2002] and Luquot and Gouze [2009], for example.
At the pore scale, one alternative to model rock-fluid interactions is to adopt a hybrid continuum-particle approach, in which reactants are represented by particles moving in a flow field. Some of the earliest results were obtained by Sallès et al. [1993] and Békri et al. [1995], who used a finite-difference method to solve for Stokes flow combined with a random walk method to move particles, simulating dissolution and deposition in threedimensional synthetic rocks. Békri et al. [1995] distinguished three main flow and reaction regimes: reaction-, diffusion-and advection-controlled. They also noticed that dissolution was concentrated along the main flow paths in the advection-controlled regime. However, c 2016 American Geophysical Union. All Rights Reserved. they did not compute an average effective reaction rate. Sadhukhan et al. [2012] solved the Navier-Stokes equation with a finite-difference method in a two-dimensional stochastically generated medium. They related the solute concentration to the porosity increase and showed qualitative agreement with experimental results. In the framework of continuous time random walks (CTRW), Edery et al. [2011] applied a particle tracking method to simulate dedolomitization and carbonate precipitation in synthetic two-dimensional systems.
In this paper, we combine and extend the modelling approaches described above. We use a finite-volume method to solve the Navier-Stokes equation for the flow field in the pore space. We then employ a streamline method to capture advective transport semianalytically, Pereira Nunes et al. [2015]. Diffusion is accommodated using a random particle motion. For reaction, as we show later, we relate the particle flux at the voidsolid boundary to the batch reaction rate.
Particle-tracking methods have also been widely used to study homogeneous reactions in porous media, illustrating the effect of pore-space heterogeneity on the mixing of reactants and how it affects the average reaction rate, e.g. Benson and Meerschaert [2008], Edery et al. [2009, Ding et al. [2013], Hansen et al. [2014] andHenri andFernàndez-Garcia [2014]. Paster et al. [2013]  [2009], while Hansen et al. [2014] expanded this treatment to relate an effective radius, within which reaction occurs, to the intrinsic reaction rate. The implementation of particle methods to systems governed by multispecies reactions and/or high-order kinetics is not straightforward and could, potentially, inhibit the application of the method to more complex systems. However, it is worth noting that there are examples applying particle methods to study both multispecies coupled reactions obeying first-order kinetics [Henri and Fernàndez-Garcia, 2014] and also homogeneous reactions described by non-linear kinetics [Tompson and Dougherty, 1992].
Intrinsic reaction rates for well-mixed fluids flowing through pure mineral suspensions or next to a pure, flat sample can be measured in the laboratory, Lasaga [1998]. In this paper the intrinsic rates (i.e. not affect by transport limitations) will be referred to as batch rates and the average rate (also called effective rate) is the rate at which a sample of porous rock reacts with solute dissolved in the fluid flowing through its pores. The direct application of batch rates to geological reservoirs is questionable. Batch reaction rates are often several orders of magnitude higher than the ones reported in field studies and core measurements, e.g. Swoboda-Colberg and Drever [1993], Salehikhoo et al. [2013], Steefel et al. [2005. This could be the result of: differences between reactive surface areas in fresh and weathered minerals and the fact that natural systems are closer to equilibrium than laboratory systems, White and Brantley [2003]; incomplete mixing of solute reactants, Li et al. [2008]; or the presence of physical and chemical heterogeneities in the field, Steefel et al. [2005], Li et al. [2006]. Another important factor contributing to the decrease of field-scale rates is their dependence on fluid residence times, since apparent reaction rates decrease with residence time, Maher [2010]. Studying silica dissolution rates in HCl c 2016 American Geophysical Union. All Rights Reserved. irrigated plots, Swoboda-Colberg and Drever [1993] reported a difference of two orders of magnitude between laboratory and field rates.
It is reasonable to expect that the same constraints (incomplete mixing of solutes, physical and chemical heterogeneities) are also present at the pore (micron) and core (millimetre/centimetre) scales, thus diminishing the observed reaction rates in comparison with the batch rates. This has been explored computationally in artificial media. Meile and Tuncay [2006] investigated the impact of pore-scale heterogeneity in diffusion-reaction processes using computer generated porous media. For the case of calcite dissolution, they reported a wide range of mismatch (up to two orders of magnitude) between true rates and upscaled rates, when arbitrarily changing the effective reactive surface. A decrease in the effective reaction rate was also predicted by Molins et al. [2012], who employed a continuum approach to study the impact of pore-space heterogeneity in a synthetic two-dimensional porous medium, but the effect -a decrease of only 20% -was much smaller than expected from the laboratory and field results quoted above. A porenetwork approach was employed by Li et al. [2007] to illustrate the possible impacts of variations in the distribution of minerals in reaction rates of sandstone reservoirs subjected to CO 2 injection. They found that the size of the clusters of reactive minerals and their position relative to the flow direction may have major effects on the upscaling of the reaction rates. Experiments by Salehikhoo et al. [2013] and Li et al. [2014] studied the influence of the spatial distribution of minerals and flow velocity on magnesite dissolution.
They found that the flow rate directly influences the centimetre-scale dissolution rate. dissolved during the flow of CO 2 -saturated brine. They reported an average reaction rate approximately 14 times smaller than the batch calcite dissolution rate.
It is still not evident whether the average reaction rate can be lower in porous rocks exclusively due to a heterogeneous flow field. Here we will simulate flow and reaction in a realistic three-dimensional system and show that one order of magnitude decrease in the reaction rate can be caused entirely by the limited transport of reactants in the pore space, a consequence of flow heterogeneity. Our results are in accord with experimental measurements. Previous work focused on simulating fluid-rock reactions in relatively simple synthetic porous media and there is still a need to develop predictive pore-scale models capable of reproducing reactive flow experiments in realistic rock samples and explaining the emergence of dissolution patterns and effective reaction rates from first principles. Our work is a contribution in this direction.
An important outcome that will follow from our simulations is that, since chemical heterogeneities are not considered in the model, physical heterogeneities alone are responsible for the observed average dissolution rate. This will lead us to the conclusion that the restrictions imposed by the characteristics of flow -a reduction of both the amount of solute in contact with the solids and the local residence time -have a major impact in the dissolution dynamics, even for fast flows where the reactants are distributed throughout the pore space. The net effect is a reduction of one order of magnitude in the observed dissolution rate in millimetre-scale rock samples.
Micro-CT imaging is a useful tool to characterize and model flow and transport phenomena at the pore scale, see Blunt et al. [2013] and Wildenschild and Sheppard [2013] c 2016 American Geophysical Union. All Rights Reserved.
for recent reviews. Recent developments have made it possible to image multiphase-flow in rocks at reservoir conditions, Andrew et al. [2013].
In this paper we make use of micro-CT images acquired during a reactive transport experiment performed at high pressure and temperature conditions [Menke et al., 2015] to present and validate a pore-scale particle tracking algorithm that simulates carbonate dissolution caused by CO 2 -saturated fluids from first principles. We start by briefly describing the advection and diffusion steps of the algorithm and follow by deriving a relationship between the intrinsic calcite dissolution rate in acidic brines and the local flux of particles into the solid grains. Next we compare the predictions of our simulations with dynamic imaging data, in which a Ketton carbonate was imaged during reactive flow of brine saturated with CO 2 in the core. We follow by showing that our method is able to accurately predict the increase in porosity and permeability and reproduce the dissolution pattern throughout the sample. To conclude, we discuss how the overall dissolution rate in the sample is controlled by flow heterogeneities and show that the average dissolution rate is one order of magnitude smaller than the batch rate. As stated before, our model assumes chemical homogeneity and such reduction is explained on the basis of a transport-controlled flux of solute to the solids at the millimetre scale.

Fluid Flow and Particle Transport
The numerical solution for the transport of solute particles in porous media involves two steps: 1) computing the velocity field in the pore-space by solving the Navier-Stokes equation and 2) transporting the particles in this field. Particle transport is, in turn, divided into two contributions: advection and diffusion. During advection the particles move along streamlines, which do not traverse the solids. Diffusion is modelled as a c 2016 American Geophysical Union. All Rights Reserved. random motion where particles jump from one streamline to another, eventually crossing the pore-solid interface. This process is shown in schematic form in figure 1. Dissolution is controlled by the local flux of particles at the interfaces.

Solution for the Flow Field
We solve the incompressible Navier-Stokes equation directly in the pore-space of segmented micro-CT images, the simulation gridblocks are the image voxels from the micro-CT image: where µ is the fluid viscosity, P is the pressure, ρ the fluid density and V the velocity. The numerical solution of the Navier-Stokes equation is obtained with the single-phase version of the finite volume method presented in Raeini et al. [2012], that uses the OpenFoam library [OpenFoam, 2013], applying constant pressure at inlet and outlet and imposing the no-slip boundary condition V = 0 at the pore walls, as described by Bijeljic et al.
[2013a]. After the velocity field is obtained by solving eq. 1, the permeability (k) of a sample of length L s and cross-sectional area A is readily determined from Darcy's law: where Q is the total flow rate calculated from the velocity field at the sample outlet and ∆P is the imposed pressure drop.
c 2016 American Geophysical Union. All Rights Reserved.

Particle Transport
The evolution of the solute concentration (C) in a time-independent velocity field is governed by the advection-diffusion equation, Sahimi [2011]: In the Lagrangian approach, instead of solving eq. 4 for the concentration, individual particle trajectories are tracked, Tompson and Gelhar [1990], Kinzelbach [1988]. Each particle follows a Brownian motion where its displacements, after a time-step (∆t), are given by the sum of two terms, an advective displacement and a (random) diffusive displacement: is the position of a particle p at a time t.
The solute concentration is recovered by calculating the average number of particles in a reference volume at a given time.

Advection
The advective term depends on the local velocity and the time-step and is usually calculated with an Euler approximation: ∆ X adv p = V ( X)∆t, which requires that the time-step is made very small, particularly in regions of high velocity, Tompson and Gelhar [1990], Sàlles et al. [1993], Békri et al. [1995]. We overcome this limitation using a semianalytical streamline method, described below.
A modified Pollock method to trace streamlines is presented in Pereira Nunes et al.
[2015]. This is a semi-analytical approach that imposes the non-slip boundary condition at the pore-solid walls and avoids numerical stability issues inherent in conventional timestepping methods, thus allowing the time-step to be constrained by diffusion only. Once c 2016 American Geophysical Union. All Rights Reserved.
the solution for the velocity field is obtained (eq. 1), we trace streamlines through the image voxels and solute particles are transported along these streamlines.

Diffusion
The diffusive term in eq. 5 is modelled with a random walk method to simulate particles jumping across streamlines. The three-dimensional mean-free path is defined by where D m is the diffusion coefficient and ∆t is the simulation time-step. At each time-step the diffusive jump is calculated according to: where λ = 24D m ∆t (7) and ξ is a random vector whose components are uniformly distributed between 0 and 1, Benson and Meerschaert [2008]. During the simulations the time-step is fixed such that the maximum diffusive jump in each direction is equal to half of the voxel size.

Calcite Dissolution in CO 2 -saturated Fluids: A Particle Tracking Approach
The kinetics of calcite (CaCO 3 ) dissolution in the presence of CO 2 -saturated brine has been described by Pokrovsky et al. [2005] and by Peng et al. [2015], who used pure crystals and batch reactors to study calcite dissolution under a range of pressures and temperatures. With a rotating disc apparatus, Peng et al. [2015] were able to determine the dissolution rate free of mass transfer effects at the surface of a single crystal at realistic reservoir conditions. When calcite dissolves, Ca 2+ ions are released in the flowing brine, and our model describes dissolution in terms of their concentration in the aqueous phase.
c 2016 American Geophysical Union. All Rights Reserved.
The solute concentration (with units of number of particles per volume) of a reactive tracer in a system controlled by one single chemical species and described by first-order kinetics satisfies the advection-diffusion-reaction equation [Lasaga, 1998]: where C eq is the equilibrium concentration, at constant pH conditions, of the solute in the liquid phase and Dissolution occurs when C < C eq , and precipitation when C > C eq . The term (C eq − C) represents a deviation from equilibrium and we rewrite eq. 8 in terms of a new variable Note the change of the sign of the last term, implying that C * decreases during dissolution.
From eq. 9 it is evident that the equilibrium concentration in terms of the new variable is C * eq = 0. A Lagrangian approach to the problem of reactive transport replaces eq. 9, governing the evolution of the solute concentration, with equations for the particle trajectories (eq. 5), as discussed in section 2.2; see also Verberg and Ladd [2002]. The reaction term of eq. 9 is defined by the flux of particles through the solid interface, section 3.1.
Under high flow rates (i.e. high Péclet number), pH gradients in the pore-space are insignificant (see section 4.1) and we describe the kinetics of calcite dissolution solely as a function of the Ca 2+ concentration in the liquid phase, which is conveniently written in terms of particles using the variable C * introduced above. Hence, the particle concentration represents the Ca 2+ undersaturation, or the deviation from the equilibrium c 2016 American Geophysical Union. All Rights Reserved.
concentration in the liquid phase. As stated before, the fluid is in equilibrium with respect to the solid when C * = 0 and no dissolution occurs, as expected, since the net flux of particles to the solids is zero. When simulating the injection of a fluid out of equilibrium with the rock, it is more natural to describe the reactions in terms of C * rather than C: particles are injected at the inlet, representing the absence of calcium ions and hence the fluid moves from equilibrium and the reaction starts; the density of particles represents the value of C * . The higher the particle concentration, the further from equilibrium the fluid is. Furthermore, associating the particles with the absence of Ca 2+ rather than with Ca 2+ itself requires fewer particles to be tracked, which is more efficient computationally.

Local Flux of Solute Particles and the Dissolution Rate
To translate the reaction rate into the particle framework, we have to establish a relationship between the particle flux and the batch dissolution rate. We propose a model that is applicable under first-order, reaction-controlled kinetics at the solid surface. Under these assumptions, it is possible to derive analytically an expression relating the flux of particles to the laboratory measured rate in terms of a priori characteristics of the system, as shown below.
The model assumes first-order kinetics, meaning that reaction is proportional to the local concentration of particles in the vicinity of the interfaces. The particle flux at the pore-solid walls is, in turn, determined by diffusion. The flux is calculated by counting the number of particles that cross from the surrounding pore voxels to a specific solid voxel.
Our goal is to represent conditions under which the fluid is constantly undersaturated with Ca 2+ , produced during calcite dissolution, thus implying that the flux of particles is from the pore to the solid. Once a particle crosses the solid-fluid interface it reacts (it c 2016 American Geophysical Union. All Rights Reserved. is removed from the simulation) and the number of hits which that particular solid voxel has received is updated. After a given number of hits (emulating the reaction rate), the solid voxel is dissolved into a pore voxel.
To establish a relationship between particle flux and the batch reaction rate we must consider the number of hits to dissolve, which is simply the cumulative flux of particles to the solid, and relate it to the diffusion length, the particle concentration, the reactive surface area and the number of moles of calcite in each solid voxel. We assume that the solid voxels are pure calcite, a very good approximation for the Ketton oolite, composed of 99% calcite, Andrew et al. [2014].
We start by considering the case of a solid with one face exposed. Taking a small volume containing a single plane surface it can be shown (see appendix A) that the average number of particles crossing the interface during a diffusive time step is: where l 2 is the cross-sectional area of the face of the voxel. The diffusive flux, with units of number of particles per unit area per unit time is then: The mass balance at the interface requires that the diffusive flux must be balanced by the It follows that one voxel composed of N voxel mol mols of calcite is fully dissolved when the c 2016 American Geophysical Union. All Rights Reserved.
cumulative flux of particles (the total number of particles that penetrated that voxel) is equal to N hits : This establishes a direct relationship between the dissolution rate measured in the laboratory batch reactors and the flux of particles necessary to dissolve one voxel. Using the definition of λ, eq. 7: Thus, N hits is explicitly written in terms of physically meaningful variables: diffusion coefficient, time-step, particle concentration, moles of calcite per voxel and, crucially, the batch reaction rate. It represents the total number of particles necessary to dissolve a given amount of calcite. In realistic 3D geometries, solids with more than one face exposed dissolve faster than the ones with only one exposed face, since they receive particle hits through any of its free faces.
Micro and intragranular porosity may constitute a relevant portion of the pore-space in carbonate rocks, Lucia [2007]. Together with the surface rugosity, they are usually below the resolution of the micro-CT scanners. However, such sub-resolution rock attributes can be taken into account in our model by introducing the grain porosity (φ grain ), which corrects the number of calcite moles present in the solid phase. Grain porosity is calculated by comparing the helium porosity (φ HP ) with the micro-CT porosity (φ), the difference between both being attributed to sub-resolution porosity located in the grains: c 2016 American Geophysical Union. All Rights Reserved.
Using the data from Andrew et al. [2014] we have φ grain = 0.12 in the Ketton carbonate and the effective number of calcite moles in the solid voxels is: where n mol is readily calculated using the voxel volume (l 3 ) divided by the calcite molar volume 3.69 × 10 −5 m 3 .mol −1 , Wenk and Bulakh [2004].

Comparison with a dynamic imaging experiment
To validate and assess the predictive power of our method we test it against the results of a dynamic imaging experiment performed at reservoir conditions, see Menke et al. [2015] for a complete description. The experiment was performed on the Ketton carbonate, a relatively homogeneous, oolitic limestone. Figure  To establish the number of hits to dissolve (eq. 14), we take the independently measured intrinsic rate of calcite dissolution at the experimental pressure and temperature conditions using a rotating disc apparatus with a single calcite crystal -

Flow regime and dimensionless numbers
Transport and chemical phenomena in fluid flow are conveniently described by dimensionless numbers. The Péclet number (P e) compares the rate of advection to the rate of diffusion and is defined as: here q is the Darcy velocity and L = πV /S A is the characteristic length, where V is the volume of the sample and S A is the area of the pore-solid interface, Mostaghimi et al. [2012]. When P e > 1, transport is advection-dominated. The experiment was conducted at a fixed flow rate of 0.5 mL/min, corresponding to q = 660 µm/s; the same was also imposed during the simulations. At the beginning of the experiment P e = 2300, decreasing to 1500 after fifty minutes due to the increase in porosity.
c 2016 American Geophysical Union. All Rights Reserved.
The (advective) Damköhler number (Da) is defined as the ratio of the advective time and the reaction time, both calculated over the same characteristic length, Lasaga [1998]. Considering the advective time (t adv = Lφ/q), the dissolution time (t react = V n(1 − φ)/rS A ), and the definition of the characteristic length, we can write: where n = 2.39 × 10 4 mols.m −3 is the number of moles of calcite per cubic metre in the Ketton grains; n = n mol (1 − φ grain )/l 3 . At the experimental conditions Da increases from 3.2 × 10 −5 at the beginning to 5.0 × 10 −5 after fifty minutes, in both the experiments and simulation. However, the no-slip condition imposes V = 0 at the interfaces, thus the transport close to the solids is diffusion-dominated. Therefore the product P e.Da, which compares the reaction rate to the diffusion rate, provides a more useful characterization of the dissolution dynamics. P e.Da stays approximately constant during the entire experiment (and simulation) -P e.Da = 7 × 10 −2 . A summary of the image, flow and reaction parameters is presented in table 1.
Heterogeneous chemical reactions require both the transport of reactants to and from the solid surfaces as well as the reaction at the surface. The slowest of these is the limiting factor in the overall kinetics. There are two limiting cases to be considered: 1) the characteristic time of reaction is smaller than the diffusive characteristic time, meaning that the reaction is transport-limited; 2) The diffusive time is smaller than the reactive time, the reaction is limited by the reaction rate at the surface and dissolution is reaction limited. Under the experimental conditions, while advection dominates over diffusion, P e > 1, diffusion is faster than reaction in the vicinity of the grain surfaces, P e.Da < 1.
When P e is large and P e.Da small, as in our particular case, the solute renewal rate is c 2016 American Geophysical Union. All Rights Reserved.
high (conversely, the residence time is short) and the fluid is constantly out of equilibrium with respect to the solid. Hence, there is not enough time for the development of pH and concentration gradients within the pore-space. To illustrate this, we perform a simulation in a subsample of the Ketton carbonate, injecting 1.2 × 10 6 particles while keeping the average flow velocity the same as in the experiment. The results are shown in figure   3, which is consistent with the two-dimensional results of Li et al. [2008]. Furthermore, advection-dominated transport and reaction-limited reaction favour the development of a mostly uniform dissolution and the experimental data -section 4.3 -show that dissolution spreads across the sample. This is in agreement with previous results in the literature, e.g. Golfier et al. [2002] or Luquot and Gouze [2009].

Simulation procedure
To initialize the simulation, 2.0 × 10 4 particles are injected with a flow weighted rule at the inlet and the flow rate is kept constant at 0.5 mL/min. When a particle exits the domain either by diffusion or advection, or crosses the pore-solid interface, a new particle is injected at the inlet, keeping the number of particles constant and representing constant sample-averaged pH conditions. This is consistent with a high velocity flow where there is a constant replenishing of solute. It should be noted that, in our specific case, the rock dissolution is not limited by the amount of solute injected, but by the specifics of flow controlling the contact of the particles with the solid phase.
Due to the evolving pore space geometry caused by the reactive flow, the velocity field must be updated. Since dissolution is distributed throughout the sample, we assume a quasi-steady regime where small local increases in porosity do not lead to large variations in the global pressure field. The velocity field is updated only after the porosity increment c 2016 American Geophysical Union. All Rights Reserved.
(∆φ) is above a certain threshold, thus saving computational time. We established that ∆φ ≥ 0.001 is a good compromise between performance and accuracy; choosing a smaller ∆φ did not affect our results. We also studied the sensitivity of the results to the number of particles and the size of the time-step. We found that the time-step is sufficiently small for the results to be independent of its exact value. The number of hits to dissolve is proportional to the particle concentration, this imposes a lower bound on the number of particles, since we demand N hits > 1 throughout the sample. However, once this condition is met, there is no need for a large number of particles. The kinetics is independent of it: an increase in the number of particles leads to an increase in N hits , which is counterbalanced by an increase in J dif f . In fact, this is one advantage of the method, it does not require many particles to obtain meaningful results.
Given that the flow field is steady between updates and that reacting particles are killed and reinjected at the inlet, the profile of particle concentration along the direction of flow is also (quasi-) steady. The particle concentration is higher at the inlet than at the outlet. To make the reaction rate uniform throughout the sample, as suggested by the low value of Da and the experimental results, and to eliminate spurious boundary effects, we introduce a correction in the number of hits to dissolve. Equation 14 is modified replacing the volumetric particle concentration (C * ) with C * (s x ), the slice-averaged concentration of solute, where s x stands for the tomographic slice number in the flow direction. On average, N ef f hits ranges from 20 at the inlet to 11 at the outlet. Using eq. 19 the number of hits to dissolve becomes a function of the position along the sample axis, which is a necessary c 2016 American Geophysical Union. All Rights Reserved. assumption in steady state situations, when there is constant particle injection at inlet and the particle concentration is equilibrated. For all other scenarios (e.g. randomly injected particles) equation 14 should be used instead.
A summary of the algorithm is as follows: 1. Inject particles at the inlet using a flow weighted rule.
2. Advance particles along streamlines traced with the modified Pollock algorithm.
3. Perform diffusion by a random walk; if a particle hits a solid, reinject a new one at the inlet.
4. When the number of hits on a particular solid voxel is equal to N ef f hits (eq. 19) dissolve solid (change the solid to void) and update the porosity.
5. Repeat steps 2 -4 for every particle and time step.
6. Update the flow field when ∆φ ≥ 0.001 and update the number of hits to dissolve (eq. 19) to reflect the new particle concentration profile.
7. Repeat steps 2 -6 until the end of the simulation.

Results and discussion
We compare our predicted reaction-altered pore-space geometries to the ones obtained experimentally after fifteen, thirty three and fifty minutes of flow and reaction. Figure   4 shows the three-dimensional distribution of the dissolution front in the experimental and simulated images. A summary of image-calculated porosities and permeabilities is presented in table 2 and in figure 5; the agreement between simulation and experiment is very good. As we show later, we find preferential dissolution in some of the narrower pores and tighter portions of the rock that are initially constraining the flow -as a result the c 2016 American Geophysical Union. All Rights Reserved.
permeability increases sharply with a relatively small increase in overall porosity. After fifty minutes the experimental and simulated porosities increased approximately 50%, while the corresponding increase in the computed permeability was much higher, ≈ 500% in the experiment and ≈ 550% in the simulation. We attribute this overestimation in the simulated permeability to small, localized overpredictions in porosity that are amplified in the permeability calculation. Fitting a power law, we find that permeability scales as k ∝ and Vialle et al. [2014].
Not only the increase in average porosity is matched but also its variation across the sample, as seen in figure 6, which compares the experimental and predicted slice-averaged micro-CT porosities along the flow axis and their changes in relation to the porosity prior to flow. Figure 7 shows the relative changes in porosity -the change in porosity divided by the initial porosity -and we note a tendency for dissolution to be higher in the tighter regions of the rock, which correspond to the troughs in the slice-averaged porosity curve.
This could be an effect of the increase in solute flux in the pore-throats. Our simulation captures this aspect of the reactive flow and we point out that the positions of the peaks in the experimental and simulated curves coincide, despite some slight disagreement in the amplitude of the change. As discussed in the previous paragraph, the preferential dissolution of tight channels manifests itself as a steep increase in permeability.
c 2016 American Geophysical Union. All Rights Reserved.
In figures 8, 9 and 10 we compare slices -parallel and perpendicular to the core axis -taken from both the experimental and predicted 3D images after reactive flow. They depict the shrinking of the calcite grains and demonstrate an excellent agreement between model and experiment, as our model reproduces the measured dissolution patterns in both the radial and axial directions. It is also worth noticing that, as mentioned earlier, in high velocity flows such as this, the dissolution front is expected to reach the whole of the porespace, which is confirmed by the experimental results and reproduced in our simulations.
The increase in porosity is attributed to the calcite cement and grains being continuously dissolved and no evidence of wormhole formation is seen.
To illustrate the relationship between the flow field and dissolution, we plot the predicted change in pore structure after thirty three minutes and superimpose the local magnitude of the flow field, figure 11. The local dissolution rate is controlled by a subtle balance between advection, diffusion and reaction in the pore space. It appears that we see generally most dissolution in portions of the pore space with the greatest gradient in flow speed perpendicular to the solid surface. Diffusion allows solute to reach the solid, while a nearby fast, or even moderate, flow provides fresh solute for reaction. As discussed above, there tends to be more dissolution in regions of higher flux, where solute is forced through narrower restrictions, allowing greater contact with the solid surface. There is little dissolution in essentially stagnant regions of the pore space -those shown in white -since reactants can only access the solid through diffusion, a much more slow process. However, in some portions of the flow domain, there is no obvious relationship between pore velocity and degree of dissolution, with little change seen in some highspeed pathways (in red) and more dissolution in pores with slower-moving flow (blue).
c 2016 American Geophysical Union. All Rights Reserved.
This could be an effect of the coupled nature of dissolution and flow, as a narrow throat is dissolved and widened, its flow velocity must drop, since the flow is mass conservative. The cumulative dissolution is the result of the availability of particles in a particular portion of the sample during a time interval. However, since the geometry, and consequently, the flow field are constantly changing, the instantaneous velocity may not be a good approximation for the velocity at previous times. A more quantitative analysis of this phenomenon is a possible subject for future work.
The increase in the average pore size is also reflected in the velocity probability distribution, figure 12. With increasing connectivity, the velocity distribution after dissolution is narrower: velocities are more uniform and closer to the average with fewer stagnant regions. The velocities are distributed across four orders of magnitude, indicating a wide range of local ratios of diffusive over advective time. This suggests that part of the porespace will not be accessible to the reactants at the time scale relevant for the experiment.
Overall, the agreement between the simulations and experiment indicates that, at the millimetre scale, the reaction is restricted by the transport of solute to and from the solids.
Our simulations have no adjustable parameters: our inputs are the initial pore geometry, as imaged using micro-CT scanning, the experimentally-imposed flow rate, the diffusion coefficient, the calcite volume in the carbonate grains and the independently-measured batch reaction rate. Also note that we have assigned an average porosity to the grains -in reality the degree of micro-porosity may vary from grain to grain -and we do not resolve small-scale roughness in the pore space. As a consequence, the ability to predict local changes in porosity -on a pore-by-pore basis -also indicates that the reaction is controlled c 2016 American Geophysical Union. All Rights Reserved.
by the availability of reactants, rather than by the local sub-resolution properties of the solid interface.
In addition to porosity and permeability, another property influenced by changes in porous structure is the tortuosity, a measure of the complexity of the pore space geometry, with higher values indicating more complex geometries. Making use of the semi-analytical method described in Pereira Nunes et al. [2015], and momentarily ignoring diffusion, we trace streamlines from the inlet to the outlet of the sample and calculate the hydraulic tortuosity (T ), defined as [Duda et al., 2011]: where L sl is the average length of the streamlines and L x is the length of the sample (3.5mm). Tortuosity was found to decrease from 1.54 prior to injection to 1.25 in the experimental image and to 1.22 in the simulated image after fifty minutes of flow and reaction. Such decrease reinforces the scenario of a permeability increase caused by the expansion of pore-throats -figures 8, 9 and 10 -and is consistent with the expected decrease in complexity of a rock experiencing an almost uniform dissolution. As the porosity increases, more pores become connected and the length of the fluid paths from inlet to outlet decrease. Of course, in the extreme case of complete dissolution, the tortuosity would decrease to unity.

Average dissolution rate vs batch rate
The average effective reaction rate in the sample can be computed from the measured surface area and the change in porosity, which is directly proportional to the number of moles of calcite that have dissolved. Between times t and t , we write the effective c 2016 American Geophysical Union. All Rights Reserved. dissolution rate as: where ∆φ = φ(t ) − φ(t) and S A (t) is the surface area at time t, which is assumed to be slowly changing, so that S A (t) ≈ S A (t ). V and n have been defined in section 4 and are the volume of the sample and the number of moles of calcite per cubic metre in the Ketton carbonate grains, respectively. Table 3 shows the measured and predicted increments in porosity and the corresponding reaction rates. We see that the effective reaction rate, calculated from the experimental images, is approximately 10 times slower than the batch rate. Analysing the predicted rates, calculated from the simulated images, we would be drawn to the same conclusion. This is in a system with negligible chemical heterogeneity and where the agreement between experiment and prediction implies that all the solid surface is available for reaction. This demonstrates that the heterogeneity of the flow field itself acts as a significant limitation on the transport of solute to the surface, even in a regime where diffusion dominates over reaction at the pore scale (P e.Da 1).
The Ketton limestone studied is, by comparison with most reservoir rocks, or indeed other quarry carbonates, relatively homogeneous [Bijeljic et al., 2013b], yet we see transport-limited reaction which is one order of magnitude slower than its batch value when averaged over a sample only a few millimetres across, even when the effects of chemical heterogeneity and effective reactive surface area are ignored. We conjecture that the origin of this effect is the wide variation in local flow speeds, as illustrated in figures 11 and 12. This variation in the velocity implies that the residence times and the access of solute to the solid surface due to diffusion is not evenly distributed throughout the pore c 2016 American Geophysical Union. All Rights Reserved.
space. In an advection-dominated flow parts of the sample are not be available to the reactants in the time frame of the experiment and simulation.
Diffusion and reaction locally deplete the particle concentration near the solid; replenishment of particles is controlled by both advection and diffusion occurring over a large range of time scales through a tortuous pore space geometry. The overall effect is a substantial decrease in mixing (here understood as the amount of solute in contact with the solids) and reaction compared to a well-mixed system. At larger scales the effective reaction rate may be several orders of magnitude lower still. This is the combined effect of several kinds of heterogeneity, both physical and chemical, at different spatial and time scales. Our results indicate that, when transitioning from the micron to the millimetre scale, flow field heterogeneity plays a significant role in decreasing the average rate. However, to upscale this effect and help explain the observed discrepancy between batch and field-averaged rates referred to in the Introduction, further work in larger systems needs to be done.

Conclusions
We have presented a pore-scale calcite dissolution method capable of predicting the evolution of petrophysical properties in carbonate rocks subjected to CO 2 -saturated brine injection at reservoir conditions.The method is applicable to advection-dominated flows and reaction-controlled kinetics, when the fluid in constantly out of equilibrium with the rock. Under these conditions we simplified the multispecies problem of calcite dissolution to a single species reaction. The method is based on a Lagrangian approach where particle advection is performed with a semi-analytical streamline tracing algorithm that obeys the no-flux condition at pore-solid walls with diffusion incorporated using a random walk c 2016 American Geophysical Union. All Rights Reserved.
method. The fact that streamlines are traced semi-analytically greatly improves the computational efficiency of the code and eliminates numerical stability issues. In contrast with previous results in the literature, we have approached the problem of modelling dissolution in porous media using an image of a real, three-dimensional carbonate rock and have validated our results against an experiment performed at in situ conditions.
The deviation from equilibrium saturation of calcium is represented by the particle concentration and the local flux of particles is the mechanism controlling calcite dissolution.
We have derived a relationship between the cumulative flux of particles required to dissolve one solid voxel and the batch dissolution rate of calcite in acidic brine at reservoir conditions. We have validated our method against experimental data from a dynamic imaging dissolution experiment performed at high temperature and pressure without any adjustable parameters. The predicted increases in average porosity were matched with a high degree of accuracy and the simulated porosity-permeability trend closely follows the experimental one. The dissolution front is reproduced throughout the sample in directions parallel and perpendicular to the flow axis. The rapid increase in permeability is explained by the preferential dissolution of tighter portions of the rock. Such alterations in pore space geometry are also reflected in the decrease of the hydraulic tortuosity. We argue that the overall dissolution effect not only increases porosity and permeability but also reduces the pore-space complexity.
We have demonstrated a decrease of one order of magnitude in the effective dissolution rate of a porous carbonate sample compared to the batch dissolution rate. There is no diffusive limitation on the reaction at the grain surfaces, suggesting that any eventual transport limitation manifests itself not at the pore scale, but rather at the sample scale.
c 2016 American Geophysical Union. All Rights Reserved.
This millimetre-scale transport limitation controls the local availability and residence time of the solute particles, and its overall result is a decrease in the sample-averaged reaction rate. This occurs even with all the surface available to reaction and in a relatively homogeneous and chemically uniform sample. Our results indicate that the very large differences observed between batch and core-scale reaction rates could, in large part, be explained by the inhomogeneity in the flow field and the consequent transport-limited flux of reactants at the solid surface. The local variations in dissolution are the result of a complex interplay between the velocity field, that carries the solute from the inlet, diffusion, responsible for the flux across the pore-solid interface, and the pore-space geometry, that may influence the local concentration of solute and controls the effective reactive surface area.
We suggest that our pore-scale results give valuable insight into the processes governing carbonate dissolution and may provide a starting point to the refinement of upscaling techniques for reactive flows. As our results clearly show, flow field heterogeneities are a potential source of the discrepancy between millimetre-and micron-scale dissolution rates.
Moreover, our methodology can be readily employed to simulate a range of dissolution patterns in different carbonates caused by different flow regimes and investigate further the mechanism by which the complexity of the flow field affects the average reaction rate in heterogeneous porous media.
c 2016 American Geophysical Union. All Rights Reserved.

Appendix A: Average Flux of Particles in a Diffusion Time-Step
Here we calculate the average number of particles crossing the solid-fluid interface during diffusion, eq. 6. A particle located at a distance s from a solid voxel face crosses it when: which implies that s ≤ ∆X dif f i ≤ λ/2 and that s ≤ λ/2. Since ξ is uniformly distributed over 0 and 1, the probability of crossing is: If the local concentration of particles is C and the cross-sectional area of the voxels is l 2 , the number of particles that cross this particular interface during the diffusive step is: Acknowledgments.  Table 3.
Image-based experimental and simulated porosity increases and the inferred av-  Pictorial representation of the advection-diffusion-reaction process. A particle moves along the analytically-traced streamlines during advection (black arrows), jumps from one to another during diffusion (red arrows) and crosses the pore-solid interface during reaction. The total displacement during a time-step is the sum of the advective and diffusive displacements.
White cells are pores and black are solid.
c 2016 American Geophysical Union. All Rights Reserved.

Figure 2.
Top, three-dimensional segmented image of a Ketton carbonate superimposed with the velocity (logarithmic colorbar) along the flow direction. The pore-space is shown in grey. Bottom, streamlines traced using a pore-scale semi-analytical method in the pore space reconstructed from the micro-CT and coloured by vertical coordinate. Solute particles move from