A global method for coupling transport with chemistry in heterogeneous porous media

Modeling reactive transport in porous media, using a local chemical equilibrium assumption, leads to a system of advection–diffusion PDEs coupled with algebraic equations. When solving this coupled system, the algebraic equations have to be solved at each grid point for each chemical species and at each time step. This leads to a coupled non-linear system. In this paper, a global solution approach that enables to keep the software codes for transport and chemistry distinct is proposed. The method applies the Newton–Krylov framework to the formulation for reactive transport used in operator splitting. The method is formulated in terms of total mobile and total fixed concentrations and uses the chemical solver as a black box, as it only requires that one be able to solve chemical equilibrium problems (and compute derivatives) without having to know the solution method. An additional advantage of the Newton–Krylov method is that the Jacobian is only needed as an operator in a Jacobian matrix times vector product. The proposed method is tested on the MoMaS reactive transport benchmark.     

Finite element methods for variable density flow and solute transport

Saltwater intrusion into coastal freshwater aquifers is an ongoing problem that will continue to impact coastal freshwater resources as coastal populations increase. To effectively model saltwater intrusion, the impacts of increased salt content on fluid density must be accounted for to properly model saltwater/freshwater transition zones and sharp interfaces. We present a model for variable density fluid flow and solute transport where a conforming finite element method discretization with a locally conservative velocity post-processing method is used for the flow model and the transport equation is discretized using a variational multiscale stabilized conforming finite element method. This formulation provides a consistent velocity and performs well even in advection-dominated problems that can occur in saltwater intrusion modeling. The physical model is presented as well as the formulation of the numerical model and solution methods. The model is tested against several 2-D and 3-D numerical and experimental benchmark problems, and the results are presented to verify the code.     

HYTEC results of the MoMas reactive transport benchmark

A specific benchmark has been developed by the French research group MoMas in order to improve numerical solution methods applied by reactive transport models, i.e., codes that couple hydrodynamic flow and mass transport in porous media with geochemical reactions. The HYTEC model has been applied to this benchmark exercise, and this paper summarizes some of the principal results. HYTEC is a general-purpose code, applied by industrials and research groups to a wide variety of domains, including soil pollution, nuclear waste storage, cement degradation, water purification systems, storage of CO₂, and valorization of stabilized wastes. The code has been applied to the benchmark test-cases without any specific modification. Apart from the benchmark imposed output, additional information is provided to highlight the behavior of HYTEC specifically and the simulation results in particular.     

Reactive transport benchmark of MoMaS

We present here the definition of the reactive transport benchmark of Groupement Mathematical Modeling and Numerical Simulation for Nuclear Waste Management Problems. The aim of this benchmark is to propose a challenging test for numerical methods used for reactive transport modeling in porous media. In order to focus on numerical methods, the problem presented here is of quite a small size, both from a hydrodynamical and from a geochemical point of view. Though the chemical coefficients used in this benchmark are not taken from a real chemical system, they are realistic, and the test case is quite challenging.     

Solution of the MoMaS reactive transport benchmark with MIN3P—model formulation and simulation results

This paper summarizes the governing equations as implemented in the MIN3P multicomponent flow and reactive transport code (Mayer et al., Water Resour Res 38:1174, 2002) and introduces the equations in discretized form. Linearization and solution methods are presented including adaptive time stepping and update modification schemes. Code-specific details for the implementation of the GdR MoMaS benchmark simulations (Carrayrou et al., Comput Geosci, 2009) are presented. The standard version of the MIN3P code was used to solve the Easy, Medium, and Hard Test Cases, in one and two spatial dimensions, for both advection- and diffusion-dominated conditions. An analysis of the sensitivity of the solution in relation to spatial and temporal discretization parameters is provided for the Easy Test Case, selected results are presented for the Medium and Hard Test Cases, and the performance of the code as a function of discretization parameters is evaluated for all test cases.     

Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems—the MoMaS benchmark case

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.     

Modeling reactive transport in non-ideal aqueous–solid solution system

