Direct and split operator approaches with ELLAM for reactive transport equations

On the behavior of approaches to simulate reactive transport

A two-grid Eulerian–Lagrangian localized adjoint method to miscible displacement problems with dispersion term

SUMMARYThe Eulerian–Lagrangian Localized Adjoint Method (ELLAM) is well adapted for advection‐dominated transport problems. The method can use large time steps because the advection term is accurately approximated without any CFL restriction.A new formulation of ELLAM was developed by Younes et al. (Adv. Water Resour. 2006; 29:1056–1074) for unstructured triangular meshes. This formulation uses only strategic numerical integration points and thus requires a very limited number of integration points (usually 1 per element). To avoid numerical dispersion due to interpolation when several time steps are used, the method continuously tracks characteristics, and only changes due to dispersion are interpolated at each time step.In this work, we study the applicability and efficiency of this formulation for highly heterogeneous domains. The results show that (i) this formulation remains efficient even for highly heterogeneous domains, (ii) special care must be taken for the approximation of the dispersive term when large time steps are used, (iii) interpolation at intermediate times can be used to reduce memory requirement for problems with a large number of elements and small time steps, and (iv) the method can be used successfully for heterogeneous problems including injection and pumping wells, and remains much more efficient than the discontinuous Galerkin finite element method even when small time steps are required. Copyright © 2011 John Wiley & Sons, Ltd.

Reactive transport is a phenomenon resulting from the interaction and coupling of solute transport and chemical reactions. A new method to solve reactive transport known as Water Mixing Approach (WMA) was introduced by Soler-Sagarra et al. (2022). The idea is to interpret solute transport as water mixing and advection, where diffusion and dispersion are simulated as water exchange instead of Fickian solute flux. The WMA has the advantage of decoupling transport and chemistry. Transport computations are restricted to the evaluation of mixing ratios. This way, reactive transport computations are restricted to reactive mixing calculations, which can be performed separately for every target (node, cell, or particle, depending on the approach adopted to simulate transport). This facilitates parallelisation. However, the original work only considered the explicit case, which is conditionally stable and therefore requires artificial values of the dispersion coefficient to avoid numerical instabilities. We present a formulation of the WMA that is implicit both in transport, to ensure stability, and in chemical reactions to be able to simulate fast reactions. The implicit formulation requires lumping the reactive term. We test the validity of the approach by comparison with analytical solutions and the Direct Substitution Approach (DSA) method in a case with 2 adjacent mineral zones in equilibrium, and in a denitrification case with two redox reactions. We find that the proposed approach is extremely efficient and accurate for small dispersion cases.

SUMMARY The study investigates thermal-, hydraulic- and chemically coupled processes of enhanced geothermal systems (EGS). On the basis of the two existing numerical codes, the finite element program fracture and the geochemical module of chemtough, frachem was developed, to simulate coupled thermal–hydraulic–chemical (THC) processes, accounting for the Soultz specific conditions such as the high salinity of the reservoir fluid and the high temperatures. The finite element part calculates the thermal and hydraulic field and the geochemical module the chemical processes. According to the characteristics of the Soultz EGS reservoir, the geochemical module was modified. (i) The Debye–Huckel approach was replaced by the Pitzer formalism. (ii) New kinetic laws for calcite, dolomite, quartz and pyrite were implemented. (iii) The porosity–permeability relation was replaced by a new relation for fractured rock. (iv) The possibility of re-injecting the produced fluid was implemented. The sequential non-iterative approach (SNIA) was used to couple transport and reactions. Sensitivity analyses proved the proper functionality of frachem, but highlighted the sensitivity of the SNIA approach to time steps. To quantify the frachem results, a comparative simulation with the code shemat was conducted, which validated frachem. Coupled THC processes in a fractured zone in the Soultz reservoir at 3 500 m (T0= 165 °C), which occur as a result of the injection of fluid (Tinj= 65 °C) at one end of the zone and the production at the other end, were modelled for 2 yr. Calcite is the most reactive mineral and therefore the porosity and permeability evolution results from the calcite reactions: near the injection point, porosity and permeability increase and near the production well they decrease. After 2 yr, the system seems to be very close to steady-state. Therefore, mineral dissolution and precipitation during the circulation of the fluid in the reservoir do not represent a limiting factor on the EGS reservoir (at 3 500 m depth) performance at the Soultz site. Finally, the numerical transmissivity was compared to the transmissivity of the 1997 circulation test. The fact that the transmissivity decreases during the circulation test, when thermomechanical effects are factored out, points to geochemical processes in the reservoir such as the precipitation of calcite. These findings highlight the importance of THC coupled EGS reservoir models. The integration of geochemical considerations is therefore indispensable for integrated simulations of EGS systems and predictions of its performance.

