Pore-space structure and average dissolution rates: A simulation study

We study the inﬂuence of the pore-space geometry on sample-averaged dissolution rates in millimeter-scale carbonate samples undergoing reaction-controlled mineral dissolution upon the injection of a CO 2 -saturated brine. The representation of the pore space is obtained directly from micro-CT images with a resolution of a few microns. Simulations are performed with a particle tracking approach on images of three porous rocks of increasing pore-space complexity: a bead pack, a Ketton oolite, and an Estaillades limestone. Reactive transport is simulated with a hybrid approach that combines a Lagrangian method for transport and reaction with the Eulerian ﬂow ﬁeld obtained by solving the incompressible Navier-Stokes equations directly on the voxels of three-dimensional images. Particle advection is performed with a semianalytical streamline method and diffusion is simulated via a random walk. Mineral dissolution is deﬁned in terms of the particle ﬂux through the pore-solid interface, which can be related analytically to the batch (intrinsic) reaction rate. The impact of the ﬂow heterogeneity on reactive transport is illustrated in a series of simulations performed at different ﬂow rates. The average dissolution rates depend on both the heterogeneity of the sample and on the ﬂow rate. The most heterogeneous rock may exhibit a decrease of up to two orders of magnitude in the sample-averaged reaction rates in comparison with the batch rate. Furthermore, we provide new insights for the dissolution regime that would be traditionally characterized as uniform. In most cases, at the pore-scale, dissolution preferentially enlarges fast-ﬂow channels which greatly restricts the effective surface available for reaction.


