Reservoir-condition pore-scale imaging of dolomite reaction with supercritical CO2 acidified brine_ Effect of pore-structure on reaction rate using velocity distribution analysis

To investigate the impact of rock heterogeneity and flowrate on reaction rates and dissolution dynamics, four millimetre-scale Silurian dolomite samples were pre-screened based on their physical heterogeneity, defined by the simulated velocity distributions characterising each flow field. Two pairs of cores with similar heterogeneity were flooded with supercritical carbon-dioxide (scCO2) saturated brine under reservoir conditions, 50 °C and 10 MPa, at a high (0.5 ml/min) and low (0.1 ml/min) flowrate. Changes to the pore structure brought about by dissolution were captured in situ using X-ray microtomography (micro-CT) imaging. Mass balance from effluent analysis showed a good agreement with calculations from imaging. Image calculated reaction rates (reff) were 538 times lower than the corresponding batch reaction rate under the same conditions of temperature and pressure but without mass transfer limitations. For both high (Péclet number = 2600-1200) and low (Péclet number = 420-300) flow rates, an impact of the initial rock heterogeneity was observed on both reaction rates and permeability-porosity relationships.


Introduction
Understanding reactive transport in subsurface formations is important in many applications including carbon capture and storage (CCS) (Bachu, 2000), well acidization (Wang et al., 1993), contaminant transport (Zheng and Bennett, 1995), and leaching (Lin et al., 2016b). In the context of CCS in carbonate rich aquifers, dissolution of minerals raises a concern regarding CO 2 storage security (Metz et al., 2005). Field-scale flow modelling that incorporates reaction processes between the CO 2 , formation brine, and rock is a necessary precaution and may help control uncertainties which affect storage security (Sifuentes et al., 2009). However this task proves challenging as there are large variations in reported values for parameters, resulting in predicted reaction rates that span orders of magnitude (Black et al., 2015). Therefore, in such models, a fundamental pore-scale understanding is key in accurately representing the behaviour of different length-scale heterogeneities and forms the basis upon which field-scale models are built (Li et al., 2008;Steefel et al., 2015;Steefel and Lasaga, 1994;Steefel and Maher, 2009). It is common to use a power law to predict permeabilityporosity relationships, calibrated by its exponent − denoted in this work by m (Carman, 1937;Kozeny, 1927). This relationship is needed to characterize the effect of dissolution on averaged flow properties as input into field-scale models. Menke et al., 2016Menke et al., , 2015Noiriel et al., 2004;Smith et al., 2017Smith et al., , 2013. Pore-scale observations can help explain the mismatch between laboratory and field scale reaction rates (Li et al., 2006;Moore et al., 2012;Noiriel et al., 2004;Swoboda-Colberg and Drever, 1993;White and Brantley, 2003) and help elucidate the relationship between permeability and porosity for different initial pore structures.
Carbonate reservoirs, most susceptible to dissolution, are complex heterogeneous systems with mixed mineralogies, each having different reaction rates and characteristics (Black et al., 2015). Rock heterogeneities, fluid transport and dissolution kinetics are all factors which influence macroscopic dissolution patterns (Noiriel et al., 2009;Steefel and Lasaga, 1990). In recent years, direct pore-scale visual observations of geologic rocks, including dissolution, have been made possible by the use of X-ray imaging (Al-Khulaifi et al., 2017;Gharbi et al., 2013;Gouze and Luquot, 2011;Lin et al., 2016a;Lin et al., 2017;Luquot and Gouze, 2009;Menke et al., 2016Menke et al., , 2015Noiriel et al., 2009Noiriel et al., , 2005Noiriel et al., , 2004. The way in which rock property alterations develop from dissolution is dependent on mineralogy, reservoir conditions, and pore structure (Al-Khulaifi et al., 2017;Noiriel et al., 2005;Wang et al., 1993). A systematic way to study reactive transport in carbonates is to start with reaction in pure minerals and then characterise dissolution dynamics in more complex mixtures. In single mineral systems, dissolution-induced changes in dynamic properties have been shown to be dependent on initial pore structure (Maheshwari et al., 2013;Menke et al., 2016;Smith et al., 2017). Thus far, most single-mineral imaging of dissolution by CO 2 acidic brine have featured calcite bearing limestones (Gharbi et al., 2013;Gouze and Luquot, 2011;Luquot and Gouze, 2009;Menke et al., 2016Menke et al., , 2015Noiriel et al., 2005Noiriel et al., , 2004. They observe dissolution regimes that vary depending on the influences of transport limitation, reaction kinetics, pore and surface heterogeneities. In the works of Noiriel et al. (2004Noiriel et al. ( ), (2005, the dissolution of a mmscale limestone core by flowing ambient condition CO 2 -enriched brine was captured with micro-CT imaging and chemical effluent monitoring. The intent was to identify the inherent relationship between fluid chemistry, hydrodynamic and structural properties by assessing porosity and reactive surface changes in parallel to permeability changes. This new type of experimentation enables the documentation of permeability-porosity evolution and its relation to structural modification. The study's main conclusions were that the initially large increase in permeability was observed by imaging and was attributed to the removal of fine grain calcite particles. Following this, the permeability increased with porosity at a lower rate because of an observable decrease in wall roughness and increase in pore connectivity. Finally, a higher flow rate gave a higher dissolution rate due to transport-limited conditions that accelerate dissolution in narrow pore pathways. In another work, Luquot and Gouze (2009) conducted a series of reservoir condition core flooding experiments on a limestone reservoir rock using reactive brine with different partial pressures of CO 2 to observe the effect of CO 2 injection away from an injection well. Non-uniform dissolution was reported for transport-controlled reaction while lower saturations of CO 2 , seen farther away from the injection well, exhibited reaction-controlled uniform dissolution. Luquot and Gouze (2009) also reported an initial pore structure memory effect the loss of which offsets the power law relationship between permeability and porosity.
Fewer studies have looked at a slower reacting dolomite (Smith et al., 2017) despite the fact that dolomite is the second most abundant carbonate mineral (Ahr, 2008), and features as the main mineral in some prominent hydrocarbon reservoirs (Al-Siddiqi and Dawe, 1999). Smith et al. (2017) conducted reactive transport type experiments on dolomite samples to constrain 3D models that predict the evolution of porosity and permeability for CO 2 storage. The work compared the experimentally determined power law exponent of the permeabilityporosity relationship and that from a model.
The approach in this work allows for the quantification of heterogeneity enabling direct comparison of its effects on dissolution patterns and perceived reaction rates. Therefore, the impact of pore structure on reaction rates and dissolution patterns in Silurian dolomite reacting with CO 2 acidified brine at reservoir conditions will be examined.
The general kinetics of any acidic brine-rock interactions at the pore scale are dictated by the interplay between hydrodynamic transport of reactants to fluid-rock the interface, reaction at the interface, and the transport of reactants away from the interface (Golfier et al., 2002;Lasaga, 1984). Studying the impact of pore structure heterogeneity for both single and mixed mineralogy systems is also possible using micro-CT imaging (Al-Khulaifi et al., 2017;Smith et al., 2013). Al-Khulaifi et al. (2017) imaged the dissolution of a millimetre-scale composite core made of Ketton calcite and Silurian dolomite arranged in series. The aim was to examine the effect of physical and chemical heterogeneity on dissolution regime and effective reaction rate; the average reaction rate of dolomite under similar conditions to those presented later in this work and a Péclet range of 3000-1000 was 6.46 × 10 −6 mol/m 2 s. Menke et al. (2016) injected super-critical CO 2 saturated brine at 50°C and 10 MPa into millimetre-scale calcite cores (Portland and Estaillades, > 97.2% calcite), with variations of physical heterogeneity between the samples and imaged the pore space changes resulting from dissolution. In more heterogeneous rock it was observed that dissolution occurs along preferential flow pathways leading to channel growth. On the other hand, homogeneous pore structures gave rise to a uniform dissolution i.e. an even growth of all flow pathways.
Dolomite exhibits a reaction rate one order of magnitude lower than that of calcite (Peng et al., 2016). In the context of CO 2 storage, Luhmann et al. (2014) performed core-scale (centimetre-scale) dissolution experiments in dolomite core plugs to investigate the effect of varying flowrate; a faster increase of permeability was observed at high rates. In a more recent study, Smith et al. (2017) recognized the importance of pore-scale observations and imaged the dissolution of dolomite cores with different permeability and physical heterogeneity at 43 μm resolution to calibrate a continuum-scale reactive transport model. It has been reported that physical heterogeneity drives the dissolution front Hao et al., 2013;Kalia and Balakotaiah, 2009;Maheshwari et al., 2013;Smith et al., 2013) that determines the porosity-permeability relationship and effective reaction rate.
The approach in this paper is to characterize reactive transport in dolomite for CO 2 storage by combining pore-scale experimental observation and modelling. We provide micro-scale observations at an image resolution of 5.2 μm to provide insights into coupled reactive phenomena unresolvable at larger scales. To make direct observations of the effect of physical heterogeneity and flow conditions on reaction rates and patterns, a pre-screening method for characterizing the rock heterogeneity prior to conducting the reaction experiments is introduced. First, dry scans of the core samples are performed. Then, on these images, direct numerical simulations using a Navier-Stokes flow solver are run Raeini et al., 2012) to obtain velocity fields as their intrinsic heterogeneity characteristics. Using this analysis, four dolomite samples are classified into two heterogeneity groups A and B, based on their pore-scale velocity distributions. Next, dynamic reactive transport experiments are performed to obtain effective reaction rates and examine the impact of pore structure and initial flow field heterogeneity at two flowrates. Reactive transport experiments are then analysed by again employing direct flow simulation on time-lapsed pore-space dissolution images acquired during experimentation to provide (a) the porosity-permeability relationship, and (b) dynamic changes in flow heterogeneity characteristics from reaction-induced alterations in pore velocity distributions. Finally, there will be a discussion the origin of the effects of flowrate and initial rock heterogeneity on dissolution rates and patterns.

Materials and methods
The Silurian dolomite ( > 99% dolomite) used in the reaction experiments was from the Thornton formation deposited between 443-416 million years ago (Al-Khulaifi et al., 2017). Silurian is homogeneous in chemical composition but spatially has a very heterogeneous pore structure even at the scale of sample size chosen for this work, 4.8 mm long by 4.8 mm diameter. As a pre-screening step, 8 samples were cored out of a single block of Silurian dolomite with the length of the cores being parallel to the bedding of the rock. After imaging the 8 samples at a resolution of 5.2 μm, flow simulation was performed directly on their segmented pore spaces (detailed description of flow modelling in Section 2.5). Characterizing the intrinsic flow heterogeneity of each sample was achieved by observing the probability density functions (PDFs) of velocities calculated in each pore voxel from the simulation. Two pairs were chosen for experimentation based on the shape of their velocity distributions reflecting their structural heterogeneity. Each individual pair has samples of matching heterogeneity (heterogeneity A and heterogeneity B) and within each pair, a low (0.1 ml/min) and high (0.5 ml/min) rate dissolution experiment was performed, giving a total of four experiments.

Experimental methodology
Loading a core into the coreholder begins by wrapping a single layer of aluminium foil tape around the core (leaving the ends open) to prevent CO 2 gas from leaching laterally out of the core when reservoir conditions are prevalent. The core is then inserted snugly into a Viton rubber sleeve of which both ends are capped by flow end pieces. Next, two more layers of aluminium tape are tightly wrapped around the Viton-end piece assembly, to further prevent CO 2 gas from leaching out. During the experiment, a 2 MPa pressure differential is applied to confine the samples. Finally, the Viton assembly is put into an X-ray transparent carbon-fibre Hassler type core holder. The core holder is wrapped in heating tape to maintain the core at a temperature of 50°C; this is controlled by a thermocouple threaded along the core. Fig. 1 shows a schematic of the core holder assembly.
To achieve reservoir conditions, a Hastelloy reactor (Parr instruments Co.) is filled with water containing 5 wt% NaCl and 1 wt% KCl. It is pressurized to 10 MPa by injecting scCO 2 via the injection pump, heated to 50°C, and left to equilibrate overnight. A mechanical stirrer ensures that the brine is fully saturated resulting in a final pH of 3.1 (Peng et al., 2013). Finally, a two position actuator with sampling loop (VICI Valco instruments Co. Inc.) is connected to the flow line downstream the core holder, just outside the confines of the micro-CT enclosure, for periodic effluent collection. Fig. 2 displays a schematic of the experimental apparatus.
The experiment starts by flooding the core with a 6 wt% potassium iodide (KI) brine by back-flowing it from the receiving pump side until it drips out of the 3-way valve at the bottom of the reactor. This brine attenuates more X-rays causing the pores to show brighter in the images. Next, the core sample and flow lines are pressurized to 10 MPa by the receiving pump. The reactor, which has been maintained at 10 MPa pressure from earlier, is opened to the flow line at the 3-way valve. Finally, the receiving pump is set to refill at the desired flowrate of the experiment to begin reactive brine flow. Experimental time began when projections showed the bright KI leaving the core, signalling the arrival of reactive brine. At that precise moment, the sequence of scans is started using the parameters shown in Table 1.
To complete each experiment, nine images are taken (one dry scan to capture the initial pore structure and eight subsequent scans to capture changes caused by dissolution) (Fig. 3). The intrinsic reaction rate of dolomite, that is the batch reaction rate of 5.1 × 10 −5 mol/m 2 s measured with scCO 2 saturated brine at the same conditions as in this Fig. 1. Carbon-fibre core holder assembly. The carbon-fibre sleeve allows for the experiments to be performed at reservoir conditions while being sufficiently x-ray transparent to enable the imaging of the core inside. Fig. 2. Experimental apparatus design of which the main parts are a micro-CT scanner, 3 pumps, a reactor vessel, an effluent sampler and carbon-fibre core holder. The only part inside the scanner's lead enclosure is the core holder; all lines to and from the core holder go through an x-ray baffle located on the door of the scanner's hutch. An effluent sampler is connected to the line downstream the core holder near the baffle on the exterior side.
work (Peng et al., 2016), is sufficiently slow for a laboratory based micro-CT (here an Xradia VersaXRM-500) to capture the dissolution dynamics. Each scan takes 50 min and an effluent sample is collected at the end of each scan for the subsequent analysis by means of inductively coupled plasma mass spectrometry (ICP-MS).

Image processing
Each raw 3D image captured during experimentation was filtered to reduce image noise using an edge preserving nonlocal-means algorithm (Buades et al., 2005). Next, the images were segmented into pore and rock phases with a seeded watershed algorithm in Avizo 9 (Schlüter et al., 2014).

Effective reaction rate and Péclet-Damköhler analysis on images
The dimensionless numbers, Damköhler (Da) and Péclet (Pe), are used to characterize and quantify reactive transport. The Damköhler number characterizes the relationship between the chemical reaction timescale to the mass transport timescale of a particular system: = reaction rate mass transport rate Da (1) The most common form of the Damköhler number relates the reaction rate to the mass transport rate through advection or mass flow; this is called the advective Damköhler number: where k is the chemical reaction rate constant (1/s), L c is the characteristic length (m), and u avg is the interstitial velocity (m/s) defined as the average velocity of the fluid flow within the pores.
u D is the Darcy velocity (m/s) and ϕ CT is the porosity calculated from the segmented micro-CT images. Da is also defined as the ratio of the time for advection to time for reaction over a characteristic length and reaction time can be interpreted as the time taken to dissolve a portion of the rock of length L. As stated earlier, the non-transport limited reaction rate of dolomite under condition of this work was measured to be r = 5.1 × 10 −5 mol/m 2 s. If Silurian dolomite was assumed to be pure dolomite of spherical grains, then the reaction rate constant can be estimated by k = πr/nL (Menke et al., 2015), where n represents moles of dolomite per unit rock volume: ρ dolomite is the grain density of dolomite measured from helium porosimetry (2820 kg/m 3 ) and M dolomite is its molecular mass (0.1844 kg/ mol). ϕ total is ϕ CT plus the unresolvable microporosity (ϕ micro ): where (ϕ micro ) is estimated from helium porosimetry using the equation: For the particular Silurian samples used in the experiments, ϕ micro was calculated to be 0.05, assuming a homogeneous grain porosity throughout the rock. With the above estimation of reaction rate constant, Eq. (2) can be rewritten as follows (Menke et al., 2015): The Péclet number compares the effectiveness of mass transport of the solute by advection to that by diffusion: This definition leads to Péclet number as an estimation of fluid movement over a characteristic length (Péclet 1827): D m is the temperature dependent molecular diffusion coefficient (7.5 × 10 −10 m 2 /s for Ca 2+ at 25°C) ). The characteristic length (L c ) of consolidated rock samples is found from the ratio of sample bulk volume V B to micro-CT image calculated surface area A s i.e. specific surface area S (m −1 ) (Mostaghimi et al., 2012): The effective reaction rate (r eff − mol/m 2 s) of dolomite dissolution, including mass transfer limitations, was calculated from porosity changes seen between each consecutive segmented micro-CT image (Al-Khulaifi et al., 2017): Although the reactive brine can invade microporosity by diffusion, it is assumed that this process is too slow and does not cause dissolution to take place within the grainsthis is confirmed by the images as no observable change in the grey scale of the rock grains. The grain porosity term ϕ micro in the equation for r eff accounts for the fact that the solid contains micro-porosity. Peng et al. (2016) reported the batch non-transport limited reaction rate of dolomite mineral to be 5.1 × 10 −5 mol/m 2 s at similar solution composition and temperature and pressure conditions as in this study.
Eq. (11) captures an average rate of mass transfer in the experiment: we divide the mass transfer by the specific surface area from image analysis to find an effective reaction rate (Li et al., 2008;Luquot and Gouze, 2009;Menke et al., 2016Menke et al., , 2015Noiriel et al., 2004;Steefel and Maher, 2009). This reaction rate is an effective value and accounts for mass transfer limitations in the dissolution process at the specific conditions of the experiment. We follow this approach to allow direct comparison with experimental measurements from batch reactors when mass transfer limitations are avoided through vigorous stirring.

Effluent analysis
To complement the assessment of dissolution from imaging, the effluent was sampled periodically over the course of the experiments and analysed with ICP-MS. A limiting factor to the sampling frequency is the scan duration. Each time a sample is taken, the flow line sees a pressure blip which may cause the core to move in the coreholder, disrupting image acquisition. Therefore, effluent samples were taken between scans to avoid disturbing the image acquisition process. Eight 2 ml samples were collected for each experiment, one at the end of each reactive brine flow scan. The samples were measured for Mg 2+ ion concentration and the number of moles of Mg 2+ in a sample is taken to be equal to the number of moles of dolomite dissolved in the sample based on dolomite's molecular formula, MgCa(CO 3 ) 2 . The dissolved dolomite volume in the period between two sampling points (50 min) was calculated assuming a linear relationship between each measured concentration: where n labels the sample number, ΔV is the dissolved volume of Silurian dolomite (m 3 ), c is the sample's dolomite ion concentration in (mol/m 3 ), t is the sampling time (s), and Q is the flowrate (m 3 /s). Using this method, the amount of dolomite dissolved during the first 50 min of dissolution cannot be accurately measured because c 0 = 0, and is therefore omitted.

Flow modelling and rock heterogeneity characterization
Simulations on the segmented micro-CT images were performed using a finite volume solver implemented in OpenFOAM Raeini et al., 2012). The simulation simultaneously solves the Navier-Stokes and volume conservation equations: u is the velocity and p the pressure calculated in every pore voxel. The solver uses the viscosity and density of fresh water for μ and ρ at standard ambient conditions. Flowrate (m 3 /s) is calculated as Q = ∫ u x dA x where A x is the cross-sectional area of the core perpendicular to the direction of flow in (m 2 ) and u x is the face velocity parallel to flow in (m/s). The Darcy velocity is calculated as u D = Q/(L y L z ) where L y and L z are the image lengths in the y and z direction. From the simulation, probability density functions (PDFs) of velocities were obtained for the pore voxels in each segmented image and these distributions are used to characterize the intrinsic flow heterogeneity of the four samples. Fig. 4 shows the velocity distributions, where velocity is sampled logarithmically, for the initial state of the four samples. Samples A are observed to have lower pore structure heterogeneity (single peak distribution) and a smaller number of low velocities, whereas samples B are seen as having a higher heterogeneity due to the double peak of their velocity distributions and a larger number of stagnant velocities. The range of flow speeds is similar, however, for all the samples and consistent with all 8 samples considered in the screening study mentioned in the start of Section 2.
For each pair of matching heterogeneities, one sample was flooded with a low rate (0.1 ml/min) and the other with a high rate (0.5 ml/ min) during their respective experiments.
The flow solver was also used to calculate the permeabilities on the segmented micro-CT images using the Darcy equation: where ∇p is the average pressure gradient across the core (Pa).