The numerical simulation of reactive mass transport processes in complex geochemical environments is an important tool for the performance assessment of future waste repositories. A new combination of the multi-component mass transport code GeoSys/RockFlow and the Gibbs Energy Minimization (GEM) equilibrium solver GEM-Selektor is used to calculate the accurate equilibrium of multiple non-ideal solid solutions which are important for the immobilization of radionuclides such as Ra. The coupled code is verified by a widely used benchmark of dissolution–precipitation in a calcite–dolomite system. A more complex application shown in this paper is the transport of Ra in the near-field of a nuclear waste repository. Depending on the initial inventories of Sr, Ba and sulfate, non-ideal sulfate and carbonate solid solutions can fix mobile Ra cations. Due to the complex geochemical interactions, the reactive transport simulations can describe the migration of Ra in a much more realistic way than using the traditional linear K_D approach only.     

A general reduction scheme for reactive transport in porous media

We present a method to transform the governing equations of multispecies reactive transport in porous media. The reformulation leads to a smaller problem size by decoupling of equations and by elimination of unknowns, which increases the efficiency of numerical simulations. The reformulation presented here is a generalization of earlier works. In fact, a whole class of transformations is now presented. This class is parametrized by the choice of certain transformation matrices. For specific choices, some known formulations of reactive transport can be retrieved. Hence, the software based on the presented transformation can be used to obtain efficiency comparisons of different solution approaches. For our efficiency tests, we use the MoMaS benchmark problem on reactive transport.     

A sequential partly iterative approach for multicomponent reactive transport with CORE2D

CORE^2D V4 is a finite element code for modeling partly or fully saturated water flow, heat transport, and multicomponent reactive solute transport under both local chemical equilibrium and kinetic conditions. It can handle coupled microbial processes and geochemical reactions such as acid–base, aqueous complexation, redox, mineral dissolution/precipitation, gas dissolution/exsolution, ion exchange, sorption via linear and nonlinear isotherms, and sorption via surface complexation. Hydraulic parameters may change due to mineral precipitation/dissolution reactions. Coupled transport and chemical equations are solved by using sequential iterative approaches. A sequential partly iterative approach (SPIA) is presented which improves the accuracy of the traditional sequential non-iterative approach (SNIA) and is more efficient than the general sequential iterative approach (SIA). While SNIA leads to a substantial saving of computing time, it introduces numerical errors which are especially large for cation exchange reactions. SPIA improves the efficiency of SIA because the iteration between transport and chemical equations is only performed in nodes with a large mass transfer between solid and liquid phases. The efficiency and accuracy of SPIA are compared to those of SIA and SNIA using synthetic examples and a case study of reactive transport through the Llobregat Delta aquitard in Spain. SPIA is found to be as accurate as SIA while requiring significantly less CPU time. In addition, SPIA is much more accurate than SNIA with only a minor increase in computing time. A further enhancement of the efficiency of SPIA is achieved by improving the efficiency of the Newton–Raphson method used for solving chemical equations. Such an improvement is obtained by working with increments of log concentrations and ignoring the terms of the Jacobian matrix containing derivatives of activity coefficients. A proof is given for the symmetry and non-singularity of the Jacobian matrix. Numerical analyses performed with synthetic examples confirm that these modifications improve the efficiency and convergence of the iterative algorithm. Changbing Yang is now at The University of Texas at Austin, USA.     

An efficient finite–discrete element method for quasi‐static nonlinear soil–structure interaction problems

An efficient finite–discrete element method applicable for the analysis of quasi‐static nonlinear soil–structure interaction problems involving large deformations in three‐dimensional space was presented in this paper. The present method differs from previous approaches in that the use of very fine mesh and small time steps was not needed to stabilize the calculation. The domain involving the large displacement was modeled using discrete elements, whereas the rest of the domain was modeled using finite elements. Forces acting on the discrete and finite elements were related by introducing interface elements at the boundary of the two domains. To improve the stability of the developed method, we used explicit time integration with different damping schemes applied to each domain to relax the system and to reach stability condition. With appropriate damping schemes, a relatively coarse finite element mesh can be used, resulting in significant savings in the computation time. The proposed algorithm was validated using three different benchmark problems, and the numerical results were compared with existing analytical and numerical solutions. The algorithm performance in solving practical soil–structure interaction problems was also investigated by simulating a large‐scale soft ground tunneling problem involving soil loss near an existing lining. Copyright © 2011 John Wiley & Sons, Ltd.     