Introduction
The flow of reactive fluids through geologic reservoirs is a complex phenomenon that occurs over a wide range of spatial and temporal scales. There is a strong coupling between flow and reaction, not directly described in terms of an ideal reactor or equilibrium models [Steefel et al., 2005]. Fluid-mineral interactions occur inherently at the pore (micron) scale, and the representation of pore-scale heterogeneities in effective (averaged) models is not straightforward. This makes the study of reactive phenomena at the pore-scale particularly relevant. It provides information, not readily available from field data, that needs to be accounted for when modeling the behavior of reactive fluids percolating through sedimentary rocks. Understanding the relation between the degree of heterogeneity of a porous rock and the average dissolution rate is fundamental to advance the development of upscaling techniques for reactive flows and to model reactive phenomena at the field scale.
Carbonate reservoirs hold the majority of the world's oil reserves [Schlumberger, 2007], and calcite is the prevalent carbonate mineral in deep-sea sediments [Morse and Arvidson, 2002]. Carbonates are very reactive minerals, and alterations in the formation fluids (either natural or human-induced) can lead to dramatic variations in the petrophysical properties of the reservoir, e.g., Schechter and Gidley [1969]; Noiriel et al. [2009]; Luquot and Gouze [2009]; Lamy-Chappuis et al. [2014]. In this work, we explore the impact of the pore-space complexity of real porous media on the average dissolution rates, and on the permeability evolution via a first principles particle tracking approach that simulates carbonate dissolution at reservoir conditions (pressure is fixed at 10 MPa and temperature at 50 o C, representative of an aquifer at 1 km depth).
Laboratory measured reaction rates are known to be, typically, several orders of magnitude higher than the ones estimated from field studies in natural systems. The observed reaction rates at the reservoir scale are affected by different physical and chemical heterogeneities that manifest themselves at different spatial scales, this explain the discrepancy between intrinsic rates (e.g., those measured in ideal reactors) and field rates. Among other factors, this discrepancy has been attributed to incomplete mixing of solutes [Li et al., 2008], variations in fluid residence times [Maher, 2010], and the presence of physical and chemical heterogeneities [Salehikhoo et al., 2013;Steefel et al., 2005;Ellis et al., 2011]. The impact of the velocity field heterogeneity on dissolution rates was studied by Molins et al. [2012] in two-dimensional systems. They showed that the average rate decreases when increasing the complexity of the pore space. Studying homogeneous reactions in different pore structures Alhashmi et al. [2015Alhashmi et al. [ , 2016 showed that the effective reaction rate depends on the degree of mixing, flow field heterogeneity, and transport conditions. In this work, we quantify the reduction of average rates at the millimeter scale for monomineralic rocks caused by pore-scale physical heterogeneities.
Time-resolved pore-scale imaging of dissolution [Menke et al., 2015] showed that, in the uniform dissolution regime, the average dissolution rate of a carbonate rock is one order of magnitude lower than the intrinsic rate. Further study [Menke et al., 2016] demonstrated that for rocks with a more heterogeneous initial pore structure, dissolution in the same regime proceeded through the growth of high-flow channels in the rock, significantly restricting the volume over which reaction occurred, with average reaction rates around 100 times smaller than the batch values. This was termed a channelled dissolution regime. In this paper, we will confirm this conceptual picture and we will use pore-scale simulations, which in previous work were used to predict the results of pore-scale imaging experiments Menke et al., 2015], to obtain a quantitative estimate on the average dissolution rates of heterogeneous rocks (a bead pack, an oolite and a bioclastic limestone). We also provide estimates for the increase in the representative elementary volume (REV) after dissolution and show that changes in permeability are sample-dependent. This paper is organized as follows. We start by describing the flow, transport, and reaction model proposed by Pereira Nunes et al.
[2016]-section 2. Next, we discuss the characterization of transport and reaction in terms of the nondimensional P eclet and Damk€ ohler numbers (section 3). A description of the samples and of the simulation strategy is found in sections 4 and 5. Section 6 presents the simulation results for the long-time behavior upon dissolution at high-flow rates. It is shown that the average dissolution rate depends on the rock heterogeneity and on the flow rate. Section 7 investigates the early time behavior and the permeability evolution at low-flow rates. Finally, section 8 summarizes our main results.

Fluid Flow, Transport, and Reaction
Here we summarize the assumptions and implementation of the particle-based dissolution method introduced by Pereira Nunes et al. [2016].
Solute transport is modeled with a particle (Lagrangian) approach. This method obtains transport characteristics by simulating individual stochastic particle trajectories, instead of solving directly the transport equations [Kinzelbach, 1988;Tompson and Gelhar, 1990;Salamon et al., 2006]. The solute concentration is recovered by calculating the average number of particles in a reference volume at a given time. We adopt a hybrid scheme: particles are transported in an Eulerian flow field resulting from the numerical solution of the incompressible Navier-Stokes equations. Transport is divided into two steps: advection and diffusion. During advection the particles move along streamlines, which do not traverse solid voxels. Diffusion is modeled as a random particle motion, allowing particles to jump across streamlines and to cross the pore-solid interfaces, thus creating a flux of particles from the pores to the solids. This defines the reaction rate, as discussed in section 2.3.

Velocity Field and Permeability
The velocity field is obtained by numerically solving the incompressible Navier-Stokes equation in the voxels of segmented micro-CT images with a finite-volume method implemented with an open source library [OpenFoam, 2013;Raeini et al., 2012;Bijeljic et al., 2013a]: where l is the fluid viscosity, P is the pressure, q the fluid density, andṼ the velocity. The boundary conditions are constant pressure at the inlet and outlet faces and no-slipṼ 50 at the interfaces.
The image-calculated permeability (k) is determined from Darcy's law: where L s is the sample length, A is the cross-sectional area, Q is the total flow rate at the sample outlet, and DP is the pressure drop.

Particle Transport
The displacement of a particle p after a time-step (Dt) has an advective and a diffusive contribution: where DX p 5X p ðt1DtÞ2X p ðtÞ andX p ðtÞ is the position of a particle p at a time t.

Streamline Tracing and Tortuosity
During advection, the particles follow streamlines traced with the semianalytical method developed by Pereira Nunes et al. [2015]. This method is based on quadratic interpolations of the velocity field within each grid block, consistent with the numerical solution of equation (1) at grid faces and with the pore-scale boundary conditions:Ṽ 50 at the solid interfaces. This semianalytical method eliminates numerical stability issues inherent in conventional time-stepping approaches and allows the time-step to be constrained by diffusion only.
Flow-weighted tracing of streamlines is used to estimate the hydraulic tortuosity (T) [Duda et al., 2011]: where hL sl i is the average length of the streamlines and L x is the sample length. The tortuosity provides a sample-scale estimative of the impact of pore-space geometry and flow field heterogeneity on the particle transport. High values correspond to more complex geometries.

Diffusion
Diffusive displacements are modeled with a random walk method. At each time-step, the diffusive jump of the particle is given by [Benson and Meerschaert, 2008]: andñ is a vector of uniformly distributed random components. This generates a three-dimensional stochastic motion with mean-free path ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 6D m Dt p , where D m is the molecular diffusion coefficient. The simulation time-step is calculated by limiting the maximum diffusive jump. Imposed to be l=2, where l is the voxel length. This leads to:

Mineral Dissolution in the Particle Approach
Mineral dissolution is determined by the flux of particles through the pore-solid walls. The model assumes first-order kinetics, meaning that the reactive flux at the mineral surface is proportional to the particle concentration, and reaction-limited dissolution, implying that no significant concentration gradients develop in the vicinity of the surface.
Once a stationary state is reached, the diffusive flux equals the reactive flux at the surface. The reactive flux is determined by the intrinsic reaction rate r ½mol:m 22 :s 21 , i.e., the rate free of mass transport limitations. Intrinsic rates are usually measured in batch reactors (stirred reactors closed to mass transfer) [Pokrovsky et al., 2005;Peng et al., 2015] and in this work we use both terms interchangeably. The diffusive flux in a fluid with reactant concentration C ½n: of particles:m 23 is determined by the average number of particles crossing the interface in a time step. By equating the number of mineral moles dissolved to the local flux of particles, it is possible to derive the number of particle hits (N hits ) required to dissolve one voxel comprised of N voxel mol moles of the dissolving mineral: This establishes a direct relationship between the intrinsic (batch) dissolution rate and the particle concentration without any adjustable parameters since the number of hits to dissolve is calculated in terms of a priori information, namely the intrinsic reaction rate (obtained from independent measurements) [Peng et al., 2015], the diffusion coefficient, the time step, and the mineral density.
In practical implementations of the method, the number of particles is fixed, and when they react (collide with the solids), they are reinjected at the inlet. This creates an imbalance in the concentration profile. The reaction rate is equalized by adjusting dynamically the number of hits to dissolve as a function of the sliceaveraged concentration. This is necessary to keep the ratio N hits =C (and hence the reaction rate) constant. This will be discussed further in section 5.
In a first-order surface reaction, dissolution is proportional to the local flux of particles at the pore-solid walls, which in turn is determined by diffusion. The flux is calculated by counting the number of particles that cross from the surrounding pore voxels to a specific solid voxel. This model represents conditions under which the fluid is constantly out-of-equilibrium with respect to the rock. Once a particle crosses the solid-fluid interface it reacts (it is removed from the simulation) and the number of hits which that particular solid voxel has received is updated. After a certain number of hits (equation (9)), the solid voxel is fully dissolved into a pore voxel.

Calcite Dissolution
The main product of calcite (CaCO 3 ) dissolution in the presence of CO 2 -saturated brines is Ca 21 , which is released in the aqueous phase. This dissolution kinetics has been extensively discussed elsewhere, e.g.,  [1989]. Throughout this work, we consider r58:1310 24 mol:m 22 :s 21 , the rate corresponding to calcite dissolution in a CO 2 -saturated brine at 10 MPa and 50 o C [Peng et al., 2015].

Water Resources Research
Pereira Nunes et al.
[2016] introduced a model that describes the dissolution kinetics in terms of the concentration of a single kind of particle, which is, in turn, interpreted as the Ca 21 undersaturation, or the deviation from the equilibrium concentration, in the liquid phase. The fluid is in equilibrium with respect to the solid when C 5 0; in this case, the flux of particles to the solids is zero and no dissolution occurs.
3. Description of Transport and Reaction: The P eclet and Damk€ ohler Numbers P eclet and Damk€ ohler numbers are used to characterize reactive transport. The P eclet number (Pe) compares the rate of mass transfer due to advection to the rate of mass transfer by diffusion over a characteristic length (L): here q5Q=A the Darcy velocity. For a porous rock, we adopt L5pV=S A , where V is the volume of the sample and S A is the area of the pore-solid interface [Mostaghimi et al., 2012]. This is a good approximation for the effective diameter of the rock grains extended to cases where identification of a typical grain size is difficult, or impossible, such as in many carbonates. When Pe > 1, transport is advection-dominated.
The (advective) Damk€ ohler number (Da) is the ratio of the advective time-scale and the reaction time-scale [Lasaga, 1998]. It is the ratio of the time to dissolve a volume V of a porous rock with molar density n ½mol:m 23 and surface area S A over the advective time t adv 5L/=q ð Þ : At the interfaces,Ṽ 50 and transport is diffusion-dominated. The product PeDa is equivalent to the diffusive Damk€ ohler number-the ratio between the reaction rate and the diffusion rate-and describes the dissolution kinetics. When PeDa < 1, the kinetics is reaction-limited, and when PeDa > 1 it is diffusion-limited.
As discussed in the introduction, behavior diagrams have been proposed to classify dissolution patterns observed in Darcy-scale simulations and experiments in terms of the Pe and Da numbers. While qualitatively correct and useful to interpret core-scale results, their direct application to pore-scale scenarios is not straightforward. This will be further explained in sections 6 and 7.

Sample Characterization
We selected three images of porous media with distinct geometries, petrophysical properties, and provenience to study the dynamics of carbonate dissolution: an artificial random packing of identical spheres, an Estaillades carbonate (a bioclastic limestone), and a Ketton oolite. The bead pack image is based on measurements of the bead centers in a monodisperse bead pack by Finney [1970], for which subsequent segmentation was performed by Prodanovic and Bryant [2006]. The Ketton and Estaillades images were acquired using a Versa XRM-500 X-Ray Microscope (Zeiss X-Ray Microscopy, Pleasanton, CA, USA) and segmented with a watershed algorithm [Andrew et al., 2014;Menke et al., 2015]. Both Estaillades and Ketton are mostly monomineralic rocks, composed of more than 95% calcite [Andrew et al., 2014]. Figure 1 shows the pore space geometry and the corresponding computed velocity fields for the three samples. The grains in the Ketton carbonate are roughly spherical and the pores are well connected. There is a significant flow throughout almost the entire pore space in both Ketton and the bead pack. In comparison, Estaillades has a less connected pore space and the flow field is highly channelized. This is reflected in its tortuosity, the highest of the three samples studied. Table 1 summarizes image-calculated sample properties: porosity, permeability, tortuosity, and characteristic pore length. Figure 2 shows the probability density function (PDF) of the horizontal velocity component. The bead pack has the narrowest distribution followed by Ketton and Estaillades, which has a much wider distribution.

10.1002/2016WR019313
Wide velocity distributions are manifestations of the presence of significant stagnant zones, regions of slow flow, with restricted solute access in advection-controlled flows [Bijeljic et al., 2013b]. Ultimately, this flow characteristic is responsible for the divergence between average rates and batch rates. Furthermore, the qualitative similarity between their corresponding PDFs explains why Ketton and bead pack exhibit a similar behavior upon dissolution, as discussed in section 6. From the geological perspective, the Ketton rock is an oolite, a limestone comprised of roughly spherical grains (ooids) of calcite albeit with significant cementation. This suggests that it retained some of its original features as a packing of spheres.
Subresolution porosity is accounted for by means of a factor that corrects the number of mineral moles in the solid voxels according to standard helium porosimetry measurements. The grain porosity (/ grain ), the porosity of the solid phase, is estimated by comparing the helium porosity (/ HP ) and the micro-CT porosity (/): Based on Andrew et al. [2014], we calculate / grain 50:12 for Ketton and / grain 50:22 for Estaillades. The grain porosity is then used to calculate the effective number of calcite moles in the solid voxels:  (5)) is higher in the Estaillades sample, which has the most heterogeneous flow field.