A comparison of two approaches for reactive transport modelling

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.

The present study is concerned with the modeling of conservative as well as non-conservative solute transport in ground water. The model is based on an operator split approach which uses an Eulerian frame work with finite volume method for the advective transport and fully implicit central difference method for the dispersive transport. This formulation helps in accurately simulating both highly advective and dispersive transport cases with less restriction on the grid size and time step. The numerical solution is compared with exact as well as the approximate analytical solution and the maximum error as percentage of peak concentration is presented for different values of decay constant at varying distances. It is found that finite volume method provides accurate solutions of both conservative and non-conservative solutes for both advection and dispersion dominated situations.

Reactive transport equations may become cumbersome to solve when a large number of species are coupled through fast (equilibrium) and slow (kinetic) reactions. Solution is further encumbered when both mobile (solutes) and immobile (minerals) species are involved. In these cases, sequential iteration approaches (SIA) may become extremely slow to converge. Direct substitution approaches (DSA), which solve all transport equations together, may become extremely large. Here we propose a formulation for optimally decoupling the reactive transport equations. The procedure is described sequentially using a paradigm system. We start by a tank paradigm in which all species are mobile and undergo equilibrium reactions. In this case, all components are fully decoupled. If some of the reactions are kinetic (canal paradigm), then we can still decouple components corresponding to equilibrium reactions. The same can be said regarding a system in which immobile species only react kinetically (river). The number of components can be reduced in the case that both types of reactions and species are present (aquifer). In short, the number of coupled equations that need to be solved simultaneously is, at most, equal to the number of kinetic reactions. This benefits both SIA and DSA solution methods. SIA should improve convergence because most components are linear and effectively decoupled, thus reducing nonlinearity to kinetic terms only in transport equations for kinetic components. DSA systems become reduced as the number of components that need to be solved together is, at most, equal to the number of independent kinetic reactions.

Application of a reactive transport model for interpreting non-conservative tracer experiments: The Rio Maior case-study

Bioremediation is a promising strategy for sustainable management and treatment of field sites with soil and groundwater contamination. Understanding the fate of contaminants is key for identifying optimal conditions and prediction of biodegradation in the subsurface. We develop the open-source Python package mibitrans for hydrogeological reactive transport modelling for simple use of field site modelling, transport prediction and management optimization. Mibitrans is developed in the frame of the EU-funded MiBiRem project.Mibitrans is based on physics-based (semi-)analytical solutions for multidimensional solute transport in groundwater. Biodegradation can be modelled in various ways, including simple first-order decay and chemical reactions based on electron acceptor concentrations. The mibitrans package is modular, fully tested, well documented and well-structured to allow for easy adaptation and use in various field situations. We will present the use with example field site data from ongoing bioremediation projects.

