The Impact of Pore Structure Heterogeneity, Transport, and Reaction Conditions on Fluid–Fluid Reaction Rate Studied on Images of Pore Space

We perform direct numerical simulation using a pore-scale fluid–fluid reactive transport model (Alhashmi et al. in J Contam Hydrol 179:171–181, 2015. doi:10.1016/j.jconhyd.2015.06.004) to investigate the impact of pore structure heterogeneity on the effective reaction rate in different porous media. We simulate flow, transport, and reaction in three pore-scale images: a beadpack, Bentheimer sandstone, and Doddington sandstone for a range of transport and reaction conditions. We compute the reaction rate, velocity distributions, and dispersion coefficient comparing them with the results for non-reactive transport. The rate of reaction is a subtle combination of the amount of mixing and spreading that cannot be predicted from the non-reactive dispersion coefficient. We demonstrate how the flow field heterogeneity affects the effective reaction rate. Dependent on the intrinsic flow field heterogeneity characteristic and the flow rate, reaction may: (a) occur throughout the zones where both resident and injected particles exist (for low Péclet number and a homogeneous flow field), (b) preferentially occur at the trailing edge of the plume (for high Péclet number and a homogeneous flow field), or (c) be disfavored in slow-moving zones (for high Péclet number and a heterogeneous flow field).