Hydrogeochemical modeling of organo-metallic colloids in the Nsimi experimental watershed,South Cameroon

This paper presents a hydrogeochemical modeling code HYDROS, which combines the multi-component transport model with equilibrium speciation module MINTEQA2. The processes of adsorption, aqueous speciation and mineral precipitation/dissolution are represented in the model. The numerical model uses a sequential iterative approach for solving the solute transport and the equilibrium geochemistry modules. Further the transport part is solved using an operator split approach wherein a finite volume method is used for solving the advective equations while a classical finite difference method is employed for solving the dispersive equations. The model performance is evaluated by comparing it with MINTOX for a literature problem. HYDROS is then applied to the case study of the transfer of transition metals with organic colloids in the swamp groundwater system of the experimental Nsimi watershed, representative of the humid tropical ecosystem of the South Cameroon Plateau. Field observations at the site swamp system suggest that the carbon is mainly transferred as organic colloids (i.e., dissolved organic carbon) produced by the slow biodegradation of the swamp organic matter. Using HYDROS, the behaviour of Al(III) and Fe(III) elements in the base flow system is simulated during inter rain events of a short rainy season (May–June 1996). The elemental time-series for Al, Fe, Cl, pH compare well with the simulation results. The colloids are found to have a strong impact on the mobilization and transfer of Al(III) and Fe(III), which are considered to have low mobility in weathering environment.     

A fast flux tube-based method for solute-transport simulation

A new method to calculate the transport of dissolved species in aquifers is presented. This approach is an extension of the stream tubes which are used for flow computation. The flux tubes defined here are conservative for solutes, but not for water mass. The flux tubes are first defined in a general domain and then calculated in a two-dimensional uniform flow field. The tubes?? computation is based on a parametric solution. The method is extended further in order to deal with heterogeneous media. A particle-tracking algorithm is used where the deviation of the flux-tube boundaries due to dispersion is accounted for. The approximate solution obtained by this approach is compared to classical numerical solutions given by a finite difference code (RT3D) and a finite element code (FEFLOW). This comparison was performed for several test cases with increasing complexity. The differences between the flux-tube approach and the other methods always remain small, even regarding mass conservation. The major advantage of the flux-tube approach is the ability to reach a solution quickly, as the method is hundreds to thousands of times faster than classical finite difference or finite element models.     

Modelling Microbial Degradation Coupled to Reactive Transport in Groundwater: A Benchmark Analysis

Microbes are ubiquitous in groundwater systems, and they play an important role in the redox state of groundwater and especially on the fate of organic contaminants. In this context, numerical simulations that couple microbial processes to reactive transport models are becoming more popular. In the present work, we revisit the mathematical ground of microbial redox reactions and perform a benchmark analysis of the simulation of aerobic benzene degradation in a shallow and oxidizing aquifer. Numerical results indicate that the two codes tested (one using the finite elements approach and the other using the finite differences approach) lead to very similar results. In addition, the coupling of heterogeneous geochemical reactions to the benchmarked example problem provides a solid basis for the understanding of the redox reactions and the changes on the carbon system triggered by the aerobic degradation of benzene.     

Determination of Zn partition coefficient and simulation of reactive transport from landfills in Mar Del Plata, Argentina

The distribution coefficient (Kd) expresses the relationship between the concentration of an element, which is adsorbed in the solid surface and its remaining concentration in the solution. The Kd is a very important factor in reactive transport, representing the source/sink term, and explaining the difference between the velocity of transport of non-conservative elements (Kd>0) and water flow velocity. In this paper, the Kd value for Zn element in loess like sediments forming the Pampeano aquifer is determined and this value is used in the modeling of reactive transport from the landfill of the city of Mar del Plata (Argentina). The determination of Kd value was done by means of batch experiments. The results obtained showed good agreement with Freundlich isotherm, with a value of K_F=300.95 ml g⁻¹ and a super index value b of 0.3961. These values were applied to reactive transport modeling using Visual Modflow code. The Zn plume obtained showed the low mobility of the element in the oxidizing conditions of the environment.     