Results
In Section 3.1, presented first are the porosity changes during dissolution and then the results of the analysis of dissolved rock volume from micro-CT images are compared with the measurements obtained from ICP-MS effluent analysis. Next, in Section 3.2 effective reaction rates from measurements are presented, based on the rate of mass transfer under these experimental conditions, followed by a pore-scale simulation analysis of the permeability-porosity relationship and dynamic changes in Pe and Da in Section 3.3. Finally, in Section 3.4, the flow fields of the four dolomite samples are visualized and will serve as a basis for interpretation of the effective reaction rates based on timedependent velocity distributions presented in the Discussion.

Evolution of porosity and material balance analysis
The increase in porosities, due to dissolution, calculated on images For samples A and B, faster increases in porosity are measured for the higher rate experiments. This implies that flowrate has a profound effect on dissolution in dolomitein both Section 3.4 and the Discussion, we will use pore-scale simulation to explain the origin of this behaviour for different rock heterogeneities. The evolution of porosity profiles along the length of the cores, shown in Fig. 6, provide more detail on how porosity increases; both porosity and change in porosity are represented.
The initial distribution of porosities along the length of the samples, shown by the blue lines in Fig. 6, is a reflection of Silurian dolomite's heterogeneous pore space. The plots of change in porosity show more dissolution at the inlets of the cores for the high rate experiments for samples A and B.
Magnesium ion (Mg 2+ ) concentrations measured in the effluent are used to estimate the rock volume dissolved and this is compared to the rock volume dissolved from the change in porosity on the micro-CT images. Fig. 7 compares the volumes calculated from both methods for all four samples.
A good agreement is seen between the results obtained by the two methods. Oscillations in the plots of volume dissolved measured from images could very well be real or caused by errors in phase segmentation. Either way, a one-to-one comparison with the effluent method is not entirely correct because it is measured by averaging dissolution over 50 min periods from sampling points of 2 ml volumes. Rather, the intent is to see that the overall trends between the two methods are the same. For the high rate conditions, approximately two times as much of the rock has dissolved than at the lower rate. For the low rate conditions, there is a very slow decrease in the amount of dolomite disolving throughout the experimental time whereas for the higher rate, the amount of dolomite dissolution decreases more markedly over time. increase in porosity from dissolution than the low rate experiments. For the low rate experiments, the porosity increases at similar rates whereas at a higher rate, the more heterogeneous sample (B) shows a more rapid increase in porosity than sample A but arrive at the same value at the end.  3.2. Effective reaction rates Fig. 8 shows effective reaction rates versus time for the four dolomite core experiments calculated from the images using Eq. (11): there is relatively little change in rate with time. The effective reaction rates for samples A and B at the higher rate are 5.47 × 10 −6 mol/m 2 s and 10.60 × 10 −6 mol/m 2 s respectively, while for the lower rate they are 3.15 × 10 −6 mol/m 2 s and 1.35 × 10 −6 mol/m 2 s respectively.
From Fig. 8 it is observed that higher flowrates result in higher effective reaction rates due to more efficient transport. Reaction is faster in the more heterogeneous sample B which is attributed to the formation of several preferential dissolution channels that were quantified as the two peaks in the high velocity regions of velocity distributions in Fig. 4 in Section 2.5. This will be explained in more detail by visualising time-dependent velocity fields in Section 3.4 and velocity distributions in the Discussion. Fig. 9 shows the permeabilities computed through direct simulation of flow by the method described in Section 2.5 against the connected porosities calculated from the micro-CT images.

Permeability simulation and Péclet-Damköhler analysis
For both samples A and B a much a greater overall increase in permeability is observed for the higher rate experiments. Irrespective of initial rock heterogeneity and initial porosity/permeability, the porosities and permeabilities at the end of core flooding in the high flowrate experiments have similar values. In contrast to this, Fig. 9 shows that at the end of the low rate experiments the permeability for the initially more heterogeneous rock B is lower than that of the sample A by a factor of almost 5. This indicates that initial rock heterogeneity and initial porosity/permeability have an appreciable impact on the dissolution rates under less advective conditions. Although smaller changes in porosity from the low flowrate experiments are induce larger jumps in permeability, the final permeabilities of the high rate experiments are greater due to the greater overall increase in porosity. Part of this phenomenon can be attributed to the effect of transport and part to the difference in pore structure between the two samples.
A power law with the equation = K cϕ m connected is used to fit the porosity-permeability relationships where K is the permeability (m 2 ), c is a constant and m is the power-law exponent. The use of this relationship has been established in earlier works to help parameterize permeability-porosity functions for modelling (Noiriel et al., 2004) and observed experimentally in the context of mineral transport limited dissolution (Al-Khulaifi et al., 2017;Luquot and Gouze, 2009;Menke et al., 2016Menke et al., , 2015Noiriel et al., 2004;Smith et al., 2017Smith et al., , 2013 with values of m in the range 0.29-21. The exponents obtained in the work presented here are in the upper range and higher. It has been observed that lower rates yield higher m values; a smaller increase in porosity yields a greater jump in permeability. As shown later, unstable dissolution fronts are seen and they are attributed to pore space heterogeneity which is the origin of this sensitive dependence of permeability on changes in porosity (Smith et al., 2017). The next section will use velocity field simulation to further examine the effects of structural heterogeneity. Fig. 10 shows the results of the Péclet-Damköhler analysis, as described in section 2.3. All experiments are considered to have relatively high Péclet number meaning that advection dominates over diffusion. However, the degree to which the dissolution is transport controlled varies between the two rates as evident by the difference in effective Fig. 7. The volume of dolomite dissolved in 50 min increments for the samples calculated from ICP-MS effluent analysis and from micro-CT images; a good match is observed between measurements of both methods. The first points from effluent analysis could not be quantified using this method because the initial effluent samples do not have representative concentrations due to the delay in the arrival of the reactive brine to the sampling point after t 0 .  All the experiments show a decrease of Péclet number through time; a faster decrease is observed in the high flowrate cases, which is due to the faster dissolution, resulting in a larger increase in porosity and hence a greater lowering of u avg for a fixed Darcy velocity. The Damköhler number is much less than 1 for all flowrates implying that the reaction rates are slow relative to advection. All experiments indicate that Da increases over time in a steady fashion − this is a consequence of the decrease in average flow speed over the experiments; the reaction rates do not change so markedly, see Fig. 8. Tables 2-5 show the flow, transport and reaction parameters (Pe, Da, ϕ CT , r eff , S, and permeability) calculated from the segmented micro-CT images for all four experiments at t = 0 min, t = 200 min, and t = 400 min.

Pore scale velocity fields
Simulation on the segmented images' connected pore voxels provides the velocity fields and the distribution of velocities. A 3D rendering of these velocities, normalized in each individual image with its respective average pore velocity (u avg ), are shown in Fig. 11-14.
The 3D visualisations of velocity fields show that regardless of initial pore heterogeneity, at high flowrates, corresponding to Péclet    numbers in the range of 1200-2600, acid reactant causes the growth of many channels. In contrast, at low flowrates, characterised by Péclet numbers in the range of 290-420, acid tends to dissolve the regions with less resistance to flow, resulting in preferential dissolution in fewer channels. Channel growth patterns of dissolution were also observed in scCO2/brine/calcite experiments by Menke et al. (2016). In this study, by pre-screening the dolomite cores for their initial flow field heterogeneity and classifying them accordingly, it is possible to discriminate the impact of initial pore structure from that of flowrate and observe different channel growth patterns for different initial heterogeneities and flowrates. In Fig. 13 for sample B in the low flowrate experiment, at times less than t = 100 min there are two equivalent channels growing at approximately the same rate. However, after 100 min, only one channel   Y. Al-Khulaifi et al. International Journal of Greenhouse Gas Control 68 (2018) 99-111 continues to grow. All the dissolution occurs predominantly in that channel alone, implying considerable mass transfer limitations to acid reactant that result in a slightly lower effective reaction rate overall (see Fig. 8). In contrast, the velocity images for sample B in the high rate experiment, visualised in Fig. 14, demonstrate that after 100 min, many dissolution-induced channels develop without favouring a particular one. This implies fewer mass transfer limitations and results in an approximately 8 times higher effective reaction rate, see Fig. 8. Visualisation of flow fields in Fig. 12 for sample A at a high rate indicates the similar development of many channels-the only difference is that for sample B it takes longer for the channels to develop at early times, as seen in Fig. 14. Finally, Fig. 11 shows the evolution of the flow field of sample A with low rate where preferential dissolution of a few channels is observed. This indicates that there is more contact of reactant with the solid surface, leading to a higher overall reaction rate when compared to the relatively more heterogeneous sample B at the same low flowrate. However, the reaction rate is still considerably lower than at the higher Pe experiments (Fig. 8).   Putting initial rock heterogeneity aside and strictly looking at the global effect of flowrate, the samples undergoing dissolution at the higher rate are associated with a narrowing of the velocity distribution over time, ultimately characterized by a single peak. This occurs in conjunction with a significant reduction in stagnant velocities, especially for dissolution in the more heterogeneous sample B. Here dissolution is more uniform with fewer mass transfer limitations. On the other hand, the velocity distributions for the lower flowrate experiments are both characterised by widening of the velocity distributions with time even though they have different initial distributions. This is a consequence of mass transfer limitations causing preferential dissolution along the single or few channels, Figs. 11 and 13, resulting in the velocity distributions characterized by more stagnant voxels compared to the initial velocity distribution.

Discussion
In a closer examination of the evolution of velocity distributions within the high rate experiments, differences are found between the responses of the two heterogeneities. The distribution of more heterogeneous sample B has a higher and narrower peak with a sharper drop towards the lower velocities, Fig. 15d. This effect resulted in a higher reaction rate than in sample A, Fig. 8, where the distribution was not as high and narrow, Fig. 15c. Furthermore, there are differences in the responses of the velocity distributions of the two heterogeneities at the lower flowrate. In comparing the two, Fig. 15a and b, the more heterogeneous sample B sees a greater broadening of the distribution than A causing reaction rates to be lower, Fig. 8. This effect is exactly opposite to how the two heterogeneities respond to high flowrates.
Tying the analyses of velocity distributions, flow fields and permeability may be used to explain the differences in effective reaction rate as a consequence of differential initial heterogeneity and flowrates for the four dolomite samples. In addition to the effective reaction rates, Fig. 16 shows PDFs of velocity, the corresponding 3D velocity fields, and permeabilities at the start (K initial ) and at the end (K final ) of the experiments for each sample.
For the low flowrate dissolution (Péclet number in the range of 300-420), the more heterogeneous sample, B, has a lower reaction rate, as the flow is confined to a narrower subset of the pore space: this is seen as a less pronounced peak in higher velocities and the 3D visualisation of the velocity field shows flow in two channels which correspond to the two peaks at high velocities in its initial PDF of velocity. This initial heterogeneity leads to only a single channel of, initially of least resistance, to be preferentially dissolved, Fig. 13. For sample A, the average effective reaction rate is higher, because there are no obvious dominating flow channels: here the initial PDF velocity plot is narrower leading to a more even dissolution. The permeabilities at the end of experiments, K final , are profoundly different (5.38 × 10 −12 m 2 versus 2.42 × 10 −11 m 2 for samples B and A respectively) indicating an appreciable impact of initial rock heterogeneity on permeability when dissolution occurs at low flowrates.
For high flowrates, the initial pore structure and its resulting PDF of velocity is not a dominant factor in creating dissolution patterns; the images show that both samples eventually undergo dissolution in many channels. The permeabilities at the end of the experiments are similar for both heterogeneities. The higher average effective reaction rate for the more heterogeneous sample B is driven by its lower initial porosity which leads to a higher pore velocity (higher Péclet number) and a more efficient penetration of reactant through the sample.
Finally, the analysis reveals an interesting difference in how velocity distributions for the low and high flowrate experiments evolve as a consequence of flowrate, Fig. 16. While the high flow dissolution results in a final velocity distribution that is more homogeneous than initially, low rate dissolution decreases the peaks of higher velocities and increases the number of stagnant velocities, causing the distribution to broaden. This is a consequence of two fundamentally different impacts of transport and dissolution. At high rates, where there is an effective penetration of reactant in the pore space, the rock dissolves rapidly and relatively uniformly. On the other hand, low flowrates, where reaction is faster in comparison to advection, cause the reactants to prefer flowing through the more permeable pathways causing these fast-flow channels to grow, increasing the flow field heterogeneity and lowering the overall reaction rates. This observation has been echoed in the work of Menke et al. (2016) during their experimentation with structurally homogeneous and heterogeneous limestone samples under similar reservoir conditions. In addition, Luquot and Gouze (2009) saw, in a heterogeneous limestone sample, a difference in dissolution regimes (using a Darcy velocity ∼ 3 × 10 −4 m/s similar to the high rate in this Fig. 15. Changes in the PDFs of velocity over the experimental time for both heterogeneities A and B flooded at low and high rates. Low flowrate distributions are shown in plots a and b, and high flowrate are shown in plots c and d. The red trends in each plot represent the initial velocity distribution of each sample, see Fig. 4. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) paper), which depended on the reactivity of the brine solution. A nonuniform regime came from equilibrating brine with pCO 2 = 10 MPa (characterized as transport controlled dissolution) and uniform dissolution from pCO 2 = 2.5 MPa reaction controlled. The conclusion there was that this observation cannot be used to model reservoir-scale processes until both the porosity changes − and the dependent permeability changes − are parameterized according to the fluid composition.

Conclusions
Dynamic micro-CT imaging is used to study the dissolution of Silurian dolomite by scCO 2 saturated brine at reservoir conditions at the millimetre scale. We investigated the effect of initial rock heterogeneity and flowrate on effective reaction rates, porosity-permeability relationship, and the nature of dissolution patterns by employing combined pore-scale imaging and direct numerical simulation modelling of flow on images. To this end, four experiments were designed for two types of dolomite rock heterogeneities, as characterized by their initial velocity distributions (samples A are more homogeneous than samples B), at two flowrates which create profoundly different transport conditions.
Firstly, it is important to emphasize that the effective reaction rates, which include mass transfer effects, are between 5-38 times lower than the batch reaction rate obtained under the same conditions of temperature and pressure but where mass transfer limitations are not present (Peng et al., 2016). Secondly, the effective reaction rates increase significantly with the flowrate, while being approximately constant over the period of the experiments. When dissolution occurs at high rates (Péclet numbers in the range of 2600-1200) the permeabilities at the end of experiments are similar for both types of heterogeneity. Here the reactant penetrates most of the pore space, widening many channels throughout the samples. However, at low flowrates (Péclet number in the range of 420-300) the reaction rates are lower, and the final permeabilities are very different. The more heterogeneous initial velocity distribution is characterized by more stagnant velocities and a lower peak of high velocities which can considerably impede transport of fresh acid reactant to the solid surface, leading to a lower reaction rate and a less significant increase in permeability.
Finally, it is observed that velocity distributions change differently from their initial states for each sample. Dissolution at the lower flowrates leads to the creation of wider velocity distributions with fewer velocities similar to the average velocity, due to the widening of a few preferential fast-flow channels. In contrast, dissolution at the higher rates results in more homogeneous velocity distributions, as dissolution is more uniform occurring over many channels. The more heterogeneous sample B exhibited a lower effective reaction rate than the less heterogeneous sample A due to the heterogeneity limiting the dissolution preferentially to one channel.
Permeability-porosity relationship parameterization must account for rock type dependent properties such as mineralogy and pore structure heterogeneity, but also dissolution regime triggered by Péclet and Damköhler conditions. Overall, this combined approach, which included analysis of velocity distribution, supported the experimental findings on the nature of dissolution regimes. Ultimately, these types of experiments are important to develop a deeper comprehension on mineral dissolution and may serve as validation for CO 2 storage predictive modelling tools which take into account geochemistry and other complexities found in carbonate rock formations.