Introduction
The study of reactive transport in the subsurface is important to understand contaminant transport (Appelo and Postma 1993), nuclear waste disposal (Bodvarsson et al. 1999), and CO 2 sequestration (Lichtner et al. 1996). Fluid-fluid reaction is encountered in a number of processes, such as the creation of acidic water as a result of a mineral spill in the soil, landfill biodegradation, and the migration of radioactive materials (Appelo and Postma 1993). One of the main complexities in understanding these phenomena stems from the fact that geological porous media are heterogeneous at all scales. The nature of reactive transport in the subsurface is affected by both the degree and the scale dependence of heterogeneity, starting from the pore scale. Hence, investigating the impact of the pore structure of the medium is key to understand the coupled transport and reaction behavior. In this study, we will identify the factors governing fluid-fluid reactive transport using a recently developed pore-scale model that performs direct numerical simulation on images of pore space.
The advection-dispersion-reaction equation (ADRE) has been proven to have weaknesses in predicting the behavior of reactive solute in porous media (Kapoor et al. 1997;Raje and Kapoor 2000;Gramling et al. 2002). It overestimates the effective reaction rate as it assumes a complete mixing between reactants. Despite that, a significant body of previous work has focused on studying the validity of the ADRE to predict fluid-fluid reactive transport behavior. These studies either applied averaged equations, which assume a complete mixing of reactants, or used fitting parameters calibrated to experimental data. Sanchez-Vila et al. (2010) developed a continuum-scale time-dependent model by which reaction rate coefficients are adjusted with time based on the ADRE. Their model successfully predicted fluid-fluid reactive transport data in a beadpack (Gramling et al. 2002). Their model involves three fitting parameters: the first is the longitudinal dispersion coefficient, which is crucial to estimate the amount of pore-scale mixing and thus the amount of reaction, while the two others define the time-dependent reaction rate. Chiogna and Bellin (2013) presented an analytical model that assumes a distribution of mixing ratio within a reference elementary volume (REV). Their model employs a calibrating parameter based on the peak concentration of the product. Hochstetler and Kitanidis (2013) upscaled a 2D porous medium comprised of randomly distributed equal-diameter circular grains into a 1D reactive transport numerical model. The model is based on an empirical relationship that was used as a calibrating parameter to find the effective reaction rate based on the ADRE. The validity of ADRE model on the basis of theoretical upscaling has been discussed elsewhere (Porta et al. 2012). Using different methods, at the Darcy scale, Edery et al. (2009Edery et al. ( , 2010 accurately reproduced the experimental data of Gramling et al. (2002) based on a continuous-time random walk (CTRW) model. The model introduced a non-local formulation of the dispersive transport term, so they were able to account for the time-dependent behavior of the dispersion coefficient. Their model requires two fitting parameters. The first is used to determine the reaction radius under which reaction may occur, and the second defines the shape of the probability density function (PDF) of concentration.
Several authors have used pore-scale simulation to avoid the use of averaged equations such as the ADRE. Willingham et al. (2008) conducted micromodel experiments of bimolecular reactive transport in different porous media structures and compared it with a lattice Boltzmann model. The porous medium was a 2D system and had different pore shapes and sizes. They demonstrated the impact of pore structure on reaction and showed that the interfacial contact area between fluids is important in determining the amount of pore-scale mixing and thus the amount of reaction. De Anna et al. (2014) presented a 2D pore-scale experiment of fluid-fluid reactive transport in a porous medium consisting of circular grains that were randomly distributed. They compared the results against a model based on imperfect mixing to estimate the concentration profile of the product of an irreversible bimolecular reaction. They demonstrated that the reaction rate is controlled by the geometry of the interface where reactants mix. Advective stretching and molecular diffusion control the geometry of this mixing zone. Their model calibrates the stretching rate and the initial concentration to fit the data. These models are based on simple 1D or 2D systems, which do not capture the complete physical and chemical behavior of reactant species in real 3D porous media. Notably, in this issue Edery et al. (2016) have studied the impact of conductivity heterogeneity on the dynamics of 2D reactive transport at the core scale and showed the importance of preferential pathways. Nevertheless, there is still lack of understanding of the impact of heterogeneity at the pore scale on effective reaction rate under different transport and reaction conditions in 3D porous media.
We use a model for simulating fluid-fluid reactive transport directly on images of the pore space (Alhashmi et al. 2015), which has been validated against benchmark experimental data conducted in a beadpack (Gramling et al. 2002). The model is based on a particle-tracking method that uses a streamline method to capture advection (Pereira Nunes et al. 2015) and a random walk to simulate diffusion with a probability of reaction (P r ) to represent the chemical kinetics. Dynamic changes in reaction rate in a beadpack were predicted accurately and shown to be affected not only by incomplete mixing as it had been demonstrated before experimentally (Raje and Kapoor 2000;Gramling et al. 2002;Jose and Cirpka 2004;Rolle et al. 2009;Edery et al. 2015) and numerically (Edery et al. 2010;Bolster et al. 2010;Dentz et al. 2011;Ding et al. 2013;Hochstetler and Kitanidis 2013), but also as a result of early-time pre-asymptotic spreading (Alhashmi et al. 2015). The characterization of the pre-asymptotic behavior through a continuum-scale model has been recently discussed by Porta et al. (2016), based on previous work by Porta et al. (2012).
The approach of Alhashmi et al. (2015) is based on using pore-scale images that capture the irregular shape of the pore structure: In this work, we study images that consist of up to 10 9 voxels, allowing an accurate description of the complex pore-scale geometry of reservoir rocks from the µm to mm scale. In addition, the method allows the study of fluid-fluid effective reaction rates in porous media based on the known independently measured batch reaction rate in a free fluid that is related to the reaction probability. We investigate the effect of pore structure heterogeneity on the dynamic effective reaction rate by performing simulations on three images for a range of transport conditions (or Péclet number, Pe) and a range of reaction conditions (or reaction rate constant, k r ). We interpret the results through analysis of the velocity field and reactive transport heterogeneity characteristics of the three samples and by comparing the dispersion characteristics of the porous media studied for the reactive and non-reactive cases.

Model Description
The model consists of four main parts: geometry, flow, transport, and chemical reaction. A more detailed description can be found in Alhashmi et al. (2015); here, we describe the salient features pertaining to this work. Table 1 shows the properties of the images used in the model and Tables 2, 3 conditions applied in the model simulations.

Pore Geometry
We use three pore-scale images: a beadpack, Bentheimer sandstone, and Doddington sandstone. The beadpack consists of a random close packing of equally sized spherical grains whose coordinates have been measured (Finney 1970). Prodanović and Bryant (2006) performed the segmentation of this beadpack into a pore-scale image. The dry-scan images of Bentheimer sandstone (provided by iRock Technologies) and Doddington sandstone (provided in-house) were acquired by Xradia Versa micro-CT scanners on cylindrical cores of 5 mm diameter and 25 mm length. The image segmentation into void and solid voxels was performed using a seeded watershed algorithm. All image processing was done using the Avizo Fire 7.0 program (Visualization Sciences Group, Burlington, MA, USA; www.vsg3d. com). All pore-scale images are provided as Cartesian grid blocks (voxels) with the beadpack consisting of 500 × 500 × 500 voxels, having a voxel size of 3 µm and a porosity (φ) of 36.23 %. Bentheimer has 1000 × 1000 × 1000 voxels, φ = 21.47 %, and a voxel size of 3.0035 µm, while Doddington has 1000 × 1000 × 1000 voxels, a voxel size of 2.6929 µm, and φ = 19.49 %. The characteristic length (l in µm) of Bentheimer and Doddington can be estimated by Mostaghimi et al. (2012): where V is the volume of the porous medium (µm 3 ) and S the total surface area between pore and solid voxels measured from micro-CT images (µm 2 ). For a packing of spheres, the beadpack, l is taken to be equal to the diameter of the grains. The value of l for the images is given in Table 1.

Flow Simulations
To obtain the flow field in the pore space, we solve incompressible Newtonian flow governed by the Navier-Stokes equations: where u is the velocity vector (µm/s), P is the pressure, ρ is the density (here we use water with ρ = 1000kg/m 3 ), and μ is the dynamic viscosity (μ = 1 × 10 −3 Pa s). Each image voxel represents a grid block in the flow simulation. The single-phase flow field on the image voxels is calculated based on an implementation of the finite volume method (Raeini et al. 2012;Bijeljic et al. 2013b) using OpenFOAM (OpenFOAM 2011). The pressure implicit with splitting of operators (PISO) algorithm is implemented to solve the velocity field. The normal and tangential components of the velocity are zero on the solid boundaries. A constantpressure boundary condition at the inlet and the outlet faces of the images is applied, while on the other faces no-slip boundary conditions are implemented. The pressure drop in the model is selected so that creeping flow with a Reynolds number 1 is studied. Figure 1 shows the pore space and computed velocity fields.
The velocity distributions of the three images are significantly different, as presented in Fig. 2. This figure shows PDF of the ratio of the absolute values of the center voxel velocity in the direction of flow (u x ) to the average pore velocity (u avg , which is equal to the ratio of Darcy velocity, q, to φ), against normalized velocity (u x /u avg ) presented in semi-log plot (Fig. 2a) and log-log plots ( Fig. 2b) (Bijeljic et al. 2013b). The velocity distributions are sampled uniformly in 256 bins of log(u x /u avg ). It is evident from Fig. 2 that the peak of velocity distributions, which is at the average velocity (u x /u avg = 1), is largest in the beadpack. Close to this velocity value, there are nearly twice as many voxels in the beadpack as compared to the two sandstone rock images implying that the beadpack has a more homogeneous flow characteristic (more velocities are similar to the average velocity) than Bentheimer and Doddington. Moreover, the distributions for Doddington and Bentheimer are spread more widely in comparison with the beadpack: they have a larger number of fast velocities and, in particular, a significantly higher number of stagnant velocities. This means that their Pore space (left) and normalized velocity field (right) for a the beadpack, b Bentheimer sandstone, and c Doddington sandstone. The normalized velocity field presented as the ratio of voxel velocity to the maximum voxel velocity, u max transport characteristics are expected to be more non-Fickian (Bijeljic et al. 2013a, b). In addition, the visual presentation of the velocity field (see Fig. 1) in the regions faster than average velocity implies that the beadpack pores are both better connected and less tortuous than the other two images representing consolidated sandstone. The remarkable feature is that for the sandstones the velocity distribution spans eight orders of magnitude. Therefore, Fig. 2 Probability density function (PDF) of velocity distributions for the beadpack (blue circles), Bentheimer (red crosses), and Doddington (green squares) presented as a semi-log and b log-log plots even for very fast average flows, there will still be significant portions of the pore space where diffusion limits transport.

Transport and Reaction Simulation
We study a second-order irreversible bimolecular reaction where A + B → C. Reaction is simulated such that initially we have equal amounts of A and B in the system. Reactant A (N tot A particles) is placed uniformly at random throughout the pore space, while all reactant B (N tot B ) is placed in the first layer of voxels using a flux-weighted rule, to represent the injection of a pulse of B. As the image size in the direction of flow of the beadpack is half the size of Bentheimer and Doddington, we placed reactant A over 2 images in the beadpack, whereas only a single image is used for the sandstones. We have the same value of the concentration of A, c A0 , in all cases: The numbers of particles used are given in Table 1. In each time step, we track particles movement by advection and diffusion: where the vector x labels the position of a particle. For advection (x Advection ), the model uses a semi-analytical streamline-based approach to move particles through the grid cell. It uses a novel method to account for zero flow at the solid wall boundaries, while using closed form expressions for the flow field within a grid block, consistent with the computed fluxes at block faces (Pereira Nunes et al. 2015). For diffusion (x Diffusion ), the particle motion is based on a random walk. It is a series of random spatial displacements based on the diffusive step (ξ in µm) which defines particle transitions (Bijeljic et al. 2013b) such that: where t is the time step size (s), D m is the molecular diffusion coefficient (µm 2 /s), and x Diffusion has components x, y, and z. ϕ and θ are random numbers in the range from 0 to 2π and 0 to π, respectively. If a particle crosses the exit face of the image domain in the diffusive part of time step, we randomly inject it at the inlet face of the next image domain according to the cross-sectional area-weighting rule. On the other hand, if it crosses during the advective part of the time step, we use a flux-weighted rule. The re-injection is performed randomly in order to avoid the possibility of introducing correlation in the system. To reduce the computational cost, we set t to the largest possible time such that ξ is equal to the voxel size. In this study, the smallest voxel size of the image studied (Doddington) was taken to find ξ (Table 1). Hence, t is equal to 1.7 × 10 −2 s. We calculate the longitudinal dispersion coefficient (D L ) that describes the spreading of particles as defined by Bijeljic et al. (2004): where σ 2 is the variance of the particles displacement. Figure 3 shows a comparison of the ratio of asymptotic D L to D m as a function of Péclet number (Pe) for transport simulations performed on our beadpack, Bentheimer, and Doddington images. In the same figure, we also provide previously published transport simulation data for a sandpack, Berea sandstone, and Portland carbonate together with the compilation of experimental data for unconsolidated beadpacks/sandpacks, as analyzed in Bijeljic et al. (2011). The results for Bentheimer and Doddington agree well with the previously published results for Berea-the sandstone data lie between the beadpack/sandpack data and the carbonate data, further confirming that transport heterogeneity increases from the beadpack to sandstones to carbonates as a result of an increase in pore structure and velocity field heterogeneity (Bijeljic et al. 2013a, b). Furthermore, Fig. 3 illustrates the difference in D L (amount of spreading) for the beadpack Fig. 3 Ratio of asymptotic D L to D m as a function of Pe for the beadpack (blue), Bentheimer (red), and Doddington (green). The black solid line represents predictions for previous transport simulations on a sandpack, while the dashed line is for Berea sandstone, and the dotted line is for Portland carbonate (Bijeljic et al. 2011). The circles are measured data for unconsolidated beadpack/sandpack (Pfannkuch 1963;Bijeljic and Blunt 2006) and different rock images, which exists as a result of different nature of transport through their pore structures. This difference is larger at higher Pe numbers. Later, in Sect. 3, we will use transport simulation results to interpret the variation in dynamic reaction rates for different pore structures as a result of different degrees of spreading. We allow reaction in every time step if it meets two conditions. First, reactants A and B must be less than a distance ξ (Eq. 5) apart. Second, the reaction occurs with a probability P r that is related to the reaction rate constant independently measured in a free fluid (k r in µm 3 /moles s) such that: where n is the number of moles that each particle in the simulation represents. The closest reactant pair is considered first, when there is more than one possible reaction within a distance ξ . After reaction, the product C replaces both reactants and is placed in the position equidistant between A and B. P r must lie between zero and 1; so, if all the parameters in Eq. (8) are kept constant, the k r will have an upper limit. In our simulations, we consider a range of rates k r that can be modeled with our approach. This method has been validated for a free fluid system as well as for experimental data in porous media (Alhashmi et al. 2015).
The advantages of our model are that: (1) it uses 3D micro-CT images that account for the irregular shapes of pore structure which is very important to capture the flow and transport properties and thus accurately predict the effective reaction rate; and (2) it can predict the effective reaction rate behavior in porous media based on known transport conditions and an independently measured batch reaction rate. Equation (8) relates the reaction rate that can be measured independently to reaction probability in the model.
The aim of our study is to look into the impact of pore structure on the effective reaction rate for a range of transport and reactive conditions. We define dimensionless numbers determining these conditions as: where t Dif is diffusive timescale and t Adv is the advective timescale. Da is Damköhler number known as the ratio between the reaction rate and the diffusion rate. It can be defined as the ratio of t Dif to the reactive timescale, t R . In each image, we perform simulations using our reactive transport model for different ranges of Pe and P r (or k r ), as shown in Table 1.

Results and Discussion
We study the effect of pore structure and the flow field heterogeneity on the dynamic effective reaction rate in the beadpack, Bentheimer, and Doddington images. Our analysis and interpretation consists of visualizing reactive plumes, presenting the reaction rates and reaction-altered dispersion coefficient and comparing it with the results for non-reactive transport, and providing velocity distributions of the voxels in which particles were injected, placed, and in which the product was formed. Figure 4 shows the plumes for the three cases presented as 3D coordinates of particle (reactants A and B and product C) positions in the pore space captured at t = 102 s for Pe = 1, 10, 100, 500, and 2240 and P r = 1. In this figure, we define particle dimensionless positions in 3D with the coordinates x D = x x avg ·t , y D = y y max , z D = z z max , where x, y, and z are the coordinates of particle displacement (µm) from the inlet. The direction of flow is in the x dimension, where y max and z max are the pore-scale image size in the other dimensions (Table 1). x D = 1 indicates that particles are located according to the mean displacement for every transport condition. In this figure, reactant A is represented with blue points, reactant B with red points, and product C with green points. As expected, the product C is located in the mixing zone between the two reactants-note that at Pe = 10 the injected reactant B appears to overlie most of the product C and reactant A, which is due to the large number of particles remaining unreacted at this time. The greater progression of the product C plume in the pore space with an increase in Pe in each pore-scale image illustrates that an increase in the amount of pore-scale mixing leads to more reaction; this can also be seen as smaller , and product C (green) at t = 100 s for Pe = 1, 10, 100, 500, and 2240 and P r = 1 reactant plumes. In addition, it demonstrates the degree of flow heterogeneity in the porescale images (described by velocity distribution in Fig. 2) as clearly seen for Pe = 100 and upward. Since the beadpack has the least heterogeneous flow characteristics, most of the particles travel close to the average pore velocity and almost no particles are trapped in the low velocity regions. On the other hand, Bentheimer as the second most heterogeneous system in terms of the velocity distribution has some particles retarded in the stagnant regions, while Doddington as the most heterogeneous has even more particles located in these immobile zones. This is especially the case for initially resident particles A that can remain unreacted and trail the injected and product particles at high Pe, as clearly observed in Fig. 2.

Effective Reaction Rates and Dispersion Coefficients
We define a product ratio, f C , as the moles of product C divided by the initial moles of reactant A. Figure 5 shows the overall rate of reaction r = d f C dt as a function of time t for Pe = 500 and P r = 0.01 to illustrate both the early-and late-time behavior. Some oscillations are seen in the reaction rate in the beadpack: As reactants travel along the converging streamlines of the narrow channels, the mixing causes higher reaction rates, followed by the lower rates as the particle trajectories diverge. Figure 6 shows a range of transport (Pe = 1, 10, 100, 500, and 2240) and reaction conditions (P r = 0.01, 0.1, and 1). Figure 7 shows the reaction rate plotted against 1/ √ t. Beyond a simulation time t Dif , between 260 and 569 s dependent on the rock studied (see Table 1), there is little change in the behavior with similar-and low-reaction rates in all cases. Once the particles have had sufficient time to diffuse a characteristic length, they are well mixed and the reaction rate is largely independent of the local pore structure.
The effective reaction rate is mainly controlled by a combination of transport and reactive factors. The transport factors are the amount of mixing between reactants due to diffusion, the amount of mixing as a result of the spreading (dispersion) of reactants, and the degree of flow/transport heterogeneity. For the same rock, higher P r means higher k r and therefore more reaction product is generated.
At late time, as shown in Fig. 7, we see r ∼ 1/ √ t for low Pe, which is expected. If we assume asymptotic dispersion and reaction-at late time-that can be modeled by the ADRE, the injected pulse of B spreads a distance √ D L t through A and it is in this volume that reaction can occur: The rate is proportional to ∂ However, for higher Pe and early time, this scaling is not followed so closely, indicating that a simplistic formulation based on a dispersive mixing with a fixed coefficient cannot capture the behavior: in these cases, there is significant dispersive spreading, but this does not mean that the reactants are Fig. 5 Reaction rate r as a function of t for Pe = 500 and P r = 0.01. Blue solid line is for beadpack, red is Bentheimer, and green is for Doddington. Shown in the inset are the results at early time Fig. 6 Reaction rate r as a function of t for different Pe and P r . Blue solid line is for beadpack, red is Bentheimer, and green is for Doddington mixed at the pore scale; instead, they may remain separated, flowing in different fast channels, or stuck in separated stagnant zones. Figure 6 shows that, generally as Pe increases, the amount of product C created increases. This is because the amount of mixing between reactants as a consequence of spreading grows. To interpret the reaction rate behavior as a function of time, we present longitudinal spreading (dispersion) for the simulations in which reaction takes place (reactive case) and compare it with the transport simulations only (non-reactive case). Figure 8 shows the longitudinal dispersion coefficient for initially resident particles A versus time Fig. 7 Reaction rate r as a function of 1/ √ t for different Pe and P r . Blue solid line is for beadpack, red is Bentheimer, and green is for Doddington for P r = 0.01, 0.1, and 1 for Pe = 1, 10, 100, 500, and 2240. The dispersion coefficients differ in the reactive and non-reactive cases for high Pe, which is indicative of a selective removal of particles by reaction, dependent on how they have moved relative to the average displacement.
At Pe = 1, the amount of fluid-fluid reaction at early time in the beadpack is larger than in the two rock images, as shown in Fig. 6. The amount of mixing due to diffusion dominates at low Pe (Sahimi 1995;Bijeljic et al. 2004) and at early times: here the more open and homogeneous pore space of the beadpack leads to a slightly larger effective dispersion coefficient and greater production of C compared to Bentheimer and Doddington. At late time, the degree of mixing must be similar, since the reaction rates are much the same, even though Bentheimer has a slightly higher dispersion coefficient. The dispersion for reactive and non-reactive conditions is similar, which implies that reactant is produced with the same likelihood throughout the region through which the particles of A and B are dispersed.
For Pe = 10, the impact of mixing due to spreading is more significant than that due to diffusion. This is shown in Fig. 8 with much higher values for the longitudinal dispersion coefficient for the two sandstones compared to the beadpack in both reactive and non-reactive cases. Again the dispersion coefficients are similar for reactive and non-reactive conditions, so reaction occurs throughout the dispersed region. Transport in the beadpack reaches the asymptotic regime sooner than for Pe = 1, but with a value of similar magnitude. At early time, shown in Fig. 6, the creation of C is largest in the beadpack, as there is more rapid diffusion-mediated mixing. Later, the reaction rate in Bentheimer sandstone is largest as a result of the highest D L (Fig. 8) followed by Doddington sandstone and then by the beadpack.
For Pe = 100, we see that there is less dispersion with reaction for Bentheimer and more for Doddington (Fig. 8). For Bentheimer, this implies that reaction is favored at the edge of the plume, removing these particles and restricting the spread of unreacted A, compared to a case with no reaction. In contrast, for Doddington, more reaction occurs for particles moving closer to the average speed with reaction less favored at the fringes. Bentheimer also has the highest D L at early times: there is more spreading and this allows reaction to occur over a larger area and the reaction rate is largest. Later, D L becomes larger in Doddington for reactive transport (Fig. 8), leading to higher reaction rates compared to the other two samples at intermediate times (Fig. 6). Here reaction does not occur preferentially at the edges of the plume: there is more pore-scale mixing allowing reaction to occur throughout the sample.
At Pe = 500 and Pe = 2240, the intermediate-and late-time production rate of C in the beadpack is greatest: this is most pronounced for low P r (k r ), as shown in Fig. 6. The dispersion coefficient does not control the reaction rate, since D L is much higher in the sandstones, Fig. 8. Moreover, the dispersion is higher for reactive transport, which indicates that particles at the edge of the plume are less likely to react. This is because the fastest moving particles of B are confined to fast flow channels and do not mix with A in the slowermoving regions of the pore space. This is evident in Fig. 4 for Pe = 500 and Pe = 2240, where slow-moving particles of A are left unreacted in the sandstones, as opposed to the beadpack, where reaction proceeds preferentially at the trailing edge of the plume, leading to less spread of A and a lower reactive dispersion coefficient. These effects will be quantified through the analysis of PDFs of the voxels in which particles reacted in the next section. In the advection-dominated regime, we have more spreading compared to mixing at the pore scale. The beadpack has less distinction between fast and slow channels and a more open pore space, which allows a greater degree of mixing locally (see Fig. 9 in the next section). This competes with the much larger extent of spreading in the sandstones, which allows the reactant to explore a bigger volume of the pore space. However, if advection dominates over diffusion, injected particles may be retained within fast channels and do not diffuse out to encounter A in the stagnant zones. At late time, the reaction rate is largest in Bentheimer and lowest in Doddington, reflecting a different pore-structure-dependent balance of mixing and spreading. Probability density function (PDF) of velocity for Pe = 1 and Pe = 500 at P r = 1 for a the beadpack, b Bentheimer, c Doddington. Blue is for the injected particles B at t = 0, black is for the placed (resident condition) particles A at t = 0, magenta is for the voxels in which product particles C were formed from injection until t = 1 s, green is for the product particles until t = 5 s, and red is for the product particles until t = 100 s Figure 9 shows the comparison of the PDFs of velocities for the voxels where particles B are injected (blue curve) and particles A are placed (resident particle condition, black curve), and the voxels where product C is formed (magenta curve from injection until t = 1 s, green curve is until t = 5 s, and red is until t = 100 s) for P r = 1, Pe = 1 and Pe = 500, for the three rocks. Since the particles B are injected according to a flux-weighted rule, their Probability density function (PDF) of velocity of voxels in which product particles were formed at Pe = 1 and Pe = 500 which are presented from injection until t = 100 s for a the beadpack, b Bentheimer, c Doddington. Black is for P r = 0.01 and blue is for P r = 1 PDF of voxel velocities shows few values below the average velocity, which is expected. Particles A, on the other hand, are placed throughout the image volume in the pore space, and thus their velocity distributions take the overall shape of the intrinsic PDFs of velocity in the entire image pore space, as shown in Fig. 2. For both Pe = 1 and Pe = 500, as the reactant particles sample more of the flow heterogeneity (with increasing time), the PDF of velocity of the formed product gradually approaches the shape of the PDF of velocity distributions of resident particles A that represents the flow heterogeneity characteristic of the porous medium. This is observed for all three media studied.

Fig. 11
Probability density function (PDF) of velocity distributions of the voxels in which product particles were formed for P r = 1 and from injection until t = 100 s for a the beadpack, b Bentheimer, c Doddington. Blue is for Pe = 1, red is for Pe = 10, green is for Pe = 100, yellow is for Pe = 500, and black is for Pe = 2240 The PDF of voxel velocities in which the product is formed is shifted to higher values for the sandstones at Pe = 500-the reaction occurs preferentially in higher-velocity regions, leaving unreacted A in slow-moving zones, as evident in Fig. 4 for Pe = 500 and revealed through the higher dispersion of A compared to a case with no reaction, as shown in Fig. 8 for the same Pe. For the beadpack, in contrast, we see that the voxel velocities in which the product was formed were preferentially sampled in the regions with a slightly lower-than-average velocity: this serves to limit the spread of A through consuming the trailing edge of the plume, leading to a lower dispersion coefficient in the reactive case.
At lower Pe, more mixing occurs between fast and the slow regions-consequently, sampling of the slow regions becomes more frequent and therefore the PDFs of the two sandstones from injection until t = 100 s at Pe = 1 are closer to the PDF of the resident particles than at Pe = 500. The full sampling of voxels of pore space is achieved most rapidly in the beadpack since it has more open pore space leading to better mixing locally as clearly shown in Fig. 9-this process is more gradual in the more heterogeneous sandstone images.
Next we illustrate the coupled effect that the reaction probability has on transport, which can lead to different effective reaction rates. Figure 10 shows the comparison of the PDF of velocity of the voxels in which product particles were formed at Pe = 1 and Pe = 500 for P r = 0.01 and 1, presented for until t = 100 s for the three images studied. For low Pe, the PDF of velocity for the two values of P r is similar. This confirms that reaction is mostly controlled by the amount of mixing due to diffusion. However, at high Pe, the PDF of velocity for different P r values differs for the more heterogeneous sandstones. In Bentheimer and Doddington sandstone images, the lower P r has the lower peak around the average velocity, complemented by the larger number of less mobile voxels in which reaction has taken place. This suggests that, under the same transport conditions, the reactant particles which for higher P r do react in the fast channels for the lower P r = 0.01 do not-in the latter case the reactant particles survived for longer and through mass transfer exchange had more time to sample the slow regions in which they subsequently react. This behavior is most profound in Doddington since it is the image with the most heterogeneous flow characteristic.
Finally, we observe the impact of transport conditions on reactions for the three images studied. Figure 11 shows the comparison of the PDF of velocity of voxels in which product particles were formed at P r = 1 for the range of Pe studied and from injection until t = 100 s. Higher Pe results in higher voxel velocities in which the product was formed, which is shown in Fig. 11 as the shifts in the PDFs of velocity to larger values. Particles B are predominantly injected in the mobile domain (see Fig. 9) and, since at higher Pe there is less time for mass transfer exchange to the immobile regions, the reaction tends to occur in higher-velocity voxels. This effect is more profound in the sandstone images as they have more fast voxels in their intrinsic PDFs of velocity of the entire pore space (Fig. 2). On the contrary, at lower Pe there is more time for mass transfer exchange so mixing due to diffusion causes better sampling of the low velocity regions-this is shown in Fig. 11 as more of the slower voxel velocities in which reaction has taken place.

Conclusions
We have implemented a fluid-fluid pore-scale reactive transport model developed by Alhashmi et al. (2015) to study the impact of pore structure and flow heterogeneity on dynamic effective reaction rates under a range of transport and reactive conditions. We studied three pore-scale images: a beadpack, Bentheimer sandstone, and Doddington sandstone.
The effective reaction rates estimated by the model vary with different pore structures that define the flow fields and therefore the transport characteristics of the system. We quantified the behavior by studying the interplay between reaction rate, the PDF of velocity, and dispersion coefficient, comparing it to non-reactive conditions. The higher the Pe and/or P r is, the larger the effective reaction rate. The reaction rate is a subtle balance of three factors whose relative contribution is different dependent on Pe, P r and pore structure: (1) the amount of mixing between reactants due to diffusion, (2) the amount of mixing between reactants due to spreading, and (3) the degree of heterogeneity in the flow field.
We found that at Pe = 1, the reaction rate in the system is predominantly controlled by the amount of mixing due to diffusion: reaction occurs throughout the plume with similar dispersion coefficients for reactive and non-reactive transport, and similar PDFs of velocity for injected, resident, and product particles. At intermediate and high Pe number, the amount of reaction is controlled by the combination of pore-scale mixing due to spreading and the degree of heterogeneity in the flow field. The higher degree of flow field heterogeneity in Bentheimer and Doddington sandstones as compared to the beadpack may allow for more mixing and faster reaction. However, for very fast flows, the injected reactant is confined in fast flow channels and does not encounter particles in more stagnant zones, which inhibits the later-time reaction rate: in such cases, a more homogeneous rock, with less spreading, can allow more reaction, since the reactants are better mixed locally. This is also seen in the timedependent dispersion coefficients, with a higher dispersion coefficient for resident particles in the sandstones compared to the non-reactive case, since reaction is inhibited in the slowmoving regions of the plume, with the resident particles bypassed by the injectant; similarly, the PDF of velocity in which the product was formed shows a shift to higher velocities compared to the resident particles. In the beadpack, the opposite is seen, with preferential reaction at the trailing edge of the plume, less spreading of resident particles, and a shift in the PDF of velocities in which the product was formed toward slightly lower-than-average flow speeds. A decrease in intrinsic reaction rate favors reaction in the lower velocity regions, since injected particles have additional time to sample more of the flow field before reacting. With an increase in Pe, there is less time for mass transfer and we observe more frequent reaction in the higher-velocity mobile regions, where most of the reactant is injected. This effect is more pronounced in media with more heterogeneous flow characteristics.