Interactions between multiphase flow, geochemical reactions, and solute transport in fractured porous media are ubiquitous in Earth's critical zone and subsurface systems, and affect the dynamics of many environmental and energy engineering applications. Examples include supercritical CO₂ (scCO₂) injection, light non-aqueous phase fluid (LNAPL) contamination, and Enhanced Oil Recovery (EOR). In geologic carbon storage systems, the injected scCO₂ displaces the native brine and can be trapped in small pores or can dissolve into the brine and react with the host rocks. The resulting mineral dissolution and precipitation can modify the porosity and permeability of the reservoirs and caprocks, affecting injectivity of the reservoirs and long-term storage security (Johnson et al., 2001, 2004; Xu et al., 2011b). In EOR, low-salinity water is injected to the oil reservoir to improve sweeping efficiency by modifying wettability via surface complexation reactions (Zhang et al., 2006; Kumar et al., 2011). In the shallow subsurface, LNAPL contamination of groundwater is a widespread environmental problem. Light non-aqueous phase fluid is insoluble in water and typically fluctuates with the water table. As a result, LNAPL is spread vertically following local drainage-imbibition cycles, which in turn affects the degradation and thus the fate of the contaminants (Ngien et al., 2012; Pan et al., 2016; Govindarajan et al., 2018). The interactions between two-phase flow and electrochemical reactions have also been identified as a research priority for the design of effective fuel cells, as gas bubbles (of e.g., H₂ and O₂) can be generated by side reactions during charging and affect power generation (Chen et al., 2017; Grunewald et al., 2021).

Continuum-scale numerical models are widely used to investigate and predict the evolution of fractured porous media caused by multiphase reactive transport processes. Models such as TOUGHREACT, OpenGeoSys, and PFLOTRAN are versatile tools that have been used to investigate scCO₂ migration in fractured formations (Xu et al., 2011a, 2019; Xiao et al., 2020), dissolution, transport and biodegradation of non-aqueous phase fluid (NAPL) in shallow aquifers (Pruess, 2004; Popp et al., 2015; Sookhak Lari et al., 2019), and heat extraction using CO₂ as a working fluid (Lichtner and Karra, 2014). These models however rely heavily on constitutive relations such as the Brooks and Corey equation and the van Genuchten equation that relates the water saturation with capillary pressure and relative permeability.

Pore-scale structural heterogeneity and processes that are not considered explicitly in the continuum description, however, could be important. For instance, microfluidic experiments from Karadimitriou et al. (2016, 2017) examined non-Fickian transport in a water-Fluorinert immiscible fluid system. The results demonstrated that the mass transfer rate between mobile–immobile zones under the two-phase flow condition is not constant, indicating that the conventional treatment of the mass transfer rate in the continuum scale Mobile–Immobile (MIM) model could introduce significant errors in the simulation. Pore-scale experiments from Jiménez-Martínez et al. (2015) using 2D microfluidic cells with homogeneous pore structures have also demonstrated that solute mixing can be significantly enhanced under multiphase flow conditions because of ramified finger flow structure, non-Fickian dispersion, and non-wetting phase clusters that limit the finger flow merging. Recent advancement in experimental techniques have also enabled direct observations of reactive transport in multiphase systems at the pore scale. Using biogenically calcite-functionalized micromodels, Song et al. (2018) observed a new microscale mechanism that affects reactive transport in the CO₂-brine-calcite system. The CO₂ gas phase can accumulate on the mineral surfaces following calcite dissolution, which protects calcite from further dissolution. Jiménez-Martínez et al. (2020) developed a high-pressure geomaterial microfluidic device and investigated mineral reactions in etched channels during co-injection of CO₂-saturated brine and scCO₂. They observed that the presence of scCO₂ bubbles significantly changed both the flow dynamics and the reaction patterns, compared to the single phase flow experiments. In the competing channels with different widths, mineral dissolution was more homogenized and carbonate precipitation was enhanced in the low-velocity regions formed as a result of the presence of the bubbles, which was also confirmed by pore-scale Lattice Boltzmann Methods (LBM) simulations.

