Reactive Transport Modeling with a Coupled OpenFOAM®-PHREEQC Platform

Upscaling geochemical reaction rates using pore-scale network modeling

The locations of a significant number of industrial facilities, landfills, and almost all mineral ore bodies are characterized by high in situ stresses and fractures, which act as flow paths for fluids underground. Toxic chemical transport is one of the main reasons for loss in groundwater quality and reactive transport modelling is an important tool necessary for predicting and controlling their fate. Modelling reactive transport has largely been achieved using the hydrodynamic transport equation, with an additional reaction and/or decay term. To date, however, most approaches to modelling reactive transport often consider an already existing mixture (multi-phase and/or multi-species) undergoing such transport in homogeneous media. The motivation of this research is to use the transport solutions of two single-phase species in homogeneous media to describe the reactive transport of two species before and after mixing, due to formation heterogeneity. This practical scenario can be observed where regional-scale fracture systems transporting pollution from spatially isolated source locations cause the mixing of chemical pollutants from different source origins due to fracture–fracture flux across two or more intersecting fractures. Here, the authors present a framework for determining the exact solution of the reactive groundwater transport model by the coupling of the homogeneous and non-homogeneous model sub-components. The solution of the homogeneous sub-component is obtained using the Laplace transform. The exact solution of the proposed non-homogeneous equation is obtained analytically using Green’s function method. Further, the two equations are numerically discretized using the Crank–Nicolson scheme and their stability conditions are also established. Instead of presenting the specific solution of a particular case study scenario, a general model to approximate the contribution of reactive processes in the system is presented. It is found that the fate and contribution of each species within the system is captured by the model for which the reaction kinetics and equilibrium partitioning between liquid and solid phases are all assumed to be linear first order.

Coupled three-phase flow and reactive transport processes are widely encountered in many scientific and engineering problems. In the present study, a pore-scale model based on the lattice Boltzmann method is developed to simulate coupled three-phase flow and reactive transport processes. The model is validated by contact angle test of droplets on a curved surface and confined reactive mass transport in a three-phase system. The pore-scale model validated is then employed to study the three-phase reactive transport in channels and porous media. The evolution of the three-phase distribution, the concentration field, and the contact line length are discussed in detail. For a two-channel structure, the result shows that as the viscosity ratio increases, the phase with higher viscosity is more difficult to be displaced. Moreover, as the surface tension force between two certain phases increases, the third phase tends to form a film between the two phases, thus suppressing the reactive transport between the two phases. Finally, pore-scale simulation results of three-phase flow in a two-dimensional porous medium show that as viscosity of the phase to be displaced increases, the recovery rate of the displaced phase decreases, and the displacing phase tends to follow the mechanism of viscous fingering. Finally, while the viscosity of the displaced phase can be reduced due to the existence of the species, the recovery rate does not necessarily increase and sometimes even reduces due to the combined bypass and lubrication effects.

In subsurface rocks, the brine can act as a carrier where non-hydrocarbon components, either natively present or introduced externally, are transported and possibly react with each other as well as with the porous medium. The proposed paper describes the integration of PHREEQC, a leading open-source software for 0D/1D geochemical modeling, with a commercial GPU-based reservoir simulator. This will provide the reservoir engineering community with a user-friendly tool incorporating reactive transport into reservoir studies. PHREEQC is widely used for batch geochemical modeling, thanks to its proven robustness as well as compatibility with extensive databases. We present a non-intrusive approach to enabling reactive transport in reservoir simulation by sequentially coupling a GPU-based reservoir simulator with PHREEQC. The former solves reservoir equations and the transport of chemical master species as aqueous tracers, and the latter computes solutes and rock composition at the end of each time step, based on the master species content in each cell. Then porosity, absolute and relative permeability are explicitly updated by the reservoir simulator as a function of the geochemical state. In this work, two application examples of reactive transport modeling (RTM) using the proposed approach are shown; first RTM is used to simulate in a mechanistic manner a low-salinity EOR process, and second to study injection impairment due to scale deposition. Note that chemical reactions of hydrocarbon/light components dissolved in the aqueous phase are not supported in the current implementation. First, results include an assessment of the method accuracy, as well as a discussion of practical advantages compared to alternative coupling approaches; because all relevant chemical reactions information is kept within PHREEQC datafiles, seamless integration between reservoir and geochemical studies is ensured. Second, results include a set of performance benchmarks to assess the computational cost of such explicit coupling, in particular considering that while the simulator runs on GPUs, chemical reactions are computed on CPUs. This paper shows the first steps towards the integration of RTM capabilities in a commercial GPU reservoir simulator, with the objective of developing the ability to account for advanced chemical processes required for both hydrocarbon recovery and energy transition problems, minimizing the impact on simulator infrastructure and leveraging capabilities of state-of-the-art CPU-based geochemical codes.