Simulation Strategy
The flow simulations are performed imposing an initially prescribed Pe, which is attained by establishing a pressure gradient that satisfies the desired Darcy velocity. 2:03 10 4 particles are injected with a flow weighted rule at the inlet and tracked until they either collide with a solid through diffusion or exit the domain. In both cases, a new particle is injected at the inlet, thus keeping the number of particles constant. The reaction rate is fixed (section 2.4), and the molecular diffusion coefficient is that of Ca 21 : D m 57:5310 210 m 2 :s 21 ]. This results in: PeDa5 1:8310 22 for the bead pack, PeDa57:0310 22 for Ketton and PeDa54:0310 22 for Estaillades. Dissolution is in the reaction-controlled regime: fast transport and slow reaction.
Following dissolution, the geometry of the pore space changes. Therefore, the velocity field must be updated. This is done at regular intervals when the bulk porosity increases by a fixed amount-D/ 0:001. This is consistent with the scenario in which localized increases in porosity do not lead to a large variation in the global pressure field. The choice of D/ is a compromise between computational performance and accuracy [Pereira Nunes, 2016]. Using smaller values of D/ does not affect the results. The reaction rate is equalized throughout the sample by introducing a correction in the number of hits to dissolve. We replace the volumetric particle concentration (C) in equation (9)  where s x stands for the tomographic slice number in the flow direction. This ensures that the reaction rate (r) is constant and eliminates spurious boundary effects caused by the higher concentration of particles at the inlet.
The main steps of the algorithm are: 1. Solve the incompressible Navier-Stokes equation. 2. Flow weighted particle injection at the inlet. 3. Particle advection along semianalytical streamlines. 4. Particle diffusion by a random walk. When particle collides with a solid, kill it and reinject at the inlet. If the number of hits on any solid voxel is equal or larger than N eff hits (equation (14)), dissolve the voxel solid (change from solid to void) and update the porosity. 5. Repeat steps 3 and 4 for every particle and time step. 6. Update the flow field (solve the Navier-Stokes equation) when D/ 0:001 and update the number of hits to dissolve-equation (14). 7. Repeat items 3-6 until the end of the simulation.