Pore-scale models provide an invaluable tool to further our understanding of pore-scale dynamics. Pore-network models (PNMs) are a computationally efficient option for simulating the reactive transport processes under two-phase flow conditions, but require simplifications of the geometric solid-fluid interface and assumptions about uniform aqueous concentrations within a pore (Xiong et al., 2016). Lattice Boltzmann Methods (Chen et al., 2013, 2015) and direct numerical simulations (DNS) (Haroun et al., 2010b; Marschall et al., 2012; Maes and Soulaine, 2018; Soulaine et al., 2018, 2021) provide alternatives that relax these restrictions. The LBM was able to reproduce experimental observations on wettability alteration during Low Salinity Water Flooding (Akai et al., 2020). It was also shown by Maes and Geiger (2018) using a DNS approach that the alternation is driven by the concentration of the potential determining ions (PDI) instead of pH. Moreover, a series of micro-continuum simulations showed that the production of CO₂ bubbles resulting from carbonate dissolution may limit the subsequent dissolution and prevent the emergence of wormholes (Soulaine et al., 2018).

In spite of the growing use of fully-resolved pore-scale models (LBM or DNS) in multiphase reactive applications, they have not been used to understand how reaction rates are influenced by flow processes in the same way as it has been done in single-phase systems (e.g., Molins et al., 2012, 2014; Deng et al., 2018). In single phase systems, reactive surface area has been typically varied to account for flow dynamics and transport limitations. Pore-scale simulations have highlighted the complex dependence of the correction factor on pore-scale geometry and flow regimes (Deng et al., 2018). In multiphase systems, however, this information does not exist. As a result, multiphase continuum scale models rarely account for the effect of the flow dynamics on reaction rates. Further, reaction rates are commonly assumed to be independent of phase saturations (e.g., Xu et al., 2011a; Lichtner et al., 2015; Águila et al., 2020; Wu et al., 2021). Only in some rare instances, the reactive surface area is allowed to vary with liquid saturation. For example, in the active fracture module implemented in TOUGHREACT, the reactive surface area and thus reaction rate follows a power law relation with respect to the water saturation and the exponent is an empirical variable (Sonnenthal et al., 2005). However, there is still a lack of pore-scale studies and mechanistic understanding that support the development of this type of constitutive relationship for considering the impacts of multiphase flow on reaction rates.

Our study aims to bridge this gap, by performing a series of well-designed pore-scale multiphase reactive transport simulations. To this end, a pore-scale multiphase reactive transport modeling framework was developed by coupling OpenFOAM and CrunchTope. A set of simulations were performed with co-injection of air and CO₂-acidified water into 2D calcite channels to examine calcite dissolution rate under a range of flow conditions. A sine wave geometry was used to introduce different levels of roughness in the channels. This simple geometry has been widely used to provide an idealized representation of surface roughness to investigate its impacts on fluid flow and chemical transport (Kitanidis and Dykaar, 1997; Bolster et al., 2009; Sund et al., 2015; Deng et al., 2018). Section Methodology details the modeling framework, the mathematical principles, and the simulation setups. The results of the numerical simulations, including analyses of the flow field and reaction rates, and observed relations between reaction rate and parameters such as saturation are presented in section Results. We discuss the broader implications of our study in section Discussion and conclude in section Summary and Conclusions.

Our approach entails the simulation of two-phase flow and reactive transport in a series of synthetic geometries. For this purpose, we develop a modeling framework by coupling two widely used and thoroughly validated codes. In this section, we present this modeling framework and we describe the setup of the simulations.

In our simulations, the recirculation zone is the dominant hydrodynamic feature affecting the mineral dissolution rate. The development of recirculation zones have been observed in single phase systems at high velocity when the contribution of inertial effect becomes significant; and roughness allows recirculation zones to form under lower velocities (Deng et al., 2018). Our results illustrate that in two phase systems, momentum transfer across the fluid–fluid interface further extends the conditions under which recirculation zones form as also reported in previous studies (Heshmati and Piri, 2018). Other multiphase flow dynamics are also observed in our results albeit the simplified setup considered. For instance, in geometry b, a few simulations at the high and low flow velocities especially with large v_w/v_nw showed coalescence of the gas bubbles, which resulted in significantly lower saturation in the reactive zone. Bubble coalescence and breakup are widely observed phenomena that are dependent on fluid properties, velocities, and geometries (Jo and Revankar, 2009; Paulsen et al., 2014; Chen et al., 2017; Mahabadi et al., 2018; Ren et al., 2020; Grunewald et al., 2021). These processes are accompanied by the re-organization of the fluid–fluid interface and can introduce perturbations in the velocity field. The simulations with gas bubble coalescence were excluded from the analyses in section Correlation Between Reaction Rate Modification and Liquid Saturation as the new hydrodynamics is not directly comparable with the other simulations. Nonetheless, this observation highlights that compared to single phase systems in which a lower velocity typically transfers to a stronger transport limitation, the two phase systems require consideration of additional hydrodynamics that may arise at a different velocity (Blois et al., 2015).

