Research topic: Reactive transport in porous media

Scientific context

Model description

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).

Publications and results

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

Objectives and projects

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.

Participants

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.

Work is now undertaken in collaboration with Andra, with the Ph-D of S. Sabit, 2010-2013.

Grants

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.

Simulations and Algorithms on Grids for Environment