Results: Average Dissolution Rates at High P eclet Numbers
This section studies the asymptotic sample-averaged dissolution rates at high-flow rates (high P eclet number). We simulate the long-time behavior of Ketton, Estaillades, and the bead pack under two initial flow rates, corresponding to Pe 5 10 and Pe 5 100. Table 2 summarizes the simulation parameters and the expected dissolution regime according to the macroscale classification from Golfier et al. [2002]; dissolution is in the uniform regime for Pe 5 100 and in the wormholing regime for Pe 5 10.
At the Darcy scale, the uniform dissolution regime is associated with slow reactions and high-flow conditions that favor the spreading of fresh reactants throughout the sample. Under these conditions, Menke et al. [2015] and Pereira Nunes et al. [2016] discussed how sample-scale transport limitations imposed by a heterogeneous flow field are responsible for a decrease of one order of magnitude in the average reaction rate when compared to the batch rate. It was shown that this flow and reaction regime, termed ''uniform'' in Darcy-scale studies, is actually a regime where the average dissolution rate is severely constrained by porescale heterogeneities.
When the flow rate is high, the average particle displacement is large enough to propagate the reaction throughout the sample, from inlet to outlet. We can then compute the average dissolution rate between times t and t 0 from the micro-CT surface area (S A ðtÞ) and the change in porosity (D/ 0 ): The average rate as a function of the number of pore volumes injected (PVI) is investigated for the three samples while varying the flow rate.
The porosity profiles (the slice-averaged porosities along the direction of flow) before and after reaction are shown in Figure 3. The three samples experience an increase in porosity across the flow axis, evidencing that the flow rate is high enough to carry the solute throughout the samples. Since the residence time is higher for Pe 5 10 than for Pe 5 100, the increase in the bulk porosity is also higher, for the same amount of fluid injected. As a consequence of the increase in  porosity, the permeability also increases. In Estaillades, it starts from 0:5 D and reaches 8:4 D and 5:1 D after 1000 PVI, for Pe 5 10 and Pe 5 100, respectively. In Ketton, it goes from 13:3 D to 128:0 D and 70:9 D, for Pe 5 10 and Pe 5 100, respectively. In the bead pack, permeability increases from 13:9 D to 61:9 D and 31:8 D, for Pe 5 10 and Pe 5 100, respectively. Figure 4 shows the ratio r eff =r for the three samples under two initial flow rates, corresponding to Pe 5 10 and Pe 5 100. The Ketton and the bead pack samples (blue and green lines, respectively) have comparable average rates, a consequence of their broadly similar velocity fields (Figure 2), while the Estaillades sample has a lower rate. We contend that this is a manifestation of its heterogeneous velocity field, which causes a very inefficient sampling of the pore space. Large portions of the Estaillades pore space remain virtually unaffected by reactive transport, even after injecting 1000 PVI, Figure 5. In contrast, practically all of the pore space of the Ketton and bead pack is modified after injection of the reactive particles. The channelized nature of the Estaillades velocity field confines the flux of solute to very localized regions, thus limiting the spread of dissolution to preferential channels, as observed experimentally by Menke et al. [2016]. The other two, more homogeneous rocks, exhibit a much more uniform dissolution pattern, Figure 5.
Dissolution increases the connectivity and enlarges the pore throats, it is thus expected that the tortuosity (equation (5)) should decrease as a consequence. In fact, this is observed in the three samples. Following the injection 1000 PVI with initial Pe 5 100, we observe that in the bead pack it decreases from T 5 1.25 to 1.16, from T 5 1.54 to 1.31 in Ketton, and from T 5 1.75 to 1.32 in Estaillades.
Transport is advection-dominated and the presence of stagnant zones inhibits the access of the reactants to a significant part of the reactive surface. Moreover, a high-flow rate (Figure 4, solid lines, Pe 5 100) results in higher average reaction rates in comparison with the lower flow rate (Figure 4, dashed lines, Pe 5 10). The increase from Pe 5 10 to Pe 5 100 results in an increase of almost one order of magnitude in the reaction rate for the three samples. Zones that have low velocity when Pe 5 10 become accessible to the particles as the velocity rises, resulting in a more even spread of the dissolution front, Figure 6.  (15)) over the batch reaction rate (r) during injection of a CO 2 -saturated brine. The average rates for the more heterogeneous rock (Estaillades, in red) are lower than the rate of the relatively homogeneous ones (Ketton in blue and bead pack in green). Higher flow rates yield higher average reaction rates -Pe 5 100, solid lines. The dashed lines are for Pe 5 10. Figure 5. Dissolved portions (shown in red) superimposed with the original pore space (light grey) after injecting 1000 pore volumes, starting with an initial P eclet number of 100. Dissolution in the Estaillades sample (a) is more channelized and restricted than in the Ketton (b) and bead pack (c) samples. This is a result of the flow field heterogeneity (Figure 1) and results in a lower average dissolution rate, as seen in Figure 4.