<p>Dissolution and precipitation processes in reactive transport in porous rocks play an important role in many contexts, such as geological CO2 storage, reactive contaminant transport, and acid injection in petroleum reservoirs. They are responsible for wormholing and alteration of the rock and transport characteristics due to feedback between the geochemical and the transport processes. A critical aspect of studying reactive transport is the influence of Peclet number on a coupled reactive process, which aims to understand the contributions of advection and diffusion as the two main transport mechanisms.</p><p>Our study investigates the influence of Peclet number on dissolution and precipitation processes in a porous calcite matrix due to reaction with the incoming low-pH flow using a particle tracking (Lagrangian) approach. The coupled reactive process is simulated in a series of computational analyses that capture the subtleties of the multiple scale heterogeneity phenomena, such as anomalous (non-Fickian) transport.</p><p>Our results show that reaction is manifested most significantly for small Peclet numbers, thus signifying the importance of diffusion in mixing processes that facilitate the reaction and increase the heterogeneity of the porous media. The evolution of the field heterogeneity due to reaction follows a similar trait, which is indirectly supported by fractal dimension estimation. Related to the field heterogeneity is the Shannon entropy, which is considered as a measure of the self-organization of the flow.</p>

Analytical solutions of partial differential equation (PDE) models describing reactive transport phenomena in saturated porous media are often used as screening tools to provide insight into contaminant fate and transport processes. While many practical modelling scenarios involve spatially variable coefficients, such as spatially variable flow velocity, v(x), or spatially variable decay rate, k(x), most analytical models deal with constant coefficients. Here we present a framework for constructing exact solutions of PDE models of reactive transport. Our approach is relevant for advection-dominant problems, and is based on a regular perturbation technique. We present a description of the solution technique for a range of one-dimensional scenarios involving constant and variable coefficients, and we show that the solutions compare well with numerical approximations. Our general approach applies to a range of initial conditions and various forms of v(x) and k(x). Instead of simply documenting specific solutions for particular cases, we present a symbolic worksheet, as supplementary material, which enables the solution to be evaluated for different choices of the initial condition, v(x) and k(x). We also discuss how the technique generalizes to apply to models of coupled multispecies reactive transport as well as higher dimensional problems.

Mixing-induced reactions are an essential feature of environmental flow and transport processes. They control many reactive transport processes, including mineral precipitation rates and contaminant remediation processes. Natural porous media are characterized by a strong structural heterogeneity, which impacts solute mixing and, therefore, the resulting chemical reaction rates. Establishing a quantitative link between pore-scale heterogeneity and mixing/reaction rates in saturated and unsaturated conditions remains an open question. Here, we study pore-scale solute mixing using high-resolution experimental measurements to quantify the overall reaction rates and product concentrations. Our goals are to study the impact of structural heterogeneity on 1) reaction rates and products during saturated flow and 2) the spatial arrangement of fluid phases during unsaturated flow and its impact on reaction rates and products. We use two-dimensional porous media consisting of circular posts in a Hele-Shaw-type flow cell. We control heterogeneity by varying the posts’ diameters disorder and correlation length; increasing this length introduces more structure in the porous medium. We utilize an irreversible oxidation reaction to produce fluorescein from its non-fluorescent form. The Damköhler number is sufficiently larger than unity, so the reaction rate is mixing-controlled. We inject a non-fluorescent tracer pulse into the porous medium sample filled with the oxidating reactant under saturated and unsaturated flow conditions. We analyze periodic fluorescence intensity images to track the evolving solute concentration field. The reaction rates and the total reaction product mass are calculated directly from the concentration images. Solute concentration images show that increasing the spatial correlation length under saturated flow conditions leads to enhanced reaction front stretching and elongation as the solute travels along preferential pathways. Due to this overall stretching, the reaction front is locally more compressed perpendicular to the elongation direction. In a non-correlated, randomly disordered porous medium, overall stretching is reduced, and the front is less compressed locally. Under unsaturated flow conditions, a main preferential flow path characterizes the correlated porous medium. In contrast, the non-correlated medium is characterized by a higher degree of branching and splitting in the velocity field. Solute pulse focusing in the correlated porous medium sample reduces reaction front stretching compared to the non-correlated porous medium, under unsaturated conditions. Under these conditions, the reaction rate increases more than the saturated case due to the unsaturated flow pattern's enhanced reaction front stretching. This effect is more pronounced for the non-correlated sample, where flow path splitting and reaction front stretching are more significant. This work shows that structural heterogeneity has a considerable effect on reactive solute transport and that this effect depends on the system’s saturation.

