Reservoir condition imaging of reactive transport in heterogeneous carbonates using fast synchrotron tomography — Effect of initial pore structure and ﬂ ow conditions

Weinvestigatetheimpact ofinitial porestructureandvelocity ﬁ eldheterogeneity on thedynamicsof ﬂ uid/solid reaction at high Péclet numbers (fast ﬂ ow) and low Damköhler number (relatively slow reaction rates). The DiamondLightsourcePinkBeamwasusedtoimagedissolutionofEstailladesandPortlandlimestonesinthepres- enceofCO 2 -saturatedbrine atreservoirconditions(10 MPa and 50°C representing~1 kmaquifer depth) attwo ﬂ ow rates for a period of 2 h. Each sample was scanned between 51 and 94 times at 4.76- μ m resolution and the dynamic changes in porosity, permeability, and reaction rate were examined using image analysis and ﬂ ow modelling. We ﬁ ndthatthe porosity caneitherincrease uniformly throughtimealong thelength ofthesamples,or may ex- hibit a spatially and temporally varying increase that is attributed to channel formation, a processthat is distinct fromwormholing,dependingoninitialporestructureand ﬂ owconditions.Thedissolutionregimewasstructure-dependent: Estaillades with a higher porosity showed more uniform dissolution, while the lower porosity Portland experienced channel formation. The effective reaction rates were up to two orders of magnitude lower than those measured on a ﬂ at substrate with no transport limitations, indicating that the transport of reactant and product is severely hampered away from fast ﬂ ow channels. is an open the BY 1 m, and reaction conditions alter complex solid and pore structures, it is necessary to measure reaction-induced changes in pore-space geometry, topology and ﬂ ow in subsurface rock systems at reservoir temperatures and pressures. High resolution is required to investigate pore-scale processes in detail. This paper describes a method of studying dynamic reactive dissolution processes in rock with complex pore structures with a focus on measuring the time and spatially dependent reaction rate between a CO 2 -acidi ﬁ ed brine and limestone rock at reservoir conditions. We observe a distinct dissolution regime, channelling, whose emergence is dependent on initial pore structure and ﬂ ow conditions. Several studies have investigated reaction in complex carbonates at due to experimentalor imagingcon-straints they have been either limited to pre- and post-reaction images, not completed at representative subsurface performed dynamic in situ imaging of reaction between a CO 2 -acidi ﬁ ed brine and Ketton limestone at the pore scale over a period several hours and at a temperature and pressure representative of an aquiferatapproximately1kmdepth.However,Kettonisarelativelyho-mogenous carbonate rock is easy to image very little time traditional μ -CT methods. Many carbonate rocks have complex pore structures that more detail to accurately resolve, can be a very time intensive process using traditional μ -CT either with a monochromatic beam at a synchrotron source, bench-top x-ray scanners. a fast method of tomography is required to see dynamic reaction-induced changes in heterogeneous carbonates. beenpreviously used by to study image Haines jumps in two-phase ﬂ ow, to study the dynamic capillary pressure changes during a reservoir condition CO 2 -brine drainage on the pore scale. We present a method thatusesfastsynchrotrontomographytostudytheeffectofini-tialpore structure and ﬂ ow conditionsonthedynamics of ﬂ uid/solid reaction in two highly heterogeneous rocks at reservoir conditions. The dynamic image sequence is ﬁ rst analysed for porosity changes and then used as inputs into ﬂ ow simulations that can elucidate changes in permeability and pore-scale velocity ﬁ effective reaction rates ow compared termsof demonstratingboth