10.1002/2016WR019313
Despite being classified as in the wormholing regime according to the Pe and Da numbers, characteristic wormhole-like channels, where macroscale paths are etched through the rock, are not present in the dissolution pattern for Pe 5 10, for any of the samples studied, Figure 6. Instead, reaction proceeds through the widening of preexisting flow channels, with some reaction throughout the sample. These channels are not wormholes, as they are a pore-scale phenomenon occurring in a slow reaction regime (Da 1) and do not retain all the reactant: dissolution proceeds along the length of the sample through the growth of fast-flow pathways, not from the etching of new sample-scale holes. Qualitatively similar behavior is seen for both Pe 5 100 and Pe 5 10.
The understanding of the reaction dynamics can be further strengthened by analysing the evolution of the velocity field distribution. As shown in Figure 5, while the long-term dissolution of Ketton and bead pack alters most of the original pore space, significant parts of Estaillades remain largely unchanged. Figure 7 shows the corresponding PDFs of velocities before and after the injection of 1000 PVI for the three samples for Pe 5 100. Estaillades evolves into a bimodal distribution, where the stagnant zones are mostly preserved by the reactions, while the high-flow regions become more pronounced. This is a markedly different behavior from the one exhibited by the Ketton and bead pack samples, in which the broad characteristics of the distributions remain. The bead pack PDF is virtually unchanged, and the Ketton PDF becomes narrower with a more pronounced peak. This narrowing of the distribution, in the case of Ketton, indicates that the pore space is becoming less complex and more well connected as a result of dissolution: more pores have velocities close to the average. The opposite is happening in Estaillades; dissolution preferentially enlarges preexisting fast flow channels, while keeping the stagnant zones unchanged. The resulting flow field is more heterogeneous and its distribution is bimodal.
Overall, these results are in qualitative agreement with the experimental results of Menke et al. [2016], who showed that dissolutioninduced changes in the velocity field depend on the initial pore structure-in more heterogeneous rocks, dissolution occurs preferentially in channels (regions of fast flow) leaving stagnant zones largely unaffected.
The representative elementary volume (REV) is defined as the minimum sample volume above which statistical variations of a given property become negligible [Bear, 1988]. Its  value is property dependent, i.e., the permeability REV is not necessarily the same as the porosity REV, the former being normally larger due to long-range correlations in the flow field. Macroscale, imagecalculated properties are meaningful only when obtained from volumes larger than the REV. As a consequence of dissolution, both the average pore-size and the connectivity increase. It is thus expected that the REV also increases with the emergence of new structures, particularly if dissolution channels are formed. A useful way of estimating the REV is provided by the empirical (semi)variogram, which describes the behavior of data correlation with distance. For a given property p, it is defined as: where h is the distance between two observations and N(h) is the number of observations such that jpðx i1h Þ2pðx i Þj5h [Doyen, 2007]. It measures the dissimilarity between data separated by a distance h. Once the variogram reaches a plateau, the data are uncorrelated. Figure 8 shows the porosity and velocity variograms before and after reaction. The variograms are plotted as a function of the horizontal distance (distance measured from inlet) divided by the initial characteristic length. All the three samples experience an increase in the correlation length for both porosity and horizontal velocity: larger REVs are necessary to characterize these rocks after dissolution. This effect is stronger in Estaillades, where we see a large increase in the velocity correlation length. This is explained by the intense channelization of the pore space, which preserves the correlation of velocity. The effect of the flow and reaction in Estaillades is to enlarge high-flow channels and decrease tortuosity, therefore, increasing the similarity between velocities at different points inside the channel. A simple analogy is to consider the flow field inside a tube: if the velocity does not change longitudinally, as each point on a streamline is perfectly correlated with the previous one.
Even though the flow and reaction regime is the same for all three samples considered-Pe 1 and PeDa 1-we observe that the average rate is both rate and sample-dependent. For Estaillades, the most heterogeneous sample, the average rates are consistently lower than the ones for the more homogeneous Ketton and bead pack, Figure 4.
The heterogeneous nature of the flow field causes an uneven distribution of the reactants throughout the pores, thus restricting dissolution to regions with significant intake of solute. Accordingly, dissolution is more uniform in the two samples with more uniform velocity fields. However, even in these cases, the average reaction rate is still one order of magnitude lower than the batch rate.