A variety of complex hydrodynamics can arise from the migration of fluid–fluid interfaces depending on the velocity (e.g., capillary numbers) and pore morphology (e.g., Berg et al., 2013; Spurin et al., 2019; Li et al., 2021; Wang et al., 2021).

From a macroscopic perspective, as capillary number (Ca) increases, interface migration transitions from the capillary fingering regime with more random local movement to the viscous fingering with more stable displacement (Toussaint et al., 2012; Li et al., 2019; Grunewald et al., 2021). However, Ca alone does not provide a good description of the fluid–fluid interface dynamics (Armstrong et al., 2015), which is more sensitive to pore-scale roughness/morphology at low capillary numbers (Toussaint et al., 2012).

From a microscopic perspective, both the magnitude and direction of flow were observed to fluctuate before the arrival of the invading front (Roman et al., 2016; Li et al., 2017). Strong instabilities of the interface can also lead to pore-scale burst events such as Haines jumps and can increase local velocities by one–two orders of magnitude (Blois et al., 2015; Li et al., 2017, 2021). These hydrodynamics are reported to promote mixing and thus reactions in the liquid phases (Jiménez-Martínez et al., 2015, 2016, 2017), and may also increase local mineral reaction rates. In the study of Jiménez-Martínez et al. (2020), calcite dissolution rate in the two phase experiment with Ca = ~10⁻⁵ and v_w/v_nw = 1.0 is 68% of the dissolution rate in the single phase experiment. This is comparable to the simulation results in the smooth channel, e.g., the reaction rate ratio is ~64 and ~72% at v_w = 0.1 m/s (i.e., Ca = ~3.5 × 10⁻⁴) for v_w/v_nw of 0.5 and 0.25, respectively (Figure 7). Their reaction rate ratio is slightly higher than what would be expected based on our simulations, which is likely because of the oscillation of gas bubbles in the presence of a system of channels and thus higher local velocity at the fluid–fluid interface, in addition to the continuous supply of CO₂ from the co-injected scCO₂ phase. The competition of different fluid pathways can result in more dynamic bubble migration and flow field rearrangement, which influences self-organization of the system, which is also illustrated in the column experiment of Ott and Oedai (2015) and the pore-scale numerical simulations of Soulaine et al. (2018). Geometry units such as pore-doublets can be used to capture such dynamics and thus offer additional insights (Mohammadi Alamooti et al., 2020; Alizadeh and Fatemi, 2021). More realistic geometries may be needed to fully examine the complexity of real systems, as flow instabilities are primarily driven by geometry when morphological heterogeneity is large (Li et al., 2017; Heshmati and Piri, 2018).

Water film is also an important contributor to the fluid–fluid interface (Li et al., 2019). While for the given Ca and channel width used in our study, water film is not expected to be well developed and thus contribute to transport significantly or affect our analyses (Roman et al., 2017), well-developed water film may become important in maintaining water-rock contact and transport pathways between isolated water parcels. For instance, water film has been observed during the drainage process in various experiments (Rücker et al., 2015; Schlüter et al., 2016; Roman et al., 2017; Moura et al., 2019). It enhances the connectivity between residual water and may cause the “snap-off” phenomenon, affecting the contact between water and the solid phase (e.g., rock) and the transport processes.

Overall, this study represents an early-stage effort in highlighting the importance of a dynamic coupling between reaction rates and multiphase flow dynamics. As such, it focuses on a specific geometric setup and a relatively limited number of flow scenarios. However, compared to the experimental studies that typically report temporally and/or spatially integrated reaction rates, our model provides information on local and instantaneous reaction rate along with detailed velocity field and solute transport, which helps to better explore microscopic mechanisms. The two mechanisms identified in our pore-scale two-phase reactive transport simulations—the surface area effect and local hydrodynamics—support recent observations from micromodel experiments (Song et al., 2018; Jiménez-Martínez et al., 2020). Moreover, the pore-scale perspective in the model enables insights into the impacts of multiphase flow dynamics on mineral reaction rate that have broader implications. Namely, results indicate that in addition to the fluid-solid interface, a direct measure of the reactive surface area, the fluid–fluid interface also plays an important role in controlling solid-phase reactions by modifying local hydrodynamics. As discussed above, local hydrodynamics that can affect mineral reaction rate is largely influenced by fluid–fluid interfaces, our study suggests that adding interfacial area in the formulation of reactive surface area in two-phase flow conditions in addition to saturation may be a logical step. In fact, it has been proposed to use the fluid–fluid interface as a measure of flow topology and thus an additional parameter in constitutive relations of relative permeability (Picchi and Battiato, 2018). Such consideration may be of particular interest in e.g., scCO₂-brine systems where mass transfer across the fluid–fluid interface also affect fluid chemistry and thus mineral reaction rate.

