GeoChemFoam: Direct modelling of multiphase reactive transport in real pore geometries with equilibrium reactions

We present the novel numerical model GeoChemFoam, a multiphase reactive transport solver for simulations on complex pore geometries, including microfluidic devices and micro-CT images. The geochemical model includes bulk and surface equilibrium reactions. Multiphase flow is solved using the Volume-Of-Fluid method and the transport of species is solved using the Continuous Species Transfer method. The reactive transport equations are solved using a sequential Operator Splitting method, with the transport step solved using our OpenFOAM-based Computational Fluid Dynamics toolbox, and the reaction step solved using Phreeqc, the US geological survey's geochemical solver. The model is validated by comparison with analytical solutions in 1D and 2D geometries. We then applied the model to simulate multiphase reactive transport in two test pore geometries: a 3D pore cavity and a 3D micro-CT image of Bentheimer sandstone. In each case, we show the pore-scale simulation results can be used to develop upscaled models that are significantly more accurate than standard macro-scale equilibrium models.


Introduction
Reactive transport in porous media is an essential field of study with broad ranging applications in a range of industries including oil and gas production, carbon dioxide (CO 2 ) and hydrogen (H 2 ) storage, geothermal energy production, nuclear waste disposal and subsurface contaminant transport [1]. These processes include fluid flow with inertia and viscous effects, advective species transport, molecular diffusion, and chemical reactions. In addition, multiple fluid phases are often present, resulting in capillary effects and interface transfer. For such complex systems, analytical solutions are restricted to very simple geometries and flow conditions [2,3]. These limitations in model complexity result in the use of experiments to investigate more complex physics with small analogue systems such as core samples [4,5] or reactive micromodels [6,7]. However, experimental studies are often time-consuming, limited in size, and hard to control. In addition, reactive transport experiments on core samples are always destructive, and since no two cores are the same, they cannot be repeated on identical natural pore structures. These studies are thus often coupled with numerical simulations, a powerful tool that can be used both during the design of the experiment to choose optimal conditions, or to augment the experimental data afterwards by providing quantities of interest that may be difficult to measure (e.g. pH) or to explore additional ranges of physical conditions (e.g. pressure, temperature) [8].
While numerical modelling of multiphase flow [9,10,11] and single-phase reactive transport [12,6,13] in pore-scale geometries have been extensively investigated independently, few studies have attempted to model the coupling between the two. Raoof et al. [14] used a pore network model to simulate reactive transport in variably saturated porous media. However, the pore network approach introduces restrictions on the transport regime and reactive surface area calculations. Chen et al. [15] employed the Lattice Boltzmann Method (LBM) to model multiphase reactive transport, with an interfacial reaction treatment rather than a direct modelling of interfacial conditions. Although this method has been used successfully in several studies [16,17], the LBM method has difficulty modelling the full range of regimes that occur during multiphase flow [11] and reactive transport [18]. Aziz et al. [19] investigated wettability alteration during low-salinity flooding using a non-reactive multiphase transport solver based on Direct Numerical Simulation (DNS). However, the model was restricted to transport in the aqueous phase with an immobile non-aqueous phase and no interfacial conditions. None of these studies include accurate modelling of interfacial conditions with phase transfer.
Recent advances in the development of DNS of multiphase transport have enabled accurate modelling of interfacial transfer. Haroun et al. [20] introduced the single-field approach to model species transport in multiphase systems with interfacial conditions. Their method is based on the Volume-Of-Fluid (VOF) method [21], where the interface between the two fluids is captured using an indicator function, which is a phase volume fraction. Although other methods such as level-set [22,23] can provide a more accurate description of the sharp interface, the VOF method is attractive due to its accuracy of mass conservation and adaptability to more complex physics. Marschall et al. [24] developed Haroun's single-field approach into a versatile and precise method for multiphase transport during bubbly flow labelled Continuous Species Transport (CST). This method was extended to problems with moving contact lines by Graveleau et al. [25] and later improved by Maes and Soulaine [26] with the introduction of interface compression. The CST method was then used to model multiphase reactive transport during low-salinity flooding [27] and mineral dissolution with CO 2 production in shale formations [28,29]. Finally, the model was extended to include local volume changes induced by interface transfer for simulating dissolution of CO 2 bubbles in liquid [30,31,32].
The objective of this paper is to present our extended model that includes multiphase reactive transport with equilibrium reactions both in the water phase and at the surface of the solid. The fully-coupled multiphase reactive transport model is presented in section 2 and validated in section 3. In particular, we show that precise representation of interfacial conditions is essential for accurate and robust modelling of reactive transport, even when the species only exist in one phase, demonstrating that the CST method can be used for reactive transport, unlike the model presented in [19]. Finally, we present the simulation and upscaling of reactive transport with two model test cases: (1) First, we simulate carbonic acid formation during dissolution of a CO 2 gas bubble in a 3D pore cavity and then (2) we introduce the first results of a multiphase reactive transport simulation on a real 3D pore space with injection of a CaCl solution into a micro-CT image of Bentheimer sandstone.