Results: Initial Pore-Structure and Dissolution at Low P eclet Numbers
While the previous section focused on the asymptotic behavior of the average dissolution rates at high P eclet number and low Damk€ ohler number conditions, this section studies the early time behavior and the influence of the initial pore geometry on the permeability evolution of Ketton and Estaillades at low-flow rates (low P eclet number), see Table 3 for the simulation parameters.
The penetration length of a reactive particle inside a sample is dependent on the P eclet number. Our results show that, for Pe 5 1 and Pe 5 2, and injecting 10PVI, dissolution is limited to regions close to the inlet. This means that, for a fixed number of injected pore volumes, if the injection rate is low enough, dissolution will not propagate throughout the sample. Figure 9 shows the slice-averaged relative porosity difference in Ketton and in Estaillades along the direction of flow for four initial injection rates, corresponding to Pe 5 1, 2, 5, and 10. In both samples, the difference is concentrated close to the inlet for the lowest rates (Pe 5 1, 2), and moves toward the outlet as the rate increases, i.e., the transport of reactants in the pore space is more effective and dissolution spreads throughout the rock. Note, however, that the penetration length is higher in Estaillades than in Ketton. This is due to the channelized characteristic of the Estaillades flow field, which favors the development of a localized dissolution around a preexistent channel, an effect that is more pronounced for higher Pe, Figure 10. In both samples, even at low Pe, distinct channels are etched through the rock.
The sample-dependent difference in the range of the dissolution (Figures 9 and 10) should have an obvious impact on the permeability behavior. This is illustrated in Figure 11, which compares the evolution of the permeability as a function of porosity for both rocks. The alterations in Ketton are concentrated at the inlet, and do not connect across the sample, even for Pe 5 10. Consequently, despite the increase in porosity, the permeability remains essentially unaffected for the flow rates considered. The opposite happens in Estaillades, where we observe the widening of a channel extending from inlet to outlet. This occurs at the lower flow rate corresponding to Pe 5 5. This localized enlargement of a highly conductive channel connecting inlet and outlet is responsible for a marked Figure 9. Slice-averaged relative porosity difference versus distance over the characteristic length (L), after injecting 10 pore volumes in (top) Ketton and in (bottom) Estaillades. Dissolution penetrates further for the higher flow rates, particularly in Estaillades, the more channelized rock. Flow is from left to right. Ketton is 8.2 L in length and Estaillades 7.9 L. increase in permeability. Furthermore, the nonuniqueness of the permeability versus porosity curve evidenced by Figure 11 has potential impacts for the modeling of dissolution in reservoir-scale models. These are currently constructed based on the assumption that the entire permeability field is a well-defined function of the porosity [Steefel et al., 2014]. However, Figure 11 shows that such assumption does not hold. Instead, the permeability behavior upon dissolution is a function not only of porosity, but also of the pore-space morphology and the flow rate. For instance, the distribution of tracer transit times, and consequently of the flow rate and P eclet number, is known to be very wide in heterogeneous reservoirs [Di Donato et al., 2003;Di Donato and Blunt, 2004]. This suggests that an adequate description of the evolution of permeability upon dissolution cannot be obtained in terms of a single permeability versus porosity relation.