In order to develop constitutive relations that faithfully reflect the coupling between mineral reaction rates and multiphase flow dynamics and are broadly applicable in natural and engineered fractured porous materials, systematic studies that explore a broader range of multiphase flow and geometric conditions are needed. While the modeling framework developed in this work can be adapted to consider these processes and thus provides a suitable modeling tool in a wide range of applications for pore-scale mechanistic investigations, it needs to be acknowledged that these investigations hinge upon further development of modeling capabilities for multiphase flow dynamics (Aboukhedr et al., 2018; Qin et al., 2020). For instance, there is still a lack of comprehensive comparison and benchmarking of pore-scale models for low Ca flow systems (Zhao et al., 2019). It has also been shown that because of the sensitivity of multiphase flow to local perturbations, a deterministic reproduction of the dynamics as observed in experiments with relatively large porous domains using numerical models is challenging, while statistical behaviors can still be captured (Ferrari et al., 2015). The consideration of reactions, however, may require higher spatial accuracy because the spatial variations (in mineral distribution and reaction rate) can be important. The impacts of multiphase hydrodynamics on aqueous reactions and mineral reactions may also need to be considered separately. Taking the recirculation zone as an example, it increases interfacial mass transfer (Maes and Soulaine, 2018), but reduces mineral reactions (Deng et al., 2018). These research needs further emphasize the need for experimental studies that underpin the development of mechanistic models for multiphase flow in reactive systems, and studies that expand our investigations from quasi-2D micromodels to 3D systems. The development of constitutive relations will also require the investigations be extended to larger scales, which calls for development of a complementary approach that leverages modeling and experiments at different scales (Blunt et al., 2013).

A pore-scale multiphase reactive transport model was established and used to investigate the dependence of mineral reaction rate on pore-scale two-phase flow dynamics. A series of numerical simulations were performed, in which CO₂-acidified water and air were co-injected at a range of velocities into 2D calcite channels with different levels of roughness as defined by a single sine wave. The simulation results showed that gas bubbles migrating through the reactive zone as a result of the two-phase co-injection caused the reaction rate to be lower than that of the single phase flow simulation with the same total injection rate. Our analyses revealed that the decrease of the mineral reaction rate is caused by a combination of two mechanisms: the reduction in the wetted surface area and the transport limitations that arise from the two phase flow. In the rough geometries, the narrow throats are reaction hotspots in the single phase flow simulations, whereas these locations are occupied by gas bubbles periodically and thus not accessible for reactions. Meanwhile, the flow field showed clear “vortices,” i.e., recirculation zones, as the gas bubbles migrate through the trough of the sine wave. As a result, more reaction products are trapped in the trough and the local thermodynamic driving force and kinetic rates are reduced. Additionally, the correlations between fluid saturation and the extent of reaction rate reduction—which is measured by the ratio between the reaction rate from the two-phase flow simulation and that from the corresponding single phase simulation—were analyzed to provide insights for continuum-scale modeling. A non-monotonic relationship was observed, which is because the contribution from wetted surface area reduction and local transport limitation show opposite dependence on the roughness of the channel, which controls the saturation for a given flow condition. Through these numerical simulations, we highlighted the complexity of reactive transport in multiphase flow systems and identified two important mechanisms through which mineral reaction rates are affected. These results highlight the need for consideration of interfacial dynamics on mineral reaction rates in multiphase flow systems, and also emphasize the need for experimental studies that underpin the development of mechanistic models for multiphase flow in reactive systems.

The raw data supporting the conclusions of this article will be made available by the authors upon reasonable request, without undue reservation.

PL contributed to developing the code, designing the study, performing the simulations and analyses, and writing the manuscript. HD designed the study, contributed to code development, simulation results analyses, and writing the manuscript. SM contributed to code development and writing the manuscript. All authors contributed to the article and approved the submitted version.

This work is supported by the Laboratory Directed Research and Development (LDRD) award from Berkeley Lab, provided by the Director, Office of Science of the U.S. Department of Energy (DOE), and by the Chemical Sciences, Geosciences, and Biosciences Division, Office of Basic Energy Sciences, DOE, under contract No. DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.