Model description 2.1 Geochemical model
We consider a multiphase system with a reactive phase p in a chemical model that includes N c and N s bulk and surface components, with N x and N y bulk and surface equilibrium reactions. Since the species are at chemical equilibrium, it is possible to partition the system into N c = N c − N x and N s = N s − N y , the primary bulk and surface species (i.e. species with independent concentrations), and N x and N y , the secondary bulk and surface species [33]. The equilibrium chemical reactions between the primary and secondary species can be written as where A j and A i are the chemical formulas of the primary and secondary species in the bulk phase, S m and S n are the chemical formulas of the primary and secondary species on the solid surface, and ν ij and ν nj are the stoichiometric coefficients. Note that on the solid surface, one secondary species is associated to one primary species only. Each equilibrium reaction provides an algebraic link between the primary and secondary species via the law of mass actions where a j,p and a i,p are the activities of primary species j and secondary species j, ω m and ω n are the activity of the primary surface species m and secondary surface species n, and K i and K n the chemical equilibrium constants. We assume that the activity of a species k in phase p is equal to where γ k,p is the activity coefficient of species k (primary or secondary) , c k,p is its concentration (kmol/m 3 ) in phase p and C 0 = 1 kmol/m 3 is the standard activity. The activity ω l of a surface species S l (primary or secondary) is equal to its mole fraction on the corresponding surface, i.e. over all surface components which share the same primary species S m . For each primary bulk species j, we also define the total concentration ψ j,p in phase p, which is the quantity conserved during chemical reactions, and can be written as where Γ is the site density (kmol/m 2 ) and A s is the specific surface area (m 2 /m 3 ) of the solid which, at the pore-scale, is calculated from the mesh. For surface reactions, the apparent stability constant K n is different from the intrinsic constant K i n due to the surface charge q where v n is the charge of the surface species n and F (= 9.649 × 10 7 C/kmol) is the Faraday constant. The double-layer surface potential ϕ is related to the surface charge by the Grahame equation [34] where ǫ (= 78.41 at 25 o C) is the dielectric constant of pure water, ǫ 0 (= 8.854 × 10 −12 C/V/m) is the vacuum permittivity, I (kmol/m 3 ) is the ionic strength of the electrolyte solution, R (= 8.314 kJ/kmol/K) is the ideal gas constant and T is the temperature. The relationship between K n and K i n is given by [34] where Z n is the net change of surface charge induced by the reaction. In this work, activity coefficients, ionic strength, surface charge, surface potential, and chemical equilibrium constants are calculated within Phreeqc [35].

