Dissolution–precipitation processes in tank experiments for testing numerical models for reactive transport calculations: Experiments and modelling

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.

The use of the subsurface for CO2 storage, geothermal energy generation, and nuclear waste disposal will greatly increase the interaction between clay(stone) and concrete. The development of models describing the mineralogical transformations at this interface is complicated, because contrasting geochemical conditions (Eh, pH, solution composition, etc.) induce steep concentration gradients and a high mineral reactivity. Due to the complexity of the problem, analytical solutions are not available to verify code accuracy, rendering code intercomparisons as the most efficient method for assessing code capabilities and for building confidence in the used model. A benchmark problem was established for tackling this issue. We summarize three scenarios with increasing geochemical complexity in this paper. The processes considered in the simulations are diffusion-controlled transport in saturated media under isothermal conditions, cation exchange reactions, and both local equilibrium and kinetically controlled mineral dissolution-precipitation reactions. No update of the pore diffusion coefficient as a function of porosity changes was considered. Seven international teams participated in this benchmarking exercise. The reactive transport codes used (TOUGHREACT, PHREEQC, with two different ways of handling transport, CRUNCH, HYTEC, ORCHESTRA, MIN3P-THCm) gave very similar patterns in terms of predicted solute concentrations and mineral distributions. Some differences linked to the considered activity models were observed, but they do not bias the general system evolution. The benchmarking exercise thus demonstrates that a reactive transport modelling specification for long-term performance assessment can be consistently addressed by multiple simulators.

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.

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

<p>The formation of alkaline magmas observed worldwide requires that low degree-melts, potentially formed in the asthenosphere, were able to cross the overlying lithosphere. Fracturing in the upper, brittle part of the lithosphere may help to extract this melt to the surface. However, the mechanism of extraction in the lower, ductile part of the lithosphere is still contentious. Metasomatic enrichment of the lithospheric mantle demonstrates that such low-degree melts interact with the lithosphere, but the physical aspect of this process remains unclear. The aim of this study is to better understand the percolation of magma in a porous viscous medium at pressure (P) and temperature (T) conditions relevant for the base of the lithosphere. We study such melt percolation numerically with a Thermo-Hydro-Chemical model of reactive transport coupled with thermodynamic data obtained via Gibbs energy minimisation. We perform Gibbs energy minimisation with Matlab using the linprog algorithm. We start with a simple ternary system of Forsterite/Fayalite/Enstatite solids and melts. All variables are a function of T, P and composition of the system (C), and are computed in both the Gibbs energy minimisation and in the reactive transport code, and can therefore vary freely.</p>

Although multicomponent reactive transport modeling is gaining wider application in various geoscience fields, it continues to present significant mathematical and computational challenges. There is a need to solve and compare the solutions to complex benchmark problems, using a variety of codes, because such intercomparisons can reveal promising numerical solution approaches and increase confidence in the application of reactive transport codes. In this contribution, the results and performance of five current reactive transport codes are compared for the 1D and 2D subproblems of the so-called easy test case of the MoMaS benchmark (Carrayrou et al., Comput Geosci, 2009, this issue). This benchmark presents a simple fictitious reactive transport problem that highlights the main numerical difficulties encountered in real reactive transport problems. As a group, the codes include iterative and noniterative operator splitting and global implicit solution approaches. The 1D easy advective and 1D easy diffusive scenarios were solved using all codes, and, in general, there was a good agreement, with solution discrepancies limited to regions with rapid concentration changes. Computational demands were typically consistent with what was expected for the various solution approaches. The differences between solutions given by the three codes solving the 2D problem are more important. The very high computing effort required by the 2D problem illustrates the importance of parallel computations. The most important outcome of the benchmark exercise is that all codes are able to generate comparable results for problems of significant complexity and computational difficulty.

One of the challenging problems in mathematical geosciences is the determination of analytical solutions of nonlinear partial differential equations describing transport processes in porous media. We are interested in diffusive transport coupled with precipitation‐dissolution reactions. Several numerical computer codes that simulate such systems have been developed. Analytical solutions, if they exist, represent an important tool for verification of numerical solutions. We present a methodology for deriving such analytical solutions that are exact and explicit in space and time variables. They describe transport of several aqueous species coupled to precipitation and dissolution of a single mineral in one, two, and three dimensions. As an application, we consider explicit analytical solutions for systems containing one or two solute species that describe the evolution of solutes and solid concentrations as well as porosity. We use one of the proposed analytical solutions to test numerical solutions obtained from two conceptually different reactive transport codes. Both numerical implementations could be verified with the help of the analytical solutions and show good agreement in terms of spatial and temporal evolution of concentrations and porosities.

X-ray microtomography characterization of porosity, permeability and reactive surface changes during dissolution

High-resolution pore-scale simulation of dissolution in porous media