Reactive transport modeling involves solving several nonlinear coupled phenomena, among them, the flow of fluid phases, the transport of chemical species and energy, and chemical reactions. There are different ways to consider this coupling that might be more or less suitable depending on the nature of the problem to be solved. In this paper we acknowledge the importance of flexibility on reactive transport codes and how object-oriented programming can facilitate this feature. We present PROOST, an object-oriented code that allows solving reactive transport problems considering different coupling approaches. The code main classes and their interactions are presented. PROOST performance is illustrated by the resolution of a multiphase reactive transport problem where geochemistry affects hydrodynamic processes.

Acid gas generation by thermochemical sulfate reduction (TSR) evolves within a complex web of petroleum-water-rock-gas interactions in reservoirs under high temperature conditions of more than ca. 100°C. The interactions lead to the formation of toxic and corrosive hydrogen sulfide (free H2S gas and dissolved H2S). Such interactions are caused by the instability of hydrocarbons in the presence of water and a reactive reservoir rock matrix containing water-soluble anhydrite. The mass conversions of the inorganic water-rock-gas interactions which are triggered by the redox degradation of hydrocarbons establish a certain, thermodynamically defined state of chemical equilibrium. Any approach to geochemically model “acid gas generation” and “H2S-risk distribution” in petroleum systems should be based on a conceptual model that adequately reproduces the intimately interconnected and interdependent nature of all isochronous hydrogeochemical reactions, whether they are kinetically controlled or establish equilibrium species distributions. Such approaches rely (1) on the thermodynamical calculation of chemical equilibrium species distribution, (2) on the coupling of kinetically controlled oil degradation and sulfate reduction by oil-derived reductants to the equilibrium calculations, and (3) on the calculation of diffusive mass transport through the free pore water network and the irreducible water film. The key to model TSR, “acid gas generation, and “H2S-risk distribution” is not to consider and model any single, isolated reaction like the kinetically controlled sulfate reduction which depends on the thermal history. The actual key to model TSR, the fate and behavior of sulfidic sulfur, and a realistic “H2S-risk distribution” in petroleum reservoirs is an overall reproduction of the hydrogeochemical reactive transport processes which temporally and spatially evolve in a complex network of oil/petroleum-water-rock-gas interactions under reservoir conditions. Consequently, we perform 3D hydrogeochemical, multi-component and multi-species reactive mass transport modeling for a semi-generic case study by using the PHAST computer code (provided by the U.S. Geological Survey) and take the following boundary conditions into account: gas reservoir; carbonate (dolomite plus calcite) reservoir rocks; anhydrite seal; 140°C; 600 bar total pressure; kinetic rate constant for sulfate reduction by CH4 = 1.08 x 10-16 mol s-1 l-1; mass transport is restricted to diffusion; modeled time span is 10 Ma

Reactive transport simulations are a common approach for the quantitative assessment of contaminant biodegradation in the subsurface. To use knowledge on microbial kinetics for the simulation of in situ biodegradation, the mass transfer processes controlling the bioavailability of the contaminants need to be described appropriately. A common approach to account for this problem is using a linear exchange model controlling the link between bulk and bioavailable concentration. Assuming that the subsequent degradation is controlled by the bioavailable concentration, only, these two steps can be combined to an analytical expression for the overall reaction rate know as the Best-Equation. In our work, we evaluate this approach for its ability to describe biodegradation kinetics limited by pore-scale mass transfer. Results from explicit numerical and analytical simulations of mass transport and reactive consumption at the pore scale are used to test the accuracy of results obtained using the Best-Equation. Our analysis shows that strictly spoken the Best-Equation is not valid. However, a good approximation can be achieved with errors of less than 6% in cases of moderate bioavailability and much lower errors in cases of either low or high bioavailability. These results support the description of mass transfer processes used in many reactive transport models. Furthermore, we present a method to obtain an accurate estimate of the unknown rate parameter controlling the diffusive mass transfer processes at the pore scale.