Multiphase flow model:VOF
In this study, the system includes two phases: the aqueous phase (phase 1) and a non-aqueous phase (phase 2), that can be either a gas or a liquid phase. In the VOF method, the interface between the two fluids is tracked using indicator functions α 1 and α 2 , where α 2 = 1 − α 1 , which are equal to the volume fractions of each phase in each grid cell. The density and viscosity of the fluid in each cell are expressed using their single-field values where ρ p (kg/m 3 ) and µ p (Pa.s) are the density and viscosity of phase p. Similarly, the velocity and pressure in the domain are expressed in term of the single-field variables where u p (m/s) and p p (Pa) are the velocity and pressure in phase p. Each phase is assumed to be Newtonian and incompressible, and fluid properties are assumed to be constant in each phase (and in particular independent of the phase composition). In this case, the single-field momentum equation [21] can be written as where g (=9.81 m/s 2 ) is the gravity acceleration and f σ is the surface tension force where σ (N/m) is the interfacial tension, n 12 is the normal vector to the fluid/fluid interface, going from phase 1 to phase 2, κ = ∇ · n 12 is the interface curvature and δ 12 is a Dirac function located at the interface. At the triple point fluid/fluid/solid, the interface forms with the normal to the solid surface a contact angle θ so that n 12 = cos θn s + sin θt s , where n s and t s are the normal and tangent vectors to the solid surface, respectively [36]. In addition, the single-field continuity equation writes whereṁ 12 (kg/m 3 /s) is the rate of mass transfer from phase 1 to phase 2 by solubility, and is calculated after solving the transport equations. To advect the indicator functions, algebraic VOF methods solve the phase transport equation where u r = u 1 − u 2 is the relative velocity, which is a consequence of mass and momentum transfer between the phases. Fleckenstein and Bothe [37] showed that u r may be neglected even in the case of very good solubility (e.g. CO 2 in water) in order to simplify Eq. (16). However, to reduce the smearing of the interface induced by numerical diffusion, an artificial compression term can be introduced by replacing u r in Eq. (16) by a compressive velocity u comp normal to the interface and with an amplitude based on the maximum of the single-field velocity [38] where c α is the compression constant (generally between 0 and 4), Φ f is the volumetric flux across a grid cell face f , and A f is the face area. In all our simulations, we choose c α = 1.0. In addition toṁ 12 , which will be calculated in the next section, the system requires models for the normal vector to the fluid/fluid interface and the surface tension force for closure. Brackbill [36] developed an approximation referred to as the Continuous Surface Force (CSF) where n 12 is calculated from α 1 and n 12 δ 12 is approximated by ∇α 1 , so that The VOF-CSF method is attractive because of its simplicity. However, many studies [39,40] have reported the presence of spurious currents in the capillary dominated regime that originate from errors in calculating the normal vector and the curvature of the interface. Spurious currents may be mitigated by a combination of smoothing and sharpening of the indicator functions [41]. Although these modifications of the CSF may reduce the magnitude of spurious currents, they do not fully eliminate them. In addition, they can potentially deteriorate contact line dynamics [42]. For these reasons, we do not apply any modifications of the CSF method in this work. Spurious currents exist in our simulations, but their impact has been shown in our previous work [26,31] to be negligible when compared to analytical solutions. Their impact in more complex geometries has yet to be understood and is a target of future research. However, in the absence of benchmark experimental data it is impossible to quantify their impact and thus for the purposes of this work we assume them to be negligible.
Multiphase flow in pore structures is generally characterised by two dimensionless numbers, the Reynolds number Re = ρ 1 UL/µ 1 and the capillary number Ca = µ 1 U/σ, where U and L are the reference velocity and length in the domain, respectively. Re describes the ratio of inertial to viscous forces and Ca the ratio of viscous to capillary forces. In this work, we concentrate our investigation to low flow rates, i.e in the creeping flow and capillary dominated regime with Re < 1 and Ca < 10 −4 .

Reactive transport model
In a multiphase system, the chemical species can be present in both fluid phases. The conservation equation is satisfied by the total concentration ψ j,p of a primary species j (Eq. 4) in phase p with where J j,p is the total diffusive flux of primary species j in phase p. We assume that the diffusive flux can be modelled using Fick's law where D j,p and D i,p are the molecular diffusion coefficients of the primary and secondary species in phase p. This is true for dilute species in a solvent, such as water, and for species in a pure or binary mixture. Chemical equilibrium in phase p is insured by the law of mass actions (Eq. 2). At the fluid/fluid interface, the jump conditions are given by the continuity of fluxes and chemical potentials, the latter described here by Henry's law [43], where H k is the Henry constant of species k (primary or secondary), while the total mass conservation at the interface is defined as The diffusion coefficients in the aqueous phase and Henry's constants used in this paper are summarized in Table 1. All species exist only in the aqueous phase, except for CO 2 that can also exist in the gas phase. In this case, the gas phase will be assumed to be pure. Therefore, the diffusion coefficient of all species in the non-aqueous phase can be assumed to be 0. In order to solve reactive transport within the VOF method, the transport equations (Eq. (19)) are integrated over a control volume using volume averaging [31], and the boundary conditions (Eq. (21) and Eq. (22)) are used to eliminate surface integrals arising from the divergence theorem [45]. Since the boundary conditions depend on the concentration of primary and secondary species, it is difficult to develop an accurate and stable transport solver for the total concentrations (ψ j ) 1≤j≤N c . Instead, our model solves directly for the concentration of the primary and secondary species, and Ion D (10 −9 m 2 /s) H (no unit) Ion D (10 −9 m 2 /s) H (no unit) Table 1: Diffusion coefficient of ions in water (obtained from [44]).
is based on a sequential non-iterative operator splitting approach [46]. The transport step solves for the single-field concentration of species k (primary or secondary) using the CST method [31]. The transport step solves the single-field transport equation where is the CST flux of species k and is the single-field diffusion coefficient of species k. At the surface of the solid, the boundary condition for the single-field concentration of species k is defined by [25] At the end of the transport step, the rate of mass transfer is calculated by [31] After the transport step is completed, the reaction step is calculated using the phase concentrations (c k,p ) 1≤k≤Nc , using the law of mass action (Eq. (2)) and the mass conservation of the primary species defined as ∂ψ j,p ∂t = 0.
In addition to the Reynolds and capillary numbers, multicomponent multiphase transport in pore structures is generally characterised using the species Péclet numbers P e j = U L/D j . The transport of a species is advection dominated if P e j > 1, and diffusion dominated if P e j <1.