Introduction
A major concern of carbon capture and storage (CCS) is long-term storage security (Herzog et al., 2003;Metz et al., 2005). Carbon dioxide, CO 2 , injected into the subsurface will dissolve in the host brine and form carbonic acid (Langmuir et al., 1997;Morse and Mackenzie, 1990;Peng et al., 2013). Carbonate host rocks have the potential to react with and be dissolved by CO 2 acidified brine (Peng et al., 2015). Dissolution of the host rock may instigate unpredictable fluid flow and can weaken carbonate cements and damage injection wells (Birkholzer et al., 2011;Mohamed and Nasr-El-Din, 2012). A complete understanding of dissolution in the brine-rock system is therefore important to predict the distribution and the rate of fluid movement and the amount and impact of dissolution in the subsurface (Daccord et al., 1993a, Daccord et al., 1993bMaheshwari et al., 2013).
However, the nature and the rate of dissolution in carbonates are dependent on both the properties of the brine (El-Maghraby et al., 2012;Fredd and Fogler, 1998;Luquot and Gouze, 2009) and the host rock (Li et al., 2006;Rötting et al., 2015). Carbonate dissolution rates are also strongly dependent on brine temperature and pressure (Peng et al., 2015), making it necessary to develop experimental techniques or measuring complex time-dependent processes at representative reservoir conditions.
Previous experiments have observed that field-scale reaction rates are typically orders of magnitude lower than experimental batch reactor measurements (Swoboda-Colberg and Drever, 1993). Mineral heterogeneity and incomplete mixing in a heterogeneous flow field are potential explanations for this phenomenon. However, it is not possible to assess the most significant factors without directly observing the evolution of the pore space during reaction. Thus, dynamic pore-scale experiments are required to provide both insights into the interplay between transport and reaction, and to validate predictive models. Traditionally, dissolution has been studied in larger (cm-scale) core floods, or in idealized two-dimensional systems (Fredd and Fogler, 1998;Marini, 2007;Ott et al., 2012). However, these experiments are unable to determine the mechanisms of dissolution within the pore space in realistic rock samples.
X-ray microtomography (μ-CT) (Berg et al., 2013;Blunt et al., 2013;Gouze and Luquot, 2011;Wildenschild and Sheppard, 2013) is an established experimental method for studying pore-scale processes in CCS applications. μ-CT has several benefits: it achieves high spatial resolutions of down to around 1 μm, it is non-invasive, and it provides three-dimensional images. To advance our understanding of how different transport and reaction conditions alter complex solid and pore structures, it is necessary to measure reaction-induced changes in pore-space geometry, topology and flow in subsurface rock systems at reservoir temperatures and pressures. High resolution is required to investigate pore-scale processes in detail.
This paper describes a method of studying dynamic reactive dissolution processes in rock with complex pore structures with a focus on measuring the time and spatially dependent reaction rate between a CO 2 -acidified brine and limestone rock at reservoir conditions. We observe a distinct dissolution regime, channelling, whose emergence is dependent on initial pore structure and flow conditions.
Several studies have investigated reaction in complex carbonates at the pore scale (Gharbi, 2014;Gharbi et al., 2013;Hao et al., 2013;Noiriel et al., 2013;Smith et al., 2013), but due to experimental or imaging constraints they have been either limited to pre-and post-reaction images, or were not completed at representative subsurface conditions. Menke et al. (2015) performed dynamic in situ imaging of reaction between a CO 2 -acidified brine and Ketton limestone at the pore scale over a period of several hours and at a temperature and pressure representative of an aquifer at approximately 1 km depth. However, Ketton is a relatively homogenous carbonate rock that is easy to image in very little time (~17 min) with traditional μ-CT methods. Many carbonate rocks have complex pore structures that require more detail to accurately resolve, which can be a very time intensive process using traditional μ-CT hours), either with a monochromatic beam at a synchrotron source, or with bench-top x-ray scanners. Thus, a fast method of tomography is required to see dynamic reaction-induced changes in heterogeneous carbonates.
Fast synchrotron tomography has been previously used by Berg et al. (2013) to study image Haines jumps in two-phase flow, and by Andrew et al. (2015b) to study the dynamic capillary pressure changes during a reservoir condition CO 2 -brine drainage on the pore scale. We present a method that uses fast synchrotron tomography to study the effect of initial pore structure and flow conditions on the dynamics of fluid/solid reaction in two highly heterogeneous rocks at reservoir conditions. The dynamic image sequence is first analysed for porosity changes and then used as inputs into flow simulations that can elucidate changes in permeability and pore-scale velocity fields. Additionally, effective reaction rates are calculated and the reactive flow dynamics is compared in terms of the initial structure and flow conditions, demonstrating both uniform and channelled dissolution structures.