Conclusions
This paper presented a modeling study to illustrate the influence of the initial pore structure on the transport of reactants in porous rocks and the associated impact on the sample-averaged fluid-solid reaction rates. A more complex pore space gives rise to a more heterogeneous flow field, which imposes limitations on the transport of reactive solute at the millimeter scale. These pore-scale constraints manifest themselves collectively as a sample-scale discrepancy of, at least, one order of magnitude between average reaction rates and batch reaction rates, even in relatively homogeneous rocks. Another indication of the transport-controlled kinetics at the sample scale is the dependency of the average rates on the flow rates-higher flow rates yield higher average reaction rates.
As illustrated in section 6, the underlying complexity of the pore geometry (and of its corresponding flow field) defines the average reaction rate in asymptotic times, even in the so-called uniform dissolution regime-Pe 1; PeDa 1. Despite being useful in several circumstances, Figure 10. Dissolution fronts following the injection of 10 pore volumes for four flow rates, corresponding to Pe51; 2; 5, and 10. For low-flow rates (Pe51; 2), reaction is contained close to the inlet. Increasing the injection rate (Pe55; 10), the channelized nature of flow in the Estaillades focuses the reaction into high-flow channels, with a greater impact on permeability. Figure 11. Permeability evolution for Ketton (dashed lines) and Estaillades (solid lines) following the injection of 10 pore volumes with four different flow rates, corresponding to Pe 5 1, 2, 5, and 10. The permeability of Ketton is largely unaffected for the range of flow rates considered: a consequence of dissolution being concentrated at the inlet. In contrast, Estaillades shows a marked increase in permeability for the higher rates: a result of fast flow channels dissolving.