Interface boundary conditions and artificial mass transfer
One of the objectives of this paper is to demonstrate that an accurate modelling of interface boundary conditions, such as carried out in the CST method, is necessary for robust modelling of multiphase reactive transport because without such modelling artificial mass transfer may arise that can critically damage the chemical equilibrium. This is true even when no interface transfer exists and the species remain in the water phase.
For this we will compare the transport model presented in this paper with the simplified transport model described in Aziz et al. [19] which only solves for the concentration of species in water (Eq. 19). This is achieved by setting the diffusion coefficient in the non-aqueous phase to 0. The single-field equation is defined as It is then generally assumed that a sharp interface between c k and α 2 will be obtained due to the absence of diffusion at the fluid/fluid interface. However, there are two sources of interface transfer that are not accounted for in Eq. (31). First, at the interface, 0 ≤ α 1 ≤ 1, so the diffusion coefficient is not 0, even though D k,2 = 0. Second, artificial mass transfer can occur due to the interface compression term in Eq. (16) if no compression is present in Eq. (31) [31]. We thus demonstrate in section 3.2 how these unaccounted-for sources of artificial mass transfer may damage the numerical solution.

Upscaling
Upscaling of multiphase transport in porous media is generally conducted in terms of the Darcy velocity U p , defined using Darcy's law where P p is the average pressure in phase p, K a is the absolute permeability of the domain and k rp is the relative permeability of phase p. Relative permeabilities are often modelled using the Brooks-Corey model [47] where S p is the macro-scale phase saturation, S wc is the critical water saturation, S nar is the residual non-aqueous saturation, k rp,max is the maximum relative permeability of phase p, and n p is the phase Corey index. The phase saturation S p can be calculated from a pore-scale simulation using where the integral is calculated over the whole domain V .
The phase velocity U p is related to the total velocity U T = U 1 + U 2 by the fractional flow function f p , such as U p = f p U T . The fractional flow functions can then be calculated from Darcy's law, and we obtain Multiphase reactive transport in porous media is usually upscaled using an equilibrium model [48], for which the phase saturation S p (Eq. 34) and the phase average concentrations C j,p are defined as and are computed using an equilibrium phase partitioning. To calculate chemical equilibrium between the species in the aqueous phase, species activities are calculated using the phase average concentrations and then the law of mass actions (Eq. (2)) is applied. However, due to the slow nature of molecular diffusion in water (D ∼ 10 −9 m 2 /s) and the variation in interfacial area due to pore-size heterogeneity ( [30]), the phase distribution is often more accurately predicted using a linear transfer model [26], for which the transfer M k (kmol/s) of species k from phase 1 to 2 is calculated as where λ k (m/s) is the mass exchange coefficient and A 12 is the interfacial area between phase 1 and phase 2, which can be calculated as In addition, equilibrium models usually overpredict the chemical reaction rates [49,50]. Instead a mixing-induced reaction rate is often introduced as where k i (kmol/m 3 /s) is the mixing-induced reaction constant and Ω i is the saturation index of reaction i. For example, for reaction i in Eq. (1) we define the saturation index as We will show in Section 4.1 how pore-scale modelling can be applied to calculate mixing reaction rates.