Abstract. Many types of geologic subsurface utilisation are associated with fluid and heat flow as well as simultaneously occurring chemical reactions. For that reason, reactive transport models are required to understand and reproduce the governing processes. In this regard, reactive transport codes must be highly flexible to cover a wide range of applications, while being applicable by users without extensive programming skills at the same time. In this context, we present an extension of the Open Source and Open Access TRANSPORT Simulation Environment, which has been coupled with the geochemical reaction module PHREEQC, and thus provides multiple new features that make it applicable to complex reactive transport problems in various geoscientific fields. Code readability is ensured by the applied high-level programming language Python which is relatively easy to learn compared to low-level programming languages such as C, C++ and FORTRAN. Thus, also users with limited software development knowledge can benefit from the presented simulation environment due to the low entry-level programming skill requirements. In the present study, common geochemical benchmarks are used to verify the numerical code implementation. Currently, the coupled simulator can be used to investigate 3D single-phase fluid and heat flow as well as multicomponent solute transport in porous media. In addition to that, a wide range of equilibrium and nonequilibrium reactions can be considered. Chemical feedback on fluid flow is provided by adapting porosity and permeability of the porous media as well as fluid properties. Thereby, users are in full control of the underlying functions in terms of fluid and rock equations of state, coupled geochemical modules used for reactive transport, dynamic boundary conditions and mass balance calculations. Both, the solution of the system of partial differential equations and the PHREEQC module, can be easily parallelised to increase computational efficiency. The benchmarks used in the present study include density-driven flow as well as advective, diffusive and dispersive reactive transport of solutes. Furthermore, porosity and permeability changes caused by kinetically controlled dissolution-precipitation reactions are considered to verify the main features of our reactive transport code. In future, the code implementation can be used to quantify processes encountered in different types of subsurface utilisation, such as water resource management as well as geothermal energy production, as well as geological energy, CO2 and nuclear waste storage.

In multi‐mineral fractured rocks, the altered porous layer on the fracture surface resulting from preferential dissolution of the fast‐reacting minerals can have profound impacts on subsequent chemical‐physical alteration of the fractures. This study adopts the micro‐continuum approach to provide further understanding of reactive transport processes in the altered layer (AL), and mass exchanges with the bordering matrix and fracture. The modeling framework couples the Darcy‐Brinkman‐Stokes (DBS) solver in COMSOL Multiphysics and the geochemical modeling capability of CrunchFlow. Three‐dimensional steady state simulations with systematically varied chemical‐physical parameters of the AL were performed to examine the impacts of individual factors and processes. Our simulation results confirm previous observations that dissolution of the fast‐reacting mineral (i.e., calcite) is largely controlled by diffusion across the AL. We also show that dissolution of the slow‐reacting mineral (i.e., dolomite), which controls AL development and fracture enlargement, increases with surface area and has a complex dependence on different local rate‐limiting processes. In particular, advection can result in evident spatial variations in the local dissolution rates of dolomite, although it does not affect the bulk chemistry significantly. The difference in the spatial patterns between simulations with and without advection in the AL is more noticeable in the locations with smaller apertures, with up to 20% difference in local reaction rates. Therefore, it is important to include a full depiction of advection, diffusion, and reactions for accurately capturing local dynamics that control long‐term fracture evolution.

Numerical investigation of pore and continuum scale formulations of bimolecular reactive transport in porous media