Abstract. We explore coupling to a configurable subsurface reactive transport code as a flexible and extensible approach to biogeochemistry in land surface models. A reaction network with the Community Land Model carbon–nitrogen (CLM-CN) decomposition, nitrification, denitrification, and plant uptake is used as an example. We implement the reactions in the open-source PFLOTRAN (massively parallel subsurface flow and reactive transport) code and couple it with the CLM. To make the rate formulae designed for use in explicit time stepping in CLMs compatible with the implicit time stepping used in PFLOTRAN, the Monod substrate rate-limiting function with a residual concentration is used to represent the limitation of nitrogen availability on plant uptake and immobilization. We demonstrate that CLM–PFLOTRAN predictions (without invoking PFLOTRAN transport) are consistent with CLM4.5 for Arctic, temperate, and tropical sites.Switching from explicit to implicit method increases rigor but introduces numerical challenges. Care needs to be taken to use scaling, clipping, or log transformation to avoid negative concentrations during the Newton iterations. With a tight relative update tolerance (STOL) to avoid false convergence, an accurate solution can be achieved with about 50 % more computing time than CLM in point mode site simulations using either the scaling or clipping methods. The log transformation method takes 60–100 % more computing time than CLM. The computing time increases slightly for clipping and scaling; it increases substantially for log transformation for half saturation decrease from 10−3 to 10−9 mol m−3, which normally results in decreasing nitrogen concentrations. The frequent occurrence of very low concentrations (e.g. below nanomolar) can increase the computing time for clipping or scaling by about 20 %, double for log transformation. Overall, the log transformation method is accurate and robust, and the clipping and scaling methods are efficient. When the reaction network is highly nonlinear or the half saturation or residual concentration is very low, the allowable time-step cuts may need to be increased for robustness for the log transformation method, or STOL may need to be tightened for the clipping and scaling methods to avoid false convergence.As some biogeochemical processes (e.g., methane and nitrous oxide reactions) involve very low half saturation and thresholds, this work provides insights for addressing nonphysical negativity issues and facilitates the representation of a mechanistic biogeochemical description in Earth system models to reduce climate prediction uncertainty.

In this study, CO2 transport in density-driven flows within an ideal model of a fractured porous medium, which contains a single or two intersecting fractures, is investigated numerically. The study employs a multi-scale modeling in which the flow in the matrix is modeled by Darcy's law, while the flow in the fracture is modeled by the Navier–Stokes equations. Our study shows that a horizontal fracture minimally impacts CO2 distribution, and depending on its length, slightly reduces dissolved CO2 during sequestration by 1.5%–2.5%. Vertical fractures play a crucial role in redirecting CO2 movement within the matrix, guiding it toward the fractures and altering its original pathway. Notably, the observed oscillations of CO2-rich water between the interfaces of the vertical fracture highlight the flow consistency with the pore scale. The domain-scale circulation induced by the vertical fracture leads to a rapid increase in flux and dissolved CO2 mass, but early convection shuts down. The results demonstrate that a longer vertical fracture leads to earlier shut down of convection and a potential decrease in storage of over 11%. The flow behaviors observed in inclined fractures are akin to those in vertical fractures, as they disrupt the fingerlike structure of CO2 around the fracture, form the circulation around the fracture, and are accompanied by vortices at the top. Additionally, intersecting fractures can lead to dynamic interactions between the fractures, with high-angle fractures dominating mixing flow. Different fracture angle combinations minimally affect dissolution mass.

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.

Reactive transport modeling has established itself as a key player to analyze sophisticated hydro-geochemical interactions that occur over spatiotemporal scales on par with subsurface applications. In this paper, we benchmarked a new reactive transport package—porousMedia4Foam at the continuum scale considering complex and well-established cases available in the reactive transport modeling community. porousMedia4Foam was born by the successful coupling of two open-source packages—OpenFOAM® and PHREEQC. The flow governing equations, transport of species and the evolution of porous media properties are handled by OpenFOAM®, whereas the geochemistry is exclusively handled by PHREEQC. We further demonstrated the capability of using porousMedia4Foam to investigate reactive transport processes considering unstructured meshes. As porousMedia4Foam is an open-source package with included benefits to account for reactive transport processes occurring at various scales—pore, hybrid, and Darcy scales, we believe that porousMedia4Foam opens a new dimension to analyze reactive transport physics for various intriguing subsurface applications.

Two types of fluids are encountered in the Onikobe geothermal reservoir (Japan): one is neutral and the other is acidic. It is hypothesized that acidic fluid might be upwelling along a fault zone from magma and that an impermeable barrier might be present between the acidic and neutral fluid zones. To test such a conceptual model and to study the geochemical behavior due to mixing of the two fluids, reactive geothermal transport simulations under both natural and production conditions were carried out using the code TOUGHREACT. Results indicate Mn‐rich smectite precipitates near the mixing front. Precipitation of sphalerite and galena occurs in a similar region as the Mn‐rich smectite. Precipitation of these minerals depends on pH and temperature. In addition, quartz, pyrite, and calcite precipitate in the shallow zone resulting in further development of caprock. The changes in porosity and permeability due to precipitation of Mn‐rich smectite are small compared with that of quartz, calcite, and pyrite. However, the smectite precipitation is likely to fill open fractures and to form an impermeable barrier between acidic and neutral fluid regions. The simulated mineral assemblage is generally consistent with observations in the Onikobe field. The numerical simulations described here provide useful insight into geochemical behavior and formation of impermeable barriers from fluid mixing. The method presented in this paper may be useful in fundamental analysis of hydrothermal systems and in the exploration of geothermal reservoirs, including chemical evolution, mineral alteration, mineral scaling, and changes in porosity and permeability.

A new numerical model for pore scale dissolution of calcite due to CO 2 saturated water flow in 3D realistic geometry: Principles and first results