Methods
The photon flux and energy of the x-ray source control the image acquisition time of the sample. The use of the polychromatic beam of a synchrotron source is one method of scanning quickly (Andrew, 2015;Andrew et al., 2015b;Berg et al., 2013). This so-called 'Pink Beam' provides orders of magnitude more intense light than bench-top sources and images can therefore be taken on the tens-of-second rather than hour time scale. Using this intense x-ray flux, however, is not without its challenges. The lower energy x-rays of the Pink Beam spectrum are absorbed by the sample as heat. This can interfere with the temperature control of the in situ apparatus and cause CO 2 to exsolve from solution. CO 2 -saturated brine is very sensitive to both heat and pressure and therefore a small change in thermal equilibrium can significantly change the pH of the in situ fluid (Peng et al., 2013). Thus, careful design and control elements for the x-ray spectra must be incorporated into the beam line equipment prior to imaging. Mirrors and filters are thus used to narrow the spectrum of light to suit experimental needs and absorb unwanted x-ray spectra. Mirrors absorb the high-energy spectrum while filters absorb the lower energies. It is therefore possible to narrow the spectrum to the desired band of radiation using only these tools. The energy of the Diamond Light Source I13-2 Pink Beam ranges from 8 to 30 keV with a maximum photon flux of 4 × 10 9 Ph/s. The low energy x-rays were filtered by passing the beam through 2 mm of Al and 0.1 mm of Au. A 250 μm-thick CdWO 4 scintillator with a 1.25× objective lens and a PCO EDGE camera were used in the detector assembly (Andrew et al., 2015b).
Tomographic scans were taken successively for 2 h at the centre point of the core with a 5 mm by 4 mm field of view. The dynamic scans had 1000 projections each with an exposure time of 0.04 s. Total acquisition time was~45 s per scan with between 51 and 94 scans taken over a period of 2 h. The time between scans was controlled by the data read-off speed of the camera, which varied between experiments.