Implementation
The numerical method has been implemented in GeoChemFoam [51], our reactive transport solver based on OpenFOAM ® [52]. The full code can be downloaded from www.julienmaes.com. The standard VOF solver of OpenFOAM ® , so-called interFoam, has been extended for this purpose into another solver called interReactiveTransferFoam. The full solution procedure is presented in Fig. 1.
Solve species transport (Equ. (25) interFoam solves the system formed by Eq. (15), (16) and (12) on a collocated Eulerian grid. A pressure equation is obtained by combining the continuity (Eq. (15)) and momentum (Eq. (12)) equations. These equations are then solved with a predictor-corrector strategy based on the Pressure Implicit Splitting Operator (PISO) algorithm [53]. Three iterations of the PISO loop are used to stabilise the system. An explicit formulation is used to treat the coupling between the phase distribution equation (Eq. (16)) and the pressure equation. This imposes a limit on the time-step size by introducing a capillary wave time scale described by the Brackbill conditions [36].
In interReactiveTransferFoam, the concentration equation (Eq. (25)) is solved sequentially after the PISO loop. The interfacial mass transfer (Eq. (29)) is then computed and re-injected in the continuity (Eq. (15)) and phase equations (Eq. (16)). The space discretization of the convection terms is then performed using the second-order vanLeer scheme [54]. For the compression terms, the interpolation of α d α c is carried out using the interfaceCompression scheme [52]. The diffusion term ∇. (D j ∇c j ) is discretized using the Gauss linear limited corrected scheme, which is second order and conservative. The discretization of the CST flux is performed using the phase upwinding scheme [32]. Finally, the chemical reaction step is solved using Phreeqc [35].

Verification
The multiphase transport solver has previously been validated by comparison with analytical and semi-analytical solution for a range of 1D, 2D and 3D problems [26,30,31]. In particular the calculation of the local volume change induced by interface transfer for a soluble phase has been validated by comparison with the analytical solution for dissolution of a gas phase in water in a 1D domain. In this study, we present the validation of the coupling between the multiphase transport and chemical reactions.

Multiphase reactive transport in 1D at equilibrium
The objective of this test case is to validate the coupling between multiphase transport and chemical reactions by comparison with a system where an analytical solution exists. For this, we consider a system with 3 components (A, B, and AB) and two phases (water and gas).
We assume that for this case that all activity coefficient γ k = 1.0. Therefore, the law of mass action can be written as where K = 10.0 is the equilibrium constant of the reaction (Eq. (41)). The domain is a 1D tube of 1mm length (Fig. 2). The gas/liquid interface is initially positioned at a distance l 0 = 0.5 mm from the left boundary. The left boundary has a constant pressure p = p 0 , with constant concentration c A,w = 0, c B,w = ρg HAMA and c AB,w = 0, while the right boundary has a no-flow condition.
Since the right boundary has a no-flow condition, and because the fluids are assumed incompressible, the velocity in the gas phase is equal to 0. Hence, the total mass conservation at the interface (Eq. (23)) can be written as which leads to u w ≈ w. Assuming that advective transport is negligible by comparison to diffusive transport, i.e.
the transport equation (Eq. (19) can be considered to be at equilibrium at the time-scale of interface displacement. Therefore Since K >> 1, c AB,w << c A,w and c AB,w << c B,w , an approximated analytical solution for the concentration in the water phase is As only the component A crosses the interface, which shows that Eq. The test case is simulated on a regular grid with 1000 grid blocks, and with a constant time-step t=0.01 s. In order to compare with the analytical solution, the local volume change is initially turned off and the concentration of A in the gas phase is kept equal to 1 kmol/m 3 until the concentrations in the water phase reaches an equilibrium. Local volume change is then turned on and the simulation is run until t=1000 s. Figure 3 show a comparison between simulated and analytical results. We obtain a very good agreement between the model and the analytical solution, and have thus validated the coupling between multiphase flow with interface transfer and chemical reactions in our model.

Injection of a CaCl solution in an oil-filled tube in 2D
The objective of this test case is to show that, unlike the CST method, the simplified model (Eq. (31)) generates artificial mass transfer that damages the numerical solution. First, we consider a 2D straight microchannel of size 300 µm× 100 µm. The fluid properties are summarized in Table  2. The channel is initially filled with oil. At t=0, we start injecting an aqueous solution of 1000 mg/L of CaCl from the left boundary at velocity U = 3 mm/s, which corresponds to Re = 0.3 and Ca = 10 −4 . The solid boundaries are assumed to be oil-wet, with a contact angle of 45 o . In addition, surface complexation occurs at the surface of the solid following the Na-montmorillonite SCM proposed by Bradbury and Bayens [55], which is summarised in  >SOH 0 +Ca 2+ ⇔ >SOCa + + H + 10 −5.9 Table 3: Surface-complexation reactions and their intrinsic stability constant on a clay surface [55].
The aqueous solution includes 4 dilute species (Ca +2 , Cl − , H + and OH − ). Each of these species only exists in the water phase, so that H k = 0 and D k,2 = 0. The diffusion coefficient of species in the water phase are obtained from Li and Gregory [44] and are summarised in Table 1. The transport of these species in the domain is strongly advection-dominated, with Péclet numbers varying from 10.2 to 127.
We assume that the surface of the solid has been previously equilibrated with the same solution of 1000 mg/L of CaCl. Therefore, the chemical equilibrium should be unchanged and the concentration in the water phase constant.
The simulations are performed on a 150×50 cartesian grid with a constant time-step ∆t = 0.5 ms. Figure 4 show the concentration maps for Ca +2 and H + obtained with each method at t=0.15 s. We see that the CST method leads to a sharp interface between species concentration and oil phase CST Simplified Model fraction, with constant concentration in the aqueous phase. No artificial mass transfer occur and the system remains at chemical equilibrium. However, the simplified method leads to a large amount of artificial mass transfer. The species concentrations in the aqueous phase appear diffused and we obtain significant concentration in the oil phase that is purely induced by numerical errors. Note that the simplified model only considers the concentration in the aqueous phase, so the error of concentration in the oil phase can be ignored. However, as a result of the concentration diffusion in the water phase, the chemical equilibrium is disturbed and the concentration of surface species on the solid boundary changes. Figure 5 shows the concentration of >SOCa + along the x-axis at t=0.15 s. We observe that the CST method leads to a constant concentration with no change of concentration by chemical reaction, while the simplified model has a decrease of 5% of >SOCa + across the interface, indicating that changes of concentration by chemical reaction have occurred.
This example demonstrates that the CST method rather than the simplified model should be employed to simulate multicomponent reactive transport in pore-scale images. Additionally, the CST method only requires the computational of two additional fluxes (species compression and CST fluxes), so the increase in CPU time is very limited. For the case presented here, the simplified model ran for 6067 s with 2 processors on an intel Xeon core, while the CST method ran for 6127 s, representing an increase in computational expense of 1%.

Applications
In this section we show how GeoChemFoam can be used to simulate and upscale various reactive processes in pore-scale geometries.
4.1 Test Case 1: CO 2 gas dissolution in a 3D pore cavity In this example, we investigate interface transfer and chemical reactions during dissolution of a CO 2 gas bubble in a pore cavity. The model domain is the same as presented in [32]. The geometry is a 6mm×1mm×1mm channel, with a 2mm×2mm×1mm cavity inserted in the middle (Fig. 6). The domain is meshed using a uniform grid with resolution 50 microns. Initially, CO 2 gas is trapped in the cavity and the rest is filled with water. The fluid properties are summarized in Table 4.   Table 1.
The system includes three chemical reactions that are summarized in Table 5. As CO 2 dissolves in the water phase, H + and HCO − 3 are created and the chemical equilibrium is modified, leading to a decrease in pH.  At t=0, we inject pure water at pH=7 from the left boundary at a flow rate of 0.1 mL/min which corresponds to a capillary number of 3.3×10 −6 . The simulation is run until t=3 min with a constant time-step ∆t = 20 µs. Fig. 7 show the concentration map of CO 2 , OH − and HCO − 3 at the mid-plane at t= 1 min, 2 min and 3 min. The gas/water interface is shown in white. The concentrations are shown with a color map on a log scale to enhance the contrast. We observe that the mixing of species in the water phase is poor. This is because, even though the flow rate is low with Re and Ca well into the creeping and capillary dominated regime, the transport of species is still advection-dominated. For example, the Péclet number for the CO 2 species is equal to 104. Therefore, there is a strong difference between the concentrations upstream and downstream of the cavity. From the inlet and up to the cavity, pH is close to 7 with no CO 2 present. Within the cavity, the water on top of the gas bubble has a pH close to 4 and a CO 2 concentration close to 0.03 kmol/m 3 . From the end of the cavity to the outlet, the pH is close to 5 and CO 2 is present at the bottom of the channel with  : Concentration map of CO 2 , OH − and HCO − 3 at the mid-plan during dissolution of a CO 2 bubble in a 3D pore cavity at t= 1 min, 2 min and 3 min. The gas/water interface is shown in white and the concentration are shown with a color map on a log scale to enhance the contrasts a concentration close to 0.004 kmol/m 3 , but no CO 2 is present in the top part of the channel.
Poor mixing has a strong impact when upscaling the chemical reactions to the larger scales. In Fig. 10, the evolution of gas saturation as well as the concentrations of CO 2 , OH − and HCO − 3 obtained in the pore-scale simulation are compared with the results obtained when using a fullymixed equilibrium model. The results diverge significantly as the concentrations in the equilibrium model trend in the opposite direction to those of the pore-scale simulation. This divergence occurs with the equilibrium model because the CO 2 dissolves instantaneously in the water phase, forming a carbonic acid that significantly reduces the pH of the water, and then the acid is slowly flushed out of the domain and the water becomes increasingly neutral.
However, in reality, the phase transfer occurs on a much larger time-scale (Fig. 10a) and thus a linear transfer model would be more appropriate to simulate this at the larger scale. Using Equ. (37), the mass exchange coefficient for CO 2 can be calculated from the results of the pore-scale simulation as .
The mass exchange coefficient is plotted as a function of the gas saturation S 2 in Fig. 8 and we observe that it can be approximated as a linear function of S 2  In addition, the incomplete mixing in the water phase induces a delay in the chemical reactions and the phase average concentration of species in the domain are not at chemical equilibrium. Mixing-induced reaction rates can be calculated during the pore-scale simulation by integrating the changes of concentrations obtained by chemical reaction (calculated by Phreeqc) over the full simulation domain. Fig. 9 shows the evolution of the reaction rates of the three reactions present in the system (Table 5). We observe that the rates of reactions 1 and 2 converge toward a plateau, which is typical of a mixing-induced reaction constant that does not depend on saturation. However, the rate of reaction 3 consistently decreases from t=0.5 min, which suggests that the mixing-induced reaction constant k 3 decreases as the gas saturation increases. The saturation indexes Ω 1 , Ω 2 and Ω 3 are calculated based on the averaged concentrations in the water obtained from the pore-scale simulation, and the mixing-induced reaction rates are calculated with constant k 1 = 1.80 × 10 −11 and k 2 = 1.04 × 10 −12 kmol/m 3 /s. These along with k 3 = 5.30×10 −7 1−α2 kmol/m 3 /s are plotted in Fig.  9 and compared with the rates obtained from the pore-scale simulation results. We observe that the mixing-induced rates are well fitted to the pore-scale simulation results after an initialisation time of about 0.5 min. These mixing-induced rates are included in the linear transfer model, and the concentration of CO 2 , OH − and HCO 3 obtained are plotted in Fig. 10 and compared to the results

Test Case 2: Injection of a CaCl solution in a micro-CT image
We now investigate multiphase multicomponent reactive transport in a micro-CT image. First we simulate aqueous CaCl injection into an oil saturated pore space with surface complexation. Then from the pore scale result we calculate volume averaged saturation and concentration and compared it to the result an upscaled equilibrium model. Then we propose a correction to the upscaled equilibrium model based on a reduced surface charge.
The image is a 1000 3 voxel micro-CT image of Bentheimer sandstone with a resolution of 2.5 microns, which can be downloaded from the Digital Rock Portal https://dx.doi.org/10.17612/f4h1-w124. A 512 3 voxel image is extracted from the center of the image for this example.
The domain is meshed using the OpenFoam ® snappyHexMesh utility [52]. First a 128 3 cartesian grid is generated. Next, each grid block that is crossed by the solid surface is refined once in each direction, leading to resolution of 5 microns. The cells in the solid phase are then removed, while the cells that intersect the rock/pore interface are replaced by hexahedral or tetrahedral cells that match the solid boundaries. calculated from the mesh and the absolute permeability K a can be estimated by solving the Stokes equation [56]. Our image has a porosity of 0.22 and a permeability of 2.9 × 10 −12 m 2 . The fluid properties (Table 2) and chemical system (Table 3) are the same as the ones used in Section 3.2. The pore space is initially filled with oil and the surface of the solid has been previously equilibrated with a solution of 1000 mg/L of CaCl. At t=0, we inject from the left boundary with a solution of 100 mg/L of CaCl at constant velocity U =3mm/s, corresponding to a capillary number Ca = 10 −4 . A constant pressure is set at the right boundary, while the top, bottom, front and back boundaries have a no-flow condition. The solid boundaries are assumed to be oil-wet, with a contact angle of 45 o . The simulation is run until t=0.5 s with a constant time-step ∆t = 1 µs with 24 processors on an intel Xeon core. The total CPU time of this simulation was 31 days. Fig. 11 shows the water phase fraction, the concentration of H + in the water in the bulk phase, and the concentration of >SOCa + and >SO − on the solid surfaces at t=0.5 s. Although the mixing of H + is not complete, it is better than the mixing in the previous test case, with most values of H + concentration close to 4×10 −8 kmol/m 3 . However, the mixing on the solid surface is very poor. The fractional flow of water at the outlet is calculated from the pore-scale results and plotted in Fig. 12. The curve is fitted to a Brooks-Corey model where k r1,max = k r2,max = 1.0, S wc = 0.24, S nar = 0.25, n 1 = 2, and n 2 = 3, which is also plotted on Fig. 12. Fractional flow is used in an upscaled model to calculate the evolution of the total water saturation in the domain. The results are plotted in Fig. 13a along with the water saturation obtained with the pore-scale simulation. The upscaled model fits the pore-scale simulation with a high degree of accuracy.
We then run a reactive transport model using an upscaled equilibrium model, where the mass action laws (Eq. (2) are calculated using the average concentration of solution species in the water phase and the average concentration of surface species on the solid surface. The results are plotted in Fig. 13b, c and d, along with the concentrations obtained in the pore-scale simulation. We observe that the equilibrium model predicts a higher concentration of >SOCa + and a lower concentration of >SO − . This suggests that the equilibrium model does not overpredict the reaction rate, like in the previous case, but underpredicts it. Therefore, the model cannot be improved by defining mixing-induced reaction rates. Instead, the chemical equilibrium itself should be modified. Since the model predicts a higher concentration of >SOCa + and a lower concentration of >SO − , the surface charge in the equilibrium model is lower in absolute value than the one obtained in the pore-scale model. This discrepancy will have a large impact on the chemical equilibrium as the equilibrium constant depends strongly on the surface charge through the surface potential (Eq. (7)).
In order to obtain a more accurate prediction, the model is corrected by multiplying the surface charge q (Eq. (5)) by 0.95 before calculating the surface potential (Eq. (6)). The results of the  corrected model are plotted in Fig. 13. The corrected model gives significantly more accurate results than the initial upscaled model. However, the errors in the surface concentrations are increasing and the concentration of H + in the bulk remains significantly lower than the one obtained in the pore-scale simulation. This suggests that the model could be further improved by defining mixinginduced reaction rates with the corrected equilibrium constant.

Conclusion
In this study, we presented a novel multiphase reactive transport model to perform direct numerical simulation of multiphase flow, multicomponent transport and geochemical reactions on pore space images. We built a model GeoChemFoam, which is based on OpenFOAM ® [52], an established library to solve partial differential equations, and Phreeqc [35], the most prevalent geochemical solver. The multiphase flow was solved using the VOF method [21], and the transport of species using the CST method [26]. The reactive transport solver was based on a sequential non-iterative operator splitting approach [46] and the chemical equilibrium was solved with Phreeqc. The method was validated successfully for simple configurations where analytical solutions exist. In particular, we showed that the CST method provides an accurate representation of interface boundary conditions free of artificial mass transfer, and it can therefore be applied to model reactive transport in multiphase systems.
We then applied the model to two test cases. In test case 1, we simulated reactive transport during dissolution of a CO 2 gas bubble in a 3D pore cavity. The liquid/gas interface was tracked as well as the concentration of each reactive species in the domain and incomplete mixing was observed. We showed that an upscaled model based on phase and chemical equilibrium could not predict accurately the evolution of average phase saturations and species concentrations in the domain. Instead, the total flux of interface transfer and the average reaction rates in the domain were calculated and we showed that an upscaled model based on linear transfer and mixinginduced reaction rates could accurately predict the evolution of average phase saturations and species concentrations in the domain.
Finally, in test case 2 the model was applied to simulate multiphase reactive transport in a micro-CT image of Bentheimer sandstone where a solution of CaCl was injected into an oil saturated domain with surface complexation at the solid surface. The concentration map of each species on the solid surface was calculated and we observed a poor mixing of charge on the surface. We then ran an upscaled model based on chemical equilibrium and observed that it was overpredicting the change of surface concentration by chemical reactions. Thus we show that surface concentrations cannot be modelled by mixing-induced reaction rates, and the chemical equilibrium need to be modified to take these into account. We then demonstrated that a corrected model that multiply the total surface charge by 0.95 was giving a significantly more accurate result.
The work presented in this paper has wide ranging applications in the oil and gas, carbon capture and storage, contaminant transport, battery, and fuel cell industries. Our simulation framework together with the upscaling methodologies proposed in this paper are an important step forward in our objective of fully characterizing multiphase reactive transport in porous media. Furthermore, this model enables the use of sensitivity analysis to understand how upscaled properties such as the mass exchange coefficient and mixing-induced reaction rates can change with respect to system properties such as injected flow rate or pore-size distribution. In addition, this numerical model can now be bootstraped to field scale multiphase reactive transport simulators using machine-learning regression models by extending work already done for single-phase flow and transport [57] with the ultimate goal of developing upscaling strategies that do not require pore-scale simulations [58].