<p>The mining industry globally produces millions of tons of waste rock every year. The weathering of exposed metal(loid)-rich waste rock can produce poor-quality effluent, and mine sites therefore need to establish water-quality management strategies that predict and mitigate environmental impacts. Technical frameworks to support drainage quality predictions and industrial waste-rock management typically combine classic static and kinetic testing procedures, field-scale experiments and sometimes geochemical equilibrium- and reactive-transport models. However, predictions of waste rock weathering and drainage processes remain challenging on relevant spatiotemporal scales, due to site-specificity in waste rock and local weathering conditions, unresolved heterogeneity in large waste-rock systems and the intricate (non-linear) coupling between chemical kinetics and mass- and heat transfer processes.</p><p>We synthesized long-term (>10 yr) hydrogeochemical field data from a multiscale experimental research program at the Antamina mine, Peru. At Antamina, various waste-rock materials have been extensively hydraulically, physically and geochemically characterized and weathered at different spatiotemporal scales. This data set provides a unique opportunity to quantitatively assess the mechanisms that affect drainage from different waste-rock types under field conditions. Monitoring of weathering rates in humidity cell tests (~1 kg), column experiments (~170 kg), field barrel kinetic tests (~350 kg), and mesoscale experimental piles (~6,500,000 kg) revealed that normalized mass loadings from different waste-rock types systematically decreased with increasing experimental scale.</p><p>We developed a process-based reactive-transport framework to reproduce the recorded waste-rock drainage trends from the various field experiments. For each of the experiments, 1-D reactive-transport models were constructed in MIN3P-HPC, all including the same formulations for, e.g., transient unsaturated flow, advective-diffusive transport of aqueous species, gas diffusion, gas-liquid partitioning and equilibrium or kinetic mineral dissolution and precipitation reactions. The models were exclusively parameterized with measured field hydrostatics (e.g., tracer testing, volumetric water contents; van Genuchten parameters), analyzed physicochemical bulk waste-rock properties (e.g., bulk geochemistry, mineral content, particle size), or adopted literature values (e.g., kinetic rate laws and constants).</p><p>At all experimental scales, the recorded drainage quality evolution could be successfully reproduced with the consistent suite of field-parameterized physical transport processes and kinetic rate laws. A comparison of fitted effective rate coefficients reveals that reduced weathering rates at increasing scales mostly originate from decreasing specific mineral surface areas (particle sizes increase with experimental scale) and possibly by surface passivation, although the effects of flow bypassing and channeling are not yet fully investigated. This work demonstrates that with efforts focused on the identification and parameterization of the relevant physicochemical processes, effective yet process-based models can be developed from readily available bulk waste-rock parameters to predict and upscale mine waste rock weathering and drainage quality trends across laboratory-to-practice-relevant scale ranges.</p>

<p>The modern advances in computing and experimental capabilities in the research of water-rock-interactions require geoscientists to routinely combine laboratory data and models to produce knowledge in order to solve pressing societal challenges connected to subsurface utilization. Data science is hence a more and more pervasive instrument also for  geochemists, which in turn demands flexible and easy to learn software adaptable to their specific needs. <br>In this contribution we showcase geochemical and reactive transport modelling with our RedModRphree [1] extension package for the GNU R environment and programming language. The new version of the package leverages the R interface to the established PHREEQC geochemical simulator maintained by its original authors [2]. R has established itself as de facto standard language for statistics and machine learning. It enjoys increasing diffusion in many applied scientific fields such as bioinformatics, chemometrics and ecological modelling. The availability of excellent third party extensions such as the thermodynamic package CHNOSZ [3], which extends the functionalities of SUPCRT92, as well as its advanced graphical and numerical capabilities, make R an attractive platform for comprehensive geochemical data analysis, experiment evaluation and modelling. <br>The aim of RedModRphree is to provide the user with an easy-to-use, high-level interface to program algorithms involving geochemical models, which are then solved using the PHREEQC engine: parameter calibration, error and sensitivity analysis, visualization, up to CPU-intensive parallel coupled reactive transport models. Among the out-of-the-box features included in RedModRphree, we highlight the computation and visualization of Pourbaix (Eh-pH) diagrams and the implementation of 1D advective reactive transport supporting the use of surrogate models replacing expensive PHREEQC calculations [4]. RedModRphree is open source and can be installed from https://git.gfz-potsdam.de/delucia/RedModRphree.</p><p><img src="https://contentmanager.copernicus.org/fileStorageProxy.php?f=gnp.72552df75cff57113730161/sdaolpUECMynit/12UGE&app=m&a=0&c=afa6bef86f4c5523b81a86ccdd579fc9&ct=x&pn=gnp.elif&d=1" alt=""></p><p>[1] De Lucia, M. and Kühn, M.: Coupling R and PHREEQC: Efficient Programming of Geochemical Models, Energy Procedia, 40, 464–471, doi.org/10.1016/j.egypro.2013.08.053, 2013.</p><p>[2] Charlton, S.R. and Parkhurst, D.L.: Modules based on the geochemical model PHREEQC for use in scripting and programming languages, Computers & Geosciences 37, 10, 1653–1663, doi.org/10.1016/j.cageo.2011.02.005, 2011.</p><p>[3] Dick, J.M.: CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry, Frontiers in Earth Science, 7, https://doi.org/10.3389/feart.2019.00180, 2019.</p><p>[4] Jatnieks, J., De Lucia, M., Dransch, D., and Sips, M.: Data-driven Surrogate Model Approach for Improving the Performance of Reactive Transport Simulations, Energy Procedia, 97, 447–453, doi.org/10.1016/j.egypro.2016.10.047, 2016.</p><p> </p>