The experimental apparatus
Estaillades and Portland Basebed carbonates were reacted using our in situ experimental apparatus [ Fig. 1]. Reaction was imaged between calcite and unbuffered scCO 2 saturated brine in 4 mm-diameter 1.2 cm-long carbonate cores. Carbonate rocks generally have a wide Fig. 1. The in situ experimental apparatus. CO 2 is pressurized by the injection pump and used to equilibrate brine in the reactor. Reactive brine is pulled through core assembly by the receiving pump. The cell is confined by deionized water in the confining pump at 12 MPa and heated using heating tape controlled by a thermocouple in the confining fluid. The experimental system is connected together using tubing and fluid flow is directed using Valves (V) and Unions (U). Modified from Menke et al. (2015). range in pore structure, heterogeneous velocity fields, and significant microporosity. Portland carbonate is a limestone oolite with a wide pore size distribution characterized by a few fast flow regions and abundant regions of stagnant flow (Bijeljic et al., 2011). Estaillades is a wellconnected carbonate with a bimodal pore-size distribution (Bijeljic et al., 2013a) that has abundant microporosity. Both carbonates are relatively pure (N 97.2%) calcite with the remaining fraction being quartz (Andrew et al., 2014;Gharbi, 2014).
The fluid was pre-equilibrated in a vigorously mixed Hastelloy reactor heated and pressurized to experimental conditions (Andrew et al., 2015a, Andrew et al., 2014. Fully scCO 2 saturated 1% KCl 5% NaCl brine with an ionic strength of 1.05 M was pulled through the core at flow rates of 0.1 and 0.5 mL/min at 50°C and 10 MPa: these flow rates were maintained throughout the experiment. The confining pressure of the Hassler cell was 12 MPa. Reservoir temperature was maintained in the core holder by wrapping the exterior of the sleeve in an x-ray transparent heating tape and inserting a thermocouple through the radial port of the cell and into the confining fluid. A Proportional Integral Derivative (PID) controller then regulated temperature to within 1°C. Pressure and flow conditions were maintained using three high-pressure syringe pumps that are accurate to a flow rate of 0.001 mL.min. Two salts were used for the experiment, a highly absorbing 25% wt KI unreactive brine, and a low absorbing 1% wt KCl, 5% wt NaCl reactive brine. The difference in attenuation between the two brines made it easy to see the arrival of reactive brine in the core making dead volume calculations unnecessary. Effluent sampling was not possible during these experiments as all pumps and lines were located inside the radiation enclosure.

Image processing
The raw projections were reconstructed using the filtered back projection method (Kyrieleis et al., 2009;Kyrieleis et al., 2011;Titarenko et al., 2010) and ring artefacts were removed using the Diamond Lightsource proprietary python code. Each image consists of 2000-3 voxels, which were then binned (2 × 2 × 2) to increase signal to noise resulting in an image of 1000 3 with a size of 4.76 μm. The images were then cylindrically cropped to the core dimensions and converted to 16-bit from 32-bit to further decrease size using ImageJ (http:// imagej.nih.gov/ij/). The images were then further processed using the image processing modules in Avizo 9.1 program (www.vsg.com).
The image processing workflow is presented for Portland carbonate in Figs. 2a-f. First the raw images [Fig. 2a] are filtered using the nonlocal means edge preserving filter (Buades et al., 2005, Buades et al., 2008 [ Fig. 2b]. The first scan in the series will hereafter be called the reference scan. All scans were registered to the reference scan using normalized mutual information and resampled using the robust Lanczos method (Lanczos, 1950). Scan histograms were compared to verify that the greyscale linear attenuation coefficient did not vary between scans. The resampled data was then subtracted from the reference scan to measure the change in the pore space as dissolution takes place [ Fig. 2c]. These 'difference images' were then filtered using non-local means to increase signal to noise [ Fig. 2d]. The reference scan and filtered difference images were then segmented using a watershed segmentation (Schlüter et al., 2014) [ Fig. 2e]. Finally, the segmented difference image was subtracted from the segmented reference scan to achieve segmented images for each time step [ Fig. 2f].

Flow modelling and network extraction
Each experiment was modelled using a flow solver directly on the same central region of the segmented images with 800 3 voxels at seven times during dissolution. The flow model solves both the volume conservation equation and the Navier-Stokes equations employing the method presented in Bijeljic et al. (2013b) and Raeini et al. (2012) which is based on the finite-volume code in OpenFoam (2011). An arbitrary fixed pressure drop is imposed across the sample in the direction of flow. All other solid boundaries are no flow boundaries. The simulated flow rate is then scaled to the flow rates used in the experiments using the experimental Darcy velocity. Each image voxel represents a grid block in the flow simulator that provides velocities at the voxel centres. In addition, the images are also input into our network extraction code, based on the maximal ball algorithm (Dong and Blunt, 2009), to characterize pore space connectivity. Spheres are grown in the pore space of a segmented image to identify the location and size of the pores and pore throats, which then provide a topological representation of the pore space as wide pores connected by narrow throats. Later in the Results section we will study dynamic changes in the velocity fields by visualizing the ratios of the magnitude of u at the voxel centres divided by the average pore velocity U av = q/ϕ, where q is the Darcy velocity in m·s −1 and ϕ is the porosity. We will also study the changes in pore space connectivity as revealed by network extraction.

Results and discussion
We study impact of initial pore structure and flow conditions for two types of carbonate: Estaillades and Portland Basebed [ Table 1]. In experiments 1 and 2, we react an Estaillades core at a low (0.1 mL·min −1 ) and high (0.5 mL·min −1 ) flow rate to gain insight into the impact of flow conditions. Then in experiments 3 and 4 we use the same high and low flow rates to react cores of Portland to investigate the impact of pore space heterogeneity.
Using image analysis we first investigate changes in porosity and surface area. Then we input the images into a flow solver to study changes in permeability and velocity fields. Furthermore, we investigate the dynamic changes in the dimensionless numbers describing flow and reaction conditions and the effective reaction rates.

Porosity, surface area, and permeability
The segmented images were analysed as a time series for porosity changes by counting the number of voxels of pore and rock. During dissolution, porosity increased with time. Fig. 3 depicts porosity averaged over each two-dimensional slice of the three-dimensional image along the axis of flow for each experiment at 30-min intervals. In all four experiments porosity changed uniformly in the direction of flow. For experiment 1, the porosity also increased uniformly through time. However, in experiments 2, 3, and 4 the majority of the porosity increase occurred during the first 30 min with a slower rise at later times. Experiments 1 and 3 had a total increase in porosity of 6.6% and 5.5% respectively, which is smaller than the total porosity increase in the higher-rate experiments 2 and 4 (11.5% and 12.7%).
Permeability was calculated directly on the segmented images using the Navier-Stokes solver that computes voxel centre velocities. The permeability was found from the ratio between the total flow rate across the whole sample to the average pressure gradient multiplied by fluid viscosity, from Darcy's law. Fig. 4 depicts the dynamic changes in porosity and permeability for Estaillades and Portland carbonate images at flow rates of 0.1 mL·min − 1 and 0.5 mL·min − 1 . The change in total image porosity is plotted as a function of time in Fig. 4a. Experiment 1 had a steadily increasing porosity throughout the experiment. Experiments 2, 3, and 4 increased in porosity quickly at first and then more slowly, as evident from Fig. 3.
Permeability increases rapidly as a function of porosity [ Fig. 4b]. The greatest change in permeability with porosity was in experiment 4, which had the highest flow rate and the lowest initial permeability [Fig. 4c]. The low flow rate experiments 1 and 3 showed the least change in permeability with porosity, while experiment 2, which had a very high initial permeability had large changes in porosity with a smaller permeability increase. We fit the data to an empirical power-law relationship between porosity and permeability. The best fit to the Estaillades experiments 1 and 2 were K = 7 × 10 − 8 ϕ 6.74 and K = 1 × 10 −6 ϕ 7.52 respectively. For the Portland experiments 3 and 4 we found K = 5 × 10 −4 ϕ 10.9 and K = 1.9 × 10 −3 ϕ 11.0 respectively. The exponents, around 7 to 11, are larger than the value of approximately 5 observed for dissolution of a homogeneous limestone, Ketton [Menke et al. (2015)] and much higher than a power of approximately 3 used in the Kozeny-Carman relation to predict the permeability of granular packings. This suggests that we see a very rapid rise in permeability with porosity for more heterogeneous samples.
The specific surface area, S, [ Fig. 4d], was computed on the images by adding the number of voxel faces separating grain from void. S is the surface area A s [m 2 ] divided by the bulk (total) volume. S showed an approximately linear increase with time for experiment 1, and a sharp initial increase for experiments 2, 3, and 4, with less change at later times. There was no clear trend with flow rate for experiments 3 and 4.

Pore-scale velocity fields
In Fig. 5 we visualize the ratios of the magnitude of u at the voxel centres divided by the average pore velocity for Portland and Estaillades carbonate images. U av [m·s −1 ] is calculated from the Darcy velocity q divided by the scan porosity ϕ CT taken from the segmented image. The Darcy velocity is constant throughout each experiment: it is imposed by the flow rate of the pumps.
The flow fields developed preferential flow paths at early times in experiments 2, 3, and 4, which then widened later. In experiment 1 the pore space dissolved more uniformly with less evident preferential flow pathways, except possibly at the latest times.
The probability density functions of velocity, Fig. 6, show a very wide range of flow speed, over around 8 orders of magnitude for Estaillades and over 10 orders of magnitude for Portland. For Estaillades there is no significant change in the distribution, normalized by the average velocity, over time, indicative of a relatively uniform dissolution. In contrast, in Portland, we initially observe a significant fraction of the pore space in stagnant regions with flow velocities that are several orders of magnitude lower than the average [ Fig. 6c,d]. As dissolution proceeds, the flow becomes concentrated in fast flowing regions and some stagnant regions disappear. In experiment 4, the extra peaks at later times are due to competition between multiple preferential flow paths as seen in Fig. 5d.
Our hypothesis, discussed further below, is that in experiments 2, 3 and 4, the fast flow channels, present initially, dissolve more quickly than the stagnant regions. This causes channels to form through the pore space located in regions where the flow was initially fast. This is distinct from wormholing, since the reactant moves in advance of the channel formation. We test this quantitatively by studying the correlation of velocity in time. This is defined as: with a normalized function: where v i (t) is the velocity in voxel i at time t. The calculation is only performed on voxels that are void in the reference (initial) image at time t = 0. C(0) is the variance of the velocity initially. If dissolution proceeds randomly throughout the pore space, we would expect C(t) to decline over time, as there is no correlation between the velocity in a particular voxel at time t and the same voxel originally. If the dissolution is uniform, expanding existing flow channels, then C(t) may retain an approximately constant value, since we keep the same flow paths. However, if we form channels, then we see the largest flow speeds in those regions of the initial image that were flowing fastest. Since the channels focus the flow, C(t) may increase over time, as the velocities located where there was void originally sample the high end of the velocity distribution. Fig. 7 shows the normalized correlation, Eq. (2). In all the experiments C(t) increases with time, although there is a more marked initial rise for experiments 2, 3 and 4. We see an extraordinary one to four order-of-magnitude growth of the correlation function which provides compelling quantitative evidence that the channels which form, see Fig. 5, are located in regions where the flow was initially fastest and where, as a consequence, dissolution was most rapid.

Dimensionless numbers and reaction rates
We use the definition of Péclet number Pe, Damköhler number Da, and effective reaction rate r eff from Menke et al., 2015. Pe [−] is the ratio of diffusive to advective times at the pore scale: where U av [m·s −1 ] is the average pore velocity, D is the molecular diffusion coefficient of Ca 2+ at 25°C (7.5 × 10 −10 m 2 ·s −1 ) and L [m] is the characteristic length where L = π/S and S [m −1 ] is the specific surface area.
Using the flow rates established by the pumps and the surface area from the scans, the characteristic length L is between 115 μm and 226 μm for the four samples, while Pe has a value between 145 and 3040 at the start of the experiments [ Table 2]. At high Péclet numbers, pore-scale advection is, on average, much faster than diffusion and fresh reactant is available throughout the core.
The Damköhler number [−] is used to characterize reactive conditions and is defined as the ratio of advective to reactive times: where r is the reaction rate constant of pure calcite (8.1 × 10 −4 mol·m −2 s −1 ) measured experimentally by Peng et al., 2015 at the same experimental conditions as here. n is calculated using n = ρ calcite [1 − ϕ total ]/M calcite , ρ calcite is the density of pure calcite (2.71 × 10 3 kg·m −3 ), and M calcite is the molecular mass of calcite (0.1 kg·mol − 1 ). ϕ total , the total porosity, is measured using helium pycnometry for Estaillades (ϕ total = 0.295) (Andrew et al., 2014) and measured using nitrogen porosimetry for Portland (ϕ total = 0.213).
Here we account for both the porosity evident at the resolution of the images, ϕ CT , and micro-porosity. As the experiments proceed, we assume that the intra-granular micro-porosity is unaffected by reaction (this is consistent with no detectable change in raw grey-scale in the grains). Hence we write ϕ total = ϕ grain (1 − ϕ CT ) + ϕ CT to find ϕ grain at the beginning of the experiment: we assume that its value remains constant during dissolution, and recalculate ϕ total as ϕ CT increases over time.   1) and (2), between the initial voxel velocities (time = 0) and those from each successive image through time for each experiment.
Both Da and PeDa, the ratio of the reaction rate to the diffusion rate, are much less than one throughout the experiments, representing conditions at which the reaction rate is slow compared to both the advection and diffusion rates.
The brine pH is estimated to be 3.1 using a geochemical equilibrium calculation where the Gibbs free energy is minimized (Leal et al., 2014) and the chemical conditions are assumed to be relatively constant throughout the sample, as the high Pe means acid flows fast,  continuously providing reactant to the entire pore space. This allows the change in overall porosity to be used in estimating the average effective reaction rate [mol·m −2 ·s −1 ]: where Δ t is the time between scans and Δ ϕ CT is the image-measured change in porosity within that time (we assume, as mentioned above, that the micro-porosity is not affected by the reaction). S [m −1 ] is the specific surface area at the beginning of that time.
Our computed values of Pe, Da, PeDa, and r eff are shown in Table 2 and displayed graphically in Fig. 8. In experiment 1, Pe remains approximately constant with a slow decrease consistent with an increase in porosity leading to a decrease in average velocity, Eq. (3). In other experiments, the initial drop in Pe is more dramatic, but is again followed by a slow, almost linear, decrease with time, again consistent with the evolution of porosity [Fig. 8a]. The magnitude of Pe in experiments 1 & 3 is lower than in 2 & 4 because of the lower imposed flow rate.
Da increases for all the experiments [ Fig. 8b], consistent with the decrease in average velocity, see Eq. (5).
Both PeDa and r eff remain steady for experiment 1, but showed a sharp decrease for experiments 2, 3, and 4 followed by steady rates [Figs. 8c & 8d]. At late time, the reaction rate of experiment 1 remained higher than experiment 3 despite having the same flow rate, while experiments 2 and 4 showed similar effective reaction rates. The most significant feature is that the reaction rates are much lower than the value Table 3 Network and direct flow simulation properties for times 0, 63, and 120 min. of 8.1 × 10 −4 mol·m −2 s −1 measured on a flat surface with no transport limitations: initially the effective rates are at least ten times lower, declining to 100 times lower at late times. This demonstrates that reaction is severely limited by the rate of transport of reactant to the solid surface: when channels form, these dominate the flow and limit reaction to the sides of the channels themselves; in comparison the transport of reactant and products in the pore space away from these channels is very slow. This discrepancy between intrinsic and effective rates is larger than the 14-fold variation observed in the more homogeneous Ketton sample (Menke et al., 2015). Local Péclet statistics generated from the velocity distributions of the direct simulations are shown in Table 3 and displayed graphically in Fig.  9. Each void voxels velocity is normalized to the experimental Darcy velocity and then used with equation 3 to calculate a local Péclet number. The mean local Péclet number (Fig. 9a) decreases almost linearly through time for experiment 1 and has a sharp decrease at earlier times for experiments 2-4 which then become nearly linear at later times. The maximum local Péclet number (Fig. 9b) follows a similar trend with the highest values at the beginning of the experiments and then decreasing in accordance with the constant experimental flow rate and the increase in porosity. The standard deviation of the local Péclet number decreases sharply for experiments 2, 3, and 4, but remains relatively constant for experiment 1 and is consistent with the PDFs of velocity in Fig. 6. Fig. 9d depicts the fraction of void voxels that have a local Péclet number of ≤1 and have diffusion-dominated flow. Portland has a larger portion of diffusion-dominated voxels than Estailladesthis is consistent with the lower permeability found in Portland as seen in Fig. 4b. In experiment 1 the fraction of voxels with diffusion-dominated flow is small and the fraction steadily increases with time. However, in experiments 2, 3, and 4 there is a decrease in the number of diffusiondominated voxels in the first 30 min of the experiment followed by a steady increase at later times. We believe this demonstrates that during the formation of the channels there is competition between preferential flow paths and as a consequence more of the pore space is advection dominated. After the channels are formed the majority of flow is directed through these mobile paths which also leaves larger immobile regions in which the number of diffusion-dominated voxels increases.
Network extraction [ Table 3] indicates that the pore size and throat diameter all increase with time, while the total number of pores and throats decrease. The coordination number (representing the connectivity) is highest initially for the Estaillades samples and increases slightly through time for experiment 1 and dramatically at first for experiments 2, 3, and 4 and then stabilizes at later times.
Traditionally dissolution has been studied experimentally and numerically in relatively homogeneous systems, for which there is a Pe and Da-dependent transition from uniform to wormholing patterns (Daccord et al., 1993a, Daccord et al., 1993bFredd and Fogler, 1998;Ott et al., 2012). In our experiments, we are in what would be considered to be the uniform regime, where the reaction rate is relatively low in comparison to the flow rate, allowing reactant to penetrate the whole sample. Indeed, under the same conditions, uniform dissolution was observed in a relatively homogeneous Ketton limestone [Menke et al., 2015]. However, the dissolution patterns we have observed are not always uniform, but clearly show the emergence of preferential flow paths. This is not the same as wormholing, where the acid etches an advancing channel through the rock, but channelling, where the acid rapidly flows throughout the sample and then widens the fast flow pathways, as illustrated schematically in Fig. 10 and confirmed quantitatively through the increase in the velocity correlation function over time, Fig. 7, and through the examination of local Péclet numbers in Fig. 9.
In wormholing, we see a sudden and very significant increase in permeability when the hole spans the core: this is not observed in our experiments, with a smoother increase in permeability with time, as shown in Fig. 4. Instead, in experiments 2, 3, and 4, the changes in porosity and permeability along with the evolution of the velocity fields (see, for instance, Fig. 5) are indicative of channel formation [ Fig. 10b] rather than wormholing [ Fig. 10a]. In experiment 1, the steady increase in porosity and permeability coupled with the steady reaction rate is similar to uniform dissolution. The lower flow rate combined with the initial high porosity, connectivity, and permeability of Estaillades means that advection dominates in the majority of the pore space and fresh reactant is continuously provided to most areas of the pore space preventing the formation of obvious preferential flow channels [ Fig. 10c].

Conclusions
We have found a clear impact of initial structure on the changes in porosity and permeability due to fluid/solid reaction in carbonate rocks. We see the emergence of a new type of dissolution pattern in heterogeneous systems: channelling, where preferential paths widen and dominate the flow and the dissolution process. The effective reaction rates in the channelling regime are up to two orders of magnitude smaller than that measured on flat substrates with no transport limitations, implying that transport of reactant and products in regions of the pore space away from the channels is significantly impeded.
Fast synchrotron tomography coupled with a reservoir-condition flow apparatus is a powerful experimental tool that can be adapted to explore a range of applications including multiphase flow processes (Ott and Oedai, 2015), advection-dispersion, and transport in chemically heterogeneous media. Such an approach allows for integration of experiment and modelling at reservoir conditions in complex systems, elucidating dynamic reaction processes on a pore-by-pore basis. Fig. 10. A schematic of wormholing (a), uniform (b), and channel formation (c) dissolution regimes. In the wormholing regime, the reactive brine (blue) floods the pore space (black) dissolving the rock (white) as it advances. Black dots denote intragranular porosity. In the uniform dissolution regime, brine floods the whole pore space rapidly and the rock then dissolves across its entire surface area. During channel formation, flow is faster than in the wormholing regime: reactive brine floods the rock very quickly, but the heterogeneity of the rock allows preferential flow paths to form which then provide the majority of the reactant to a more limited portion of the surface of the rock. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)