The transport equations are mass conservation with advection and dispersion. The velocity field comes from previous steady-state flow computations. The dispersion tensor is assumed independent of chemical species. Chemistry equations can be at equilibrium or kinetic. They include aqueous reactions, sorption, ion exchanges, precipitation or dissolution. Secondary species are related to components by mass action laws and total concentrations are expressed by means of mass conservation. Transport equations can be written for the total concentrations of mobile components. Equations are a set of coupled nonlinear Partial Differential Algebraic Equations (PDAE), where algebraic chemistry equations are written at each point of the domain.

The transport equations are spatially discretized by a finite element method with an unstructured mesh ; advection is discretized by an upwind scheme and dispersion is discretized by a centered scheme. In the framework of methods of lines, it allows to use any ODE solver after spatial discretization.

Three types of numerical couplings are compared :

- a so-called Sequential Non Iterative Approach; time
discretization is done with the explicit Euler scheme. The transport
equations gives the total concentrations, which are used to solve the
algebraic chemistry equations everywhere in the domain.

- a so-called Sequential Iterative Approach; time discretization is done with the implicit Euler scheme. The nonlinear problem is solved at each timestep by a block nonlinear relaxation, solving first the transport equations for mobile species with immobile species as a reaction term, then the chemistry equations with an updated total of mobile and immobile species.
- a so-called global approach; it follows the method of lines
with a set of discrete DAE to solve. We use the solver Sundials, with
adaptive timesteps and orders. Moreover, the nonlinear system is solved
by a modified Newton method at each timestep.

In a simplified model, the number of fixed species is known in
advance and does not change during the simulation. In a more general
model, this number can vary and the model includes a nonlinear
complementarity formulation, solved with a semismooth Newton method.

The simplified model with only equilibrium reactions and the different couplings are implemented in the software GRT3D.

The global approach has been successfully applied to the Momas benchmark on reactive transport (easy test case, 1D and 2D). It has also been applied to Alliances test cases (1D and 2D).

- Erhel, J. and Sabit, S. and de Dieuleveult, C.

B.H.V. Topping and P. Ivanyi (ed)

Developments in Parallel, Distributed, Grid and Cloud Computing for Engineering

Solving Partial Differential Algebraic Equations and Reactive Transport Models

Saxe-Coburg Publications, 2013 - J. Carrayrou and J. Hoffmann and P. Knabner
and S. Kräutle and C. de Dieuleveult and J. Erhel and J. Van der Lee
and V. Lagneau and K.U. Mayer and K.T.B. MacQuarrie

Comparison of numerical methods for simulating strongly non-linear and heterogeneous reactive transport problems – the MoMaS benchmark case

*Computational Geosciences,***2010,**Volume 14, 483-502 - C. de Dieuleveult and J.
Erhel

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

*Computational Geosciences,***2010,**Volume 14,451-464 -
C. de Dieuleveult, J. Erhel and M. Kern

A global strategy for solving reactive transport equations

*Journal of Computational Physics,***2009,**Volume 228, 6395-6410 -
C. de Dieuleveult and J. Erhel

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

*Computational Geosciences,***2009** - de Dieuleveult, C. & Erhel, J.

A numerical model for coupling chemistry and transport

*International Conference on SCIentific Computation And Differential Equations, SciCADE 2007 (contributed talk),***2007** - Erhel, J. & de Dieuleveult, C.

Nonlinear methods for reactive transport simulations

*Int. conf. on Approximation and Iterative Methods (invited talk),***2006**

The main objective is now to reduce this CPU time. A first approach is to reduce the size of the linear system by a subsitution technique. A second approach is to deal with the tolerance and convergence parameters of the DAE solver. Another objective is to implement the semismooth Newton method of the general model and to implement a model with kinetic reactions. Also, a posteriori error estimations will allow to refine adaptively the mesh and the timestep.

- transport in porous media
- flow in porous media
- software GRT3D

This topic started in the Sage team (participant J. Erhel) with the Hydrogrid
project (2002-2005), in collaboration with J. Carrayrou, from IMFS at
Strasbourg and M. Kern, from the team Estime at
INRIA-Rocquencourt.

The research continued in the Sage team (participant J. Erhel), with the
Ph-D of C. de Dieuleveult, 2005-2008, and still in collaboration with J. Carrayrou, M. Kern and with CEA, in 2004-2012.

This work was supported by the ACI Grid with the project Hydrogrid in 2002-2005.

It was supported by a grant from ANDRA in 2005-2008 and by projects of the GdR Momas in 2004-2012.

C. de Dieuleveult was hired by ANDRA during her Ph-D. She has now a research position at Ecole des Mines de Paris.

The Ph-D of S. Sabit is funded by a grant from ANDRA, 2010-2013.

The work is supported by ANR, with the project H2MNO4, 2013-2017.