Mixed finite element discretization and Newton iteration for a reactive contaminant transport model with nonequilibrium sorption: convergence analysis and error estimates

We present a numerical scheme for reactive contaminant transport with nonequilibrium sorption in porous media. The mass conservative scheme is based on Euler implicit, mixed finite elements, and Newton method. We consider the case of a Freundlich-type sorption. In this case, the sorption isotherm is not Lipschitz but just Hölder continuous. To deal with this, we perform a regularization step. The convergence of the scheme is analyzed. An explicit order of convergence depending only on the regularization parameter, the time step, and the mesh size is derived. We give also a sufficient condition for the quadratic convergence of the Newton method. Finally, relevant numerical results are presented.     

无压渗流问题分析的多节点有限元方法

采用基于平均值插值的多节点有限元方法分析有自由面渗流问题。在自由面附近采用多节点单元逼近自由面,利用平均值插值建立多节点单元的形函数,在远离自由面的区域采用四边形/三角形单元剖分。给定一个初始自由面位置,通过对渗流控制方程的多节点有限元求解,根据自由面上节点水头值判断自由面节点的调整方向和大小,最终迭代求出自由面的位置。土坝渗流问题的数值分析表明了所提方法的有效性和足够的计算精度。     

Two-dimensional numerical finite volume modeling of processes controlling distribution and natural attenuation of BTX in the saturated zone of a simulated semi-confined aquifer

This paper presents a two-dimensional finite volume model to predict multi-species reactive transport processes in the saturated zone of a simulated semi-confined aquifer. A multipurpose commercial software called PHOENICS was used to solve model equations numerically. Capability of the present model was first confirmed using experimental data and the results obtained by a published one-dimensional finite element reactive transport model by other researchers taking different scenarios into consideration. The model was then expanded to a two-dimensional case to simulate reactive transport of BTX compounds with discontinuous source in the saturated zone of the groundwater flow system. In addition to the physical transport processes, the two-dimensional model also incorporates linear and nonlinear adsorption isotherms, first order and Monod kinetics. The two-dimensional model considers both static and dynamics modes into account. The results show that considering chemical reactions during reactive transport of contaminants could successfully predict the contaminated zone. The results of such studies can be used for monitoring of contaminated areas, designing methods to control pollution transport, and minimize its harmful effects on aquifer systems.     

求解非稳定对流占优水质模型的非标准Galerkin有限元法

在求解非稳定地下水溶质运移模型时,若对流项占优,则模型表现出双曲方程的特性。针对这种特性,采用非标准Galerkin有限元方法进行求解是解决这类问题的有效途径。分别采用Wavelet-Galer-kin有限元方法、迎风有限元方法和特征有限元方法对强对流溶质运移模型进行了求解,并将其结果与标准Galerkin有限元和解析解进行对比。结果表明:标准Galerkin有限元方法会产生强烈的数值振荡;Wavelet-Galerkin有限元方法的时空定位效果好;迎风有限元方法能够有效降低数值振荡现象,但迎风因子对解的影响较大,而且会带来时间延迟;特征有限元方法能够提高解的精度,故可以认为特征有限元方法是求解强对流地下水溶质运移模型的首选方法。     

A global approach to reactive transport: application to the MoMas benchmark

Water resource management involves numerical simulations in order to study contamination of groundwater by chemical species. Not only do the aqueous components move due to physical advection and dispersion processes, but they also react together and with fixed components. Therefore, the mass balance couples transport and chemistry, and reactive transport models are partial differential equations coupled with nonlinear algebraic equations. In this paper, we present a global method based on the method of lines and differential algebraic system (DAE) solvers. At each time step, nonlinear systems are solved by a Newton-LU method. We use this method to carry out numerical simulations for the reactive transport benchmark proposed by the MoMas research group. Although we study only 1D computations with a specific geochemical system, several difficulties arise. Numerical experiments show that our method can solve quite difficult problems, get accurate results and capture sharp fronts.