Water Resources Research
10.1002/2016WR019313 this commonplace definition (''uniform'') does not describe what is happening at the pore scale. Specifically, in heterogeneous rocks, we observe the preferential dissolution of fast-flow channels, and the flow and dissolution regime is more accurately described as channelled, rather than uniform. This was also reported in an imaging experiment tracking the pore-space evolution during reactive transport [Menke et al., 2016].
The distinct evolution of the velocity PDFs of rocks with different pore characteristics is a strong indication that a description of dissolution patterns solely based on the nondimensional numbers Pe and PeDa misses important features of the reaction dynamics in real rocks. Previous work has focused on reproducing core experiments and adopted a Darcy-scale continuum representation of the porous medium. This approach fails to account for pore-scale heterogeneities and does not fully describe the dynamics of dissolution in terms of the pore structure [Daccord et al., 1993;Golfier et al., 2002;Maheshwari et al., 2013]. This paper is a contribution in this direction: through computational modeling we investigated the interdependency of the dissolution process and the geometry of the pore space, which, in turn, impacts the average reaction rate and the dependence between permeability and porosity. Furthermore, while our focus was on the dissolution of calcitic carbonates, reservoirs of commercial interest are usually formed by a combination of calcite and dolomite (and other minerals in smaller concentrations). The simulation method employed here can be readily expanded to treat chemically heterogeneous rocks by using a multi-phase segmentation (e.g., pores 1 calcite 1 dolomite) and attributing different reaction rates to each kind of solid voxel.
Fluid-solid reactions are broadly characterized as either surface or transport-controlled. This a useful description in batch reactors, closed systems in which the reactants are readily available in the solution in contact with the mineral surface. We argue that such description is incomplete when considering reactive transport in porous media. There is one extra factor to be taken into account: the geometry of the pore space, which may limit the transport of solute to the solid surface at the sample scale. This means that a reaction can be surface-controlled at the pore scale (a rough analogy would be to consider every pore as a batch reactor) but transport-controlled at the sample scale. This was the case exemplified here.