Flow and solute transport in groundwater are primarily controlled by the hydraulic conductivity (K) field. Variations in the K field greatly influence subsurface flow rates and solute migration and dispersion. Due to the importance of aquifer heterogeneity, several approaches have been proposed in the literature for modelling K field spatial distribution based on borehole data. In this work, we evaluate two methods to generate spatial realisations of the K field of a heterogeneous aquifer. Results are analysed in terms of the connectivity of high-K cells, and solute transport results are compared and discussed. The methodology consists of three main steps. First, borehole PSD data are used to characterise hydrofacies and generate 3D stochastic realisations of those hydrofacies using the Markov-Chain/Transition Probability approach (MC/TP). Resultant realisations are used to fill a 3D grid. Then, two methods are used to generate hydraulic conductivity fields by assigning a K value to each cell on those grids: (1) using the geometric mean of each hydrofacies (GMEAN) and (2) using a value sampled from the KDE probability functions of each hydrofacies (KDE). Two different porosity scenarios are considered (high porosity, HP; and low porosity, LP). Finally, solute transport estimates were computed using MODFLOW and MT3DMS.Connectivity analysis of the stochastic realisations shows a higher degree of connectivity of high-K cells on the GMEAN mode than in the KDE. The latter leads to overall higher average K values, but also to slower flow regions with lower K values than the GMEAN.Solute transport runs result in slower travel times and lower peak concentrations on the KDE realisations than in the GMEAN case. Breakthrough curves at different observation points show that, when concentrations fall after peaking (for a time-limited input pulse), both GMEAN and KDE curves tend to converge. The field-scale longitudinal dispersivity implied from the ensemble probability plumes is similar between the two K field realisation approaches modelled, but the vertical dispersivity is higher in the KDE realisations. The high porosity scenario shows higher dispersivities and considerably longer travel times.Results show that, for the same porosity scenario, contaminant plumes behave similarly for the KDE and GMEAN approaches at longer times and distances from the source. This suggests that, on the studied site, the small-scale heterogeneity has a reduced effect on the on the long term, field scale macrodispersivity and solute transport behaviour.

Iteration Scheme for Implicit Calculations of Kinetic and Equilibrium Chemical Reactions in Fluid Dynamics

A finite element formulation is developed for the simulation of groundwater flow and solute transport in multilayer systems of several aquifers and aquitards. This formulation is general, flexibile, and capable of taking full advantage of the nature of flow in such multilayer systems. A fully three‐dimensional spatial representation can be performed for certain aquifers or for the entire flow system if needed. Those parts of the system that require three‐dimensional spatial discretization are handled effectively by combining two‐dimensional basis functions in the x‐y space and one‐dimensional basis functions in the z space. Furthermore, the formulation has a desirable option to perform one‐dimensional representation of flow and transport in aquitards and areal representation of flow and transport in certain aquifers that do not require the fully three‐dimensional discretization. When the one‐dimensional representation of an aquitard is used, coupling of adjacent aquifer and aquitard layers is handled using a convolution integral approach. A general solution strategy is also developed to allow systematic time stepping and cost‐effective matrix handling schemes. For those parts of the system that require three‐dimensional discretization, an algorithm referred to as the alternate sublayer and line sweep procedure is presented for decomposing three‐dimensional matrix equations. This algorithm can accomodate several thousand nodal unknowns without requiring excessive core storage and CPU time. Four examples of transient flow and transport problems are provided. These examples show verification of numerical results and demonstrate that the present finite element models are far more cost‐effective than earlier three‐dimensional models that are based on the conventional Galerkin approach and direct matrix solution algorithms.

This work presents a new reactive transport framework that combines a powerful geochemistry engine with advanced numerical methods for flow and transport in subsurface fractured porous media. Specifically, the PhreeqcRM interface (developed by the USGS) is used to take advantage of a large library of equilibrium and kinetic aqueous and fluid-rock reactions, which has been validated by numerous experiments and benchmark studies. Fluid flow is modeled by the Mixed Hybrid Finite Element (FE) method, which provides smooth velocity fields even in highly heterogenous formations with discrete fractures. A multilinear Discontinuous Galerkin FE method is used to solve the multicomponent transport problem. This method is locally mass conserving and its second order convergence significantly reduces numerical dispersion. In terms of thermodynamics, the aqueous phase is considered as a compressible fluid and its properties are derived from a Cubic Plus Association (CPA) equation of state. The new simulator is validated against several benchmark problems (involving, e.g., Fickian and Nernst-Planck diffusion, isotope fractionation, advection-dispersion transport, and rock-fluid reactions) before demonstrating the expanded capabilities offered by the underlying FE foundation, such as high computational efficiency, parallelizability, low numerical dispersion, unstructured 3D gridding, and discrete fraction modeling.