Dynamic Three-Dimensional Pore-Scale Imaging of Reaction in a Carbonate at Reservoir Conditions

Quantifying CO2 transport and average effective reaction rates in the subsurface is essential to assess the risks associated with underground carbon capture and storage. We use X-ray microtomography to investigate dynamic pore structure evolution in situ at temperatures and pressures representative of underground reservoirs and aquifers. A 4 mm diameter Ketton carbonate core is injected with CO2-saturated brine at 50 °C and 10 MPa while tomographic images are taken at 15 min intervals with a 3.8 μm spatial resolution over a period of 2/2 h. An approximate doubling of porosity with only a 3.6% increase in surface area to volume ratio is measured from the images. Pore-scale direct simulation and network modeling on the images quantify an order of magnitude increase in permeability and an appreciable alteration of the velocity field. We study the uniform reaction regime, with dissolution throughout the core. However, at the pore scale, we see variations in the degree of dissolution with an overall reaction rate which is approximately 14 times lower than estimated from batch measurements. This work implies that in heterogeneous rocks, pore-scale transport of reactants limits dissolution and can reduce the average effective reaction rate by an order of magnitude.


■ INTRODUCTION
Reactions at the fluid-rock interface play a large role in many of the processes associated with fluid flow in porous media, including diagenesis, contaminant transport, and well acidization. 1−3 The focus of this paper is on reaction in carbon capture and storage (CCS), where a major concern is long-term storage security. 4,5 The injected carbon dioxide, CO 2 , will dissolve in the host brine, forming an acidic solution. 6−8 This acid in turn may react with the host rock, causing dissolution, particularly in carbonates, and−over longer time-scales− precipitation of carbonate. 9 Dissolution may compromise the integrity of geologic seals, allowing CO 2 to escape to the surface, while precipitation leads to long-term storage security. 10 Understanding dissolution in the rock-brine system, and the magnitude and distribution of fluid movement within the reservoir, is therefore imperative for accurate predictive modeling. 11−15 However, the rate and nature of carbonate dissolution are dependent on both the intrinsic properties of the rock such as porosity, permeability, and pore-size distribution, 16 as well as the properties of the brine. 17 −20 In particular, the dissolution rates are strong functions of fluid temperature and pressure, making it crucial to develop accurate experimental techniques capable of measuring complex timedependent physical and chemical phenomena at representative reservoir conditions.
One key observation of reactive transport experiments is that the average effective reaction rates measured in the field are typically several orders of magnitude lower than implied from batch reactor laboratory measurements. 21,22 Many reasons have been advanced to explain this, including mineral heterogeneity, weathering, and incomplete mixing in a heterogeneous flow field. However, without direct observation of pore-scale reaction it is difficult to assess which of these effects is most significant. Thus, time-resolved experiments are needed to provide both the insights into dynamic interplay between transport and reaction and to validate reactive transport models.
X-ray microtomography (μ-CT) is now an established experimental method for studying pore-scale processes for both petroleum and carbon storage applications. 23−25 μ-CT offers several advantages: it is noninvasive, achieves high spatial resolutions of down to around 1 μm, gives good contrast between fluid and solid phases, and provides three-dimensional images relatively quickly.
Imaging of dissolution in limestone rock has been performed at the core (∼cm) scale 26 where, for several rock types and flow regimes, it was found that reaction of the dissolved CO 2 with the solid matrix increases physical heterogeneity. Ellis et al. 27 used core scale imaging to show that exposure to CO 2 progressively reduced the sealing capacity of fractured claystone caprocks. Smith et al. 28 and Hao et al., 29 combined μ-CT and modeling of samples from the Weyburn Oil field in Saskatchewan, Canada. Multiple dolomitic cores were imaged before and after core flooding with CO 2 -saturated brine at 60°C and 12.4 MPa. The results were subsequently used to validate a continuum-scale reactive transport model. Additionally, Luquot and Gouze 20 and Gouze et al. 30 performed flow experiments with a CO 2 -rich brine on limestone reservoir samples at 100°C and 12 MPa. Samples were analyzed with μ-CT prior to and after dissolution and changes in the permeability, porosity, tortuosity, and hydraulic radius were characterized according to dissolution regime. Noiriel et al. 31 used μ-CT to measure the three-dimensional changes in a limestone fracture at several points during acid dissolution at ambient conditions. Although these experiments have provided valuable information on pore-space geometry and flow, such as the porosity-permeability relationship, the analysis was limited to comparisons between pre-and post-reaction images of the rocks studied, some were not conducted at representative subsurface conditions, while others did not have the resolution required for pore-scale predictive modeling. Measuring dynamic reaction-induced changes in pore-space geometry, topology and flow in subsurface rock systems at reservoir temperatures and pressures furthers our understanding of how different transport and reaction conditions alter the complex solid and pore structures, and is the main goal of our study. Moreover, our aim is to provide a series of time-evolving X-ray image based data sets that can serve as a benchmark experimental study to calibrate pore-scale reactive transport models.
Studying in situ reactive transport at reservoir conditions is challenging. Andrew et al. 32,33 have designed an experimental apparatus to image pore-scale multiphase displacement for the determination of residual saturation, contact angle, and interfacial curvature at elevated temperatures and pressures. These studies were performed at low flow rates and in the absence of chemical reaction to isolate the impact of capillary forces. Since CO 2 dissolves in water creating acid, it is also of interest to examine the dynamic nature of chemical reaction between the acid and the host rock at elevated temperatures and pressures. 34 Therefore, in this study we examine reactive dissolution processes and focus on measuring the time-dependent reaction rate between the CO 2 -acidified brine and limestone rock. We select this fluid/solid system since it is important in several environmental applications including subsurface contaminant transport, geological CO 2 storage, well acidization, and ocean acidification. 35−37 We perform dynamic in situ imaging of reaction between a CO 2 -acidified brine and Ketton limestone at the pore scale at an initially fixed Pećlet, Pe, and Damkoḧler, Da, numbers over a period of several hours and at a temperature and pressure representative of an aquifer at approximately 1 km depth. We then use pore-scale modeling methods on μ-CT images to characterize the evolving geometrical, topological and flow features, providing physical insights into the dynamics of dissolution, and specifically how to relate average effective reaction rates to pore-scale flow and transport.

■ MATERIALS AND METHODS
Experimental Apparatus. Our experimental apparatus, see Figure 1, is an adaption of that presented in Andrew et al. 32 A brine solution representative of a typical saline aquifer consisting of 1% KCl and 5% NaCl by weight was preequilibrated with scCO 2 in a Hastelloy reactor (Parr instruments Co., IL) at 10 MPa and 50°C and mixed using an entrainment stirrer. The core was confined using deionized water at 13 MPa. A separate, highly X-ray absorbing brine of 6% KI by weight was loaded into the receiving pump and flowed backward through the sample until the core was thoroughly saturated. Excess KI brine was ejected from the system before reaching the reactor to avoid contamination.
The Hassler cell consisted of a carbon fiber sleeve with steel end-caps. The carbon fiber sleeve allowed for X-ray transparency during the experiment. Heating tape was wrapped Figure 1. The experimental apparatus consists of three pumps: the injection pump controlling pressure; the receiving pump controlling flow−with a Hassler cell and reactor in between; and a confining pump. The core is wrapped in aluminum foil and then placed into a Viton sleeve and attached to flow lines within the Hassler cell. The receiving pump pulls reactor-equilibrated brine through the core while the injection pump maintains pressure. Additionally, the third pump confines the core using a radial pressure boundary condition on a Viton confining sleeve within the carbon fiber sleeve of the Hassler cell. All components are outside the radiation-controlled zone of μ-CT scanner except the Hassler cell, which is mounted vertically on the rotation stage within the scanner. around the sleeve for temperature control. Within the cell, a 4 mm diameter, 1.2 cm long Ketton core was placed inside a Viton sleeve connected to flow lines. A thermocouple was inserted into a radial port of the cell and adhered to the outside of the Viton sleeve with aluminum tape. The thermocouple in combination with the heating tape allowed for accurate control of the sample temperature, while still permitting the X-ray transparency needed for in situ imaging.
Reactor fluidthe 5% NaCl, 1% KCl brinewas then pulled through the core by the receiving pump for 2.5 h at a constant rate of 0.5 mL/min throughout the experiment, representing a Darcy velocity of 6.6 × 10 −4 m.s −1 , with the injection pump acting as a front pressure regulator at 10 MPa. The significant difference in attenuation values between the two brines allowed visual confirmation using two-dimensional image projections of the exact time that acidic brine started to flood the core, making cumbersome dead-volume calculations unnecessary.
Experimental Method. We imaged dynamic changes in the fluid/solid interface during reaction between unbuffered CO 2 -saturated brine and Ketton carbonate using the Versa XRM-500 X-ray Microscope (www.zeiss.com). Ketton is a 99.1% calcite 33 Limestone Oolite from the Upper Lincolnshire Limestone Member that was deposited approximately 169−176 million years ago. Mercury intrusion experiments on Ketton show a bimodal pore structure that allows for clear resolution of the macropores by μ-CT, while leaving the micropores unresolved but able to be quantified by analysis of the change in absolute gray scale. 38 The images were taken at a 3.8 μm spatial resolution with 400 projections. A 4 mm vertical section of the core was imaged, starting 2 mm away from the inlet. This location was chosen to avoid boundary effects. This 4 mm section will hereafter be referred to as the sample. The pressure drop across the sample was not measured because the sample field of view was much smaller than the overall core, and the pressure drops in any event were very small. In complex carbonate pore structures there is always a trade-off between the sample size and resolution, and we have designed the experiment to maximize the reactive transport information available from the experiment at the μm to mm scale. Images were taken successively as fast as our scanner allowed, with ∼15 min between the start of each scan ( Figure 2). The images are thus the time-averaged change over a single scanning period. We then reconstructed, filtered, and segmented the images to generate a binarized representation of the pore space. where U av [m.s −1 ] is the average pore velocity (Darcy velocity q divided by the scan porosity ϕ CT , measured directly from the segmented image−its initial value was 0.17, rising to 0.32 by the end of the experiment), L [m] is the characteristic length and D is the molecular diffusion coefficient of Ca 2+ at 25°C (7.5 × 10 −10 m 2 .s −1 ). 20 Surface area A S [m 2 ] was computed on the images by adding the number of voxel faces separating void from grain. We use the definition of characteristic length L = π/ S from Mostaghimi et al. 39 where the specific surface area [m −1 ] S = A S /V B and V B is the bulk (total) volume of the sample. S is 7.72 × 10 3 m −1 and L is 4.07 × 10 −4 m at the start of the experiment (Table 1).
Lai et al. 40 measured a BET surface area (A BET ) value of 2.5 × 10 3 m 2 .kg −1 on a similar Ketton core before reaction. Using Andrew et al.'s 33 determination of the intragranular porosity of Ketton grains, ϕ grain = 0.12, we calculate the BET specific surface area to be 5.1 × 10 6 m −1 using S BET = ρ calcite [1 − ϕ total ] A BET where the total porosity ϕ total = ϕ grain [1 − ϕ CT ] + ϕ CT = 0.27 and ρ calcite = 2.71 × 10 3 kg·m −3 . The specific surface area, measured at the resolution of the μ-CT scan is approximately 650 times smaller than the molecular-level BET measured value; the large difference is due to the inclusion of the area of the microporosity in the BET measurement and because the μ-CT cannot fully resolve surface roughness down to the molecular scale.
Using the surface area from the μ-CT scan, the characteristic length L = 407 μm with Pe = 2100 at the start of the experiment; at the end L = 392 μm and Pe = 1100. The high Pe number means that, on average, advection occurs much faster where k is a reaction rate constant measured with units of 1/ time. Physically this represents the ratio of the time for advection to the time for reaction over a characteristic length L.
We interpret the reaction time to be the time taken to dissolve a section of the rock of length L. The dissolution rate for a flat, pure sample of calcite, with no transport limitations, measured by Peng et al. 42 at the same experimental conditions as here, is r = 6.9 × 10 −4 mol.m −2 s −1 . If we assume that Ketton is pure calcite then an estimate of the reaction rate constant is k = rS/n = πr/nL, where n is the number of moles of calcite per unit volume of Ketton. n is calculated using n = ρ calcite [1 − ϕ total ]/ M calcite where M calcite is the molecular mass of calcite (0.1 kg mol −1 ). We find n = 2.0 × 10 4 mol.m −3 at the start of the experiment. We can therefore write eq 2 as Our estimated Da increases from approximately 2.8 × 10 −5 to 6.5 × 10 −5 during the experiment (U av decreases as the porosity increases for a fixed Darcy flow rate); however, it is always much less than 1. This represents the conditions at which reaction rate is slow compared to the advection rate and is thus representative of the uniform dissolution regime. 26 We only see significant dissolution once several hundreds of pore volumes have been injected, consistent with the low value of Damkoḧler number. In total, ∼2000 pore volumes were injected through the core in our experiment.
Over the scale of a characteristic length, the advective time is smaller than the diffusive time (Pe ≫ 1), while the ratio of reaction rate to diffusion rate PeDa increases slightly from 5.9 × 10 −2 to 7.0 × 10 −2 during the experiments. This implies that at the pore scale, reaction is still slow compared to diffusion. However, as we discuss later, diffusion may limit the average effective reaction rate through delaying the transport of reactants to and from the solid surface.
The estimated brine pH of 3.1 is computed using the chemical equilibrium calculation methods in Leal et al. 43,44 that has been experimentally confirmed by Peng et al. 8 Modeling Methods used in the Analysis of the Experiments. We have used network extraction for characterizing pore topology, and numerical simulation of flow directly on the image voxels. Our network extraction method is based  . Surface area to the bulk volume ratio (squares) and the imaged-averaged porosity (crosses) as a function of time measured on 10 timeevolving μ-CT images (a). Image 7 (at 100 min) was of poor quality due to X-ray source fluctuations, which may explain why the measured specific surface area is lower than the general trend. Permeability plotted as a function of porosity through time (b). A power-law relationship with the equation K = 1.4 × 10 −7 ϕ 5.16 is seen between porosity and permeability as dissolution progresses.  Figure 3b shows corresponding change in porosity for each slice. Porosity is determined by the ratio of the number of "pore voxels" to total voxels in the segmented slices. Figure 3b shows that a change in porosity occurs throughout the sample, which eliminates the possibility of face dissolution as the dominant dissolution regime for which the change would occur only at the inlet. The imaged-based average porosity increases from 17 to 32%, as seen in Figure 3a; we see no significant difference in dissolution rate over the length of the sample.

Environmental Science & Technology
We initially see a rapid increase in porosity, which slows with time; similarly we see a fast increase in the surface area to volume ratio, followed by a more gradual decline at later times, Figure 4a. It is possible that this could be due to the preferential dissolution of cement between grains, rather than the grains themselves (see Figure 5) or simply because dissolution forms new flow paths, resulting in an increasingly heterogeneous flow field, with diffusive transport further limiting the average effective reaction rate.
To derive average effective reaction rate, we assume that chemical conditions are relatively constant: at the high Pe in our experiment, the acid flows fast throughout the sample, so that fresh reactant is continuously provided. We can then estimate the reaction rate from the change in porosity during the experiment. We calculate an average effective reaction rate: where Δt is the time between scans [s], Δϕ CT is the corresponding change in porosity and S is the image surface area per unit volume at the beginning of the time period [m −1 ]. The term 1 − ϕ grain = 0.88 accounts for the fact that the grains themselves have microporosity, as discussed earlier. Initially r eff = 8.8 × 10 −5 mol.m −2 s −1 (measured between scans 1 and 2), which decreases to 2.3 × 10 −5 mol.m −2 s −1 (measured between scans 9 and 10) by the end of the experiment. The average reaction rate over the life of the experiment is 5.0 × 10 −5 mol.m −2 s −1 . This reaction rate is approximately 14 times lower than the laboratory-measured batch reaction rate of 6.9 × 10 −4 mol.m −2 s −1 . Had instead we used the BET-measured surface area to define r eff , our average effective reaction rate would be around 4 orders of magnitude lower than the batch reactor value. It is evident from the constant gray scale in the grains (even as they shrink) that there is no discernible reaction in the microporosity and so it is inappropriate to use the BET area. However, even at a pore-scale resolution of 3.8 μm, the effective surface area for reaction is considerably lower than that measured using the μ-CT scan. As discussed further below in the pore scale analysis section, we suggest that this is caused by diffusion-limited transport of reactants to and away from the surface in relatively stagnant regions of the pore space where the local Pećlet number is low. It should be noted that, although we did not analyze the effluent, some fine particle movement could be present 49,50 that would further decrease the measured reaction rate. The grain size distribution over the course of the experiment is studied to test how dissolution affects characteristic length− this is presented in Figure 5. The grain diameter is calculated by performing a watershed segmentation of the three-dimensional distance map of the pore space using the separate object module in Aviso 8. Features with a diameter smaller than 100 μm are not included in the analysis because this method of segmentation is less accurate with very small objects, which are unlikely to be whole grains of Ketton. We see that for large grains the average grain diameter decreases in size from 650 to 500 μm [ Figure 5a] implying dissolution. Note that the grain diameter is similar to the characteristic length used to define Pe, which decreases from around 410 to 390 μm during the experiment. Additionally, we see that there is no discernible difference in the change in grain diameter over the length of the sample (each section changes by approximately 10%) confirming that−at the core scale−dissolution is uniform [ Figure 5b-d].
Pore-Scale Analysis. Figure 6 depicts the pore-networks extracted from the images acquired at time 0 (a), 67 (b), and 149 (c) minutes. The pore size and throat diameter all increase with time, while the total number of pores and throats decreases ( Table 1). The coordination number (representing the connectivity) increases slightly.
We simulated flow directly on the voxels of the segmented images to quantify changes in porosity, permeability, and the flow field. From Figure 4b we see an increase in porosity and permeability at each measured time as the reaction progresses. Permeability increases by almost 1.5 orders of magnitude due to dissolution, which makes the pores larger and better connected ( Table 1). The relationship between permeability and porosity can be fitted to a power law with an equation of K = 1.4 × 10 −7 ϕ 5. 16 . We observe that the permeability increases with no stepwise jumps that would be expected in the wormhole regime at the times when a wormhole breaks through the outlet; 51 this, together with the earlier observation eliminating face dissolution (which would give more dissolution at the entrance of the sample), leads us to infer that we observe the uniform dissolution regime.
In Figure 7 we assess the changes in velocity in fast flow channels by plotting the velocity fields in the pore voxels where advection is dominant for times 0 (a) and 149 (b) minutes. It can be observed in Figure 7 that, as the reaction evolves, the fast flow paths become less tortuous as a consequence of consumption of the spherically shaped grains of Ketton and surrounding cement. The flow also tends to become concentrated through the center of the core; this may be explained by the nonreactivity of the confining Viton sleeve driving flow toward the center as the carbonate reacts and forms larger flow paths. However, Figure 2 shows that there is still dissolution near the boundaries of the core.
The distribution of flow speed for the first and last acquired image is shown in Figure 7c, as in the analysis by Bijeljic et al. 46,52 in which the ratios of the magnitude of u at the voxel centers divided by the average pore velocity U av = q/ϕ are presented. Note that the normalized average pore velocity decreases with time, since the porosity increases.
As suggested by Figure 7a,b, the fast flow paths become more uniform; this is demonstrated in Figure 7c, which shows that the spread in velocity distribution for the fastest voxels (|u|/U av > 1) becomes narrower over time. However, the more significant finding is that the focusing of the flow paths leaves an increasing fraction of the pore space with much lower flow speeds: the tail of the velocity distribution becomes more spread. In these regions, at the micron scale, diffusion may limit the transport of reactants to and from the solid surface. As a consequence, the average effective reaction rate is lower than that seen in the batch reactor, with advectively well-mixed fluids. Furthermore, the spatial distribution of stagnant regions−for instance neighboring fast channels−may be such to suppress transport to and from other portions of the pore space with faster flow speeds. As discussed above (see Table 1), this is a significant effect, with a decrease in the reaction rate by a factor of 3.8 during the experiment and an overall rate approximately 14 times lower than measured in the batch reactor. Explicit reactive transport modeling in the pore space, which is beyond the scope of this work, is necessary to quantify the impact of transport limitations on reaction and to determine if this is indeed the reason for the low average effective reaction rates encountered in this system. 53,54 Our approach allows us to integrate experimental in situ reactive transport measurements with corresponding modeling analysis on a voxel by voxel basis for supercritical CO 2 -brinerock systems at elevated temperatures and pressures. We find an increase in porosity and permeability due to fluid/solid reaction that is more rapid at the initial stage of the experiment. Permeability increases with increasing porosity as a power law,  with a flow field that shows a decrease in tortuosity and flow focusing of the fast flow channels. Analysis of the reactioninduced changes in the experimental porosity profiles and grain distribution indicate that at the flow conditions studied, Ketton dissolved predominantly in the uniform dissolution regime. Detailed analysis of the flow fields revealed that, as the fluid/ solid reaction progresses, the flow becomes more focused in fast regions; a large fraction of the pore space then experiences lower flow speeds. We suggest that this leads to diffusionlimited transport of reactants to the surface and explains the 14fold decrease in the average effective reaction rate compared to measurements in the batch reactor at the same conditions. This in situ experimental method for dynamic reservoir condition μ-CT imaging could be extended to study the impact of heterogeneity, flow rate, and reaction rate on the mechanisms controlling dissolution. The applications of this technique are far reaching in the petroleum, CCS, and subsurface hydrology areas including subsurface contaminant transport, 35 ocean acidification, 36 and well acidization. 37 Moreover, the results from this study can serve as a benchmark for fluid/solid reaction (similar to Gramling et al. 55 for fluid/ fluid reaction) in the verification of pore-scale reactive transport models 53 and elucidate potential scaling factors for resolving the discrepancy between laboratory and field reaction rates in continuum scale simulations.

* S Supporting Information
Further information on image acquisition and processing and modeling methods and segmented μ-CT scans at 10 time steps. This material is available free of charge via the Internet at http://pubs.acs.org.

Author Contributions
The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
We would like to thank the Qatar Petroleum, Qatar Science and Technology Park, and Shell as sponsors of the Qatar Carbonates and Carbon Storage Research Centre (QCCSRC) for generously funding this research.

■ ABBREVIATIONS
CCS carbon capture and storage CO 2 carbon dioxide Da Damkoḧler number μ-CT X-ray microtomography PDF probability density function PéPećlet number sc supercritical