Laila Amir
email: l.amir@uca.ma

Michel Kern
email: michel.kern@inria.fr

PRECONDITIONING A COUPLED MODEL FOR REACTIVE TRANSPORT IN POROUS MEDIA *

Keywords: Reactive transport, porous media, precondtioning, non-linear systems AMS subject classifications. 76S05, 65H10, 65F10, 65N22

We study numerical methods for solving reactive transport problems in porous media that allow a separation of transport and chemistry at the software level, while keeping a tight numerical coupling between both subsystems. We give a formulation that eliminates the local chemical concentrations and keeps the total concentrations as unknowns, then recall how each individual subsystem can be solved. The coupled system is solved by a Newton-Krylov method. The block structure of the model is exploited both at the nonlinear level, by eliminating some unknowns, and at the linear level by using block Gauss-Seidel or block Jacobi preconditioning. The methods are applied to a 1D case of the MoMaS benchmark.

1. Introduction. Reactive transport in porous media studies the coupling of reacting chemical species with mass transport in the subsurface, see [5, sec. 7.9] for an introduction, or [START_REF] Appelo | Geochemistry, Groundwater and Pollution[END_REF][START_REF] Bethke | Geochemical and Biogeochemical Reaction Modeling[END_REF] for more comprehensive references. It plays an important role in several applications when modeling subsurface flow [START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF][START_REF] Steefel | Reactive transport modeling: An essential tool and a new research approach for the earth sciences[END_REF][START_REF]Groundwater Reactive Transport Models[END_REF] :

• chemical trapping is one of the mechanism by which the safety of geological sequestration of carbon dioxide in deep saline aquifers can be ensured [START_REF] Audigane | Two-Dimensional reactive transport modeling of CO2 injection in a saline aquifer at the Sleipner site, North sea[END_REF][START_REF] Fan | A fully-coupled flow-reactive-transport formulation based on element conservation, with application to CO 2 storage simulations[END_REF][START_REF] Saaltink | A method for incorporating equilibrium chemical reactions into multiphase flow models for CO2 storage[END_REF][START_REF] Trenty | A sequential splitting strategy for CO2 storage modelling[END_REF] • nuclear waste storage is based on a multiple barrier concept so as to delay the arrival of radionuclides in the bio-sphere. The concrete barrier that seals the repository may be attacked by oxidized compounds, and again its safety needs to be assessed [START_REF] De Windt | Coupled modeling of cement/claystone interactions and radionuclides migration[END_REF][START_REF] Marty | Influence of reaction kinetics and mesh refinement on the numerical modelling of concrete/clay interactions[END_REF][START_REF] Samper | Testing Kd models of Cs+ in the near field of a HLW repository in granite with a reactive transport model[END_REF]; • bio-geochemistry involves interactions with organic chemicals, and is important for studies of soil pollution, and also for implementing bio-remediation policies [START_REF] Chilakapati | An analysis of complex reaction networks in groundwater modeling[END_REF][START_REF] Marzal | Modeling multiple reactive solute transport with adsorption under equilibrium and nonequilibrium conditions[END_REF].

The numerical simulation of reactive transport has been the topic of numerous work. The survey by Yeh and Tripathi [START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF] has been very influential in establishing a mathematical formalism for setting up models, and also for establishing the "operator splitting" approach (see below) as a standard. More recent surveys, detailing several widely used computer codes and their applications, can be found in the book [START_REF]Groundwater Reactive Transport Models[END_REF] and the survey article [START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF].

One traditionally distinguishes two main approaches for solving reactive transport problems: Sequential Iterative Approach (SIA) in this family of methods (also known as operator splitting), transport and chemistry are solved alternatively [START_REF] Carrayrou | Operator-splitting procedures for reactive transport and comparison of mass balance errors[END_REF][START_REF] Lagneau | HYTEC results of the MoMas reactive transport benchmark[END_REF][START_REF] Mayer | Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions[END_REF][START_REF] Samper | A sequential partly iterative approach for multicomponent reactive transport with CORE2D, Computat[END_REF][START_REF] Yeh | A model for simulating transport of reactive multispecies components: Model development and demonstration[END_REF]. This has the advantage that a code for reactive transport can be built from pre-existing transport and chemistry codes [START_REF] Heredia | TReacLab: an object oriented implementation of non intrusive operator splitting methods to couple independent transport and geochemical software[END_REF], and that no global system of equations needs to be solved. On the other hand, the splitting between chemistry and transport may restrict the time step in order to ensure convergence of the method, and the splitting may also introduce mass errors that need to be controlled [START_REF] Valocchi | Accuracy of operator splitting for advection-dispersionreaction problems[END_REF]. As the above references show, the method can be quite successful if it is implemented carefully. Globally implicit Approach (GIA) Here one solves the fully coupled system involving transport and chemistry is solved in one shot (see for instance [START_REF] Fahs | On the efficiency of the direct substitution approach for reactive transport problems in porous media[END_REF][START_REF] Saaltink | On the behavior of approaches to simulate reactive transport[END_REF] and the references below). The balance is the opposite of what it was for SIA: the method does not introduce spurious mass errors, and can converge with large time steps, but it leads to a large and difficult to solve system of nonlinear equations coupling all the chemical species at all grid points. The system can be solved by substituting the mass action laws into the conservation equations, a variant known as Direct Substitution approach (DSA) as in [START_REF] Fan | A fully-coupled flow-reactive-transport formulation based on element conservation, with application to CO 2 storage simulations[END_REF][START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physics-based preconditioning to biogeochemical transport[END_REF]. Recently, methods that solve the coupled problem without direct substitution have been introduced. In [START_REF] Hoffmann | A parallel global-implicit 2-D solver for reactive transport problems in porous media based on a reduction scheme and its application to the momas benchmark problem[END_REF][START_REF] Kräutle | A new numerical reduction scheme for fully coupled multicomponent transport-reaction problems in porous media[END_REF][START_REF]A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions[END_REF], local chemical systems are solved, and the solution becomes an additional term in a transformed form of the transport. In [START_REF] De Dieuleveult | A global approach to reactive transport: application to the MoMas benchmark[END_REF][START_REF] Erhel | Analysis of a global reactive transport model and results for the momas benchmark[END_REF][START_REF] Erhel | Solving partial differential algebraic equations and reactive transport models[END_REF], the overall coupled system is solved as a Differential Algebraic System. In addition to the surveys mentioned above, the methods have been compared in several studies or benchmarks [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF][START_REF] Marty | Benchmarks for multicomponent reactive transport across a cement/clay interface[END_REF][START_REF] Saaltink | A comparison of two approaches for reactive transport modelling[END_REF]. It is fair to say that no method or code emerges as a clear winner for all situations.

We finish this short (and far from exhaustive) review of the literature by noting that most of the work cited deal with one-phase flow and transport. The methods have recently been extended to the case of two-phase flow [START_REF] Ahusborde | Numerical simulation of two-phase multicomponent flow with reactive transport in porous media: application to geological sequestration of CO2[END_REF][START_REF] Fan | A fully-coupled flow-reactive-transport formulation based on element conservation, with application to CO 2 storage simulations[END_REF][START_REF] Saaltink | A method for incorporating equilibrium chemical reactions into multiphase flow models for CO2 storage[END_REF][START_REF] Sin | Integrating a compressible multicomponent two-phase flow into an existing reactive transport simulator[END_REF] .

In a previous paper [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF], the authors introduced a method that belongs to the GIA family without sacrificing the ease of implementation of the SIA methods, due to the separation of software modules for chemistry and transport. The method solves the nonlinear coupled system by a Jacobian-free Newton-Krylov method.

This paper presents several improvements to [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF]:

• a systematic study of preconditioning methods for the linearized coupled problem, and their relationship to elimination methods. This results in a method where the mobile concentrations are eliminated, and one solves for the total fixed concentrations. We give both heuristic arguments and experimental evidence that the convergence rate for both the Newton and GMRES iterations is bounded independently of the mesh size; • a formulation of the chemical equilibrium problem that does not need the a priori separation of the chemical species between primary and secondary components, but that keeps the distinction between mobile and immobile species. The method originally introduced in [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF] approximated the Jacobian matrix by vector product by finite difference, so that only solvers for transport and chemistry were needed, and they could be applied as black boxes. In this work, we slightly "open the boxes", as we propose to compute exactly the Jacobian matrix by vector products, so that access to the Jacobians of the transport and chemical solvers are needed. We feel the added accuracy and reduced cost (for the chemical part, the inverting the jacobian is much less expensive than solving the whole system) more than make up for the additional requirement on the software.

Since the convergence of Krylov solvers can be slow, we devote a significant part of the paper to the analysis of preconditioning methods. We show that block preconditioning is equivalent to the elimination of variables, and that the elimination can be performed directly on the non-linear system.

In a different direction, this paper also relaxes the reliance on the Morel formula-tion that requires that one identifies a priori a set of principal and secondary species. We show that the coupled problem can be formulated in the same way as before, but that the various totals involved can be defined in a more intrinsic way. This may be seen as a particular case of the general reduction mechanism of Knabner et al [START_REF] Hoffmann | A parallel global-implicit 2-D solver for reactive transport problems in porous media based on a reduction scheme and its application to the momas benchmark problem[END_REF][START_REF]A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions[END_REF], but we believe the simpler formulation may be of interest to practitioners. The formulation also leads to a method for solving the chemical system that makes uses of orthogonal matrices and may lead to better conditioned Jacobians (this remains to be investigated).

In this work, we consider only a simplified physical and chemical setup (one phase flow, no mineral reactions, no kinetic reactions) so as to concentrate on the numerical issues related to the coupling between transport and chemistry. Generalization to more realistic situations (2D and 3D geometries, kinetic reactions, presence of mineral species) will be the topic of future studies.

An outline for the rest of this paper is as follows: we set up the mathematical model in section 2, and show how to reduce the problem by eliminating the chemical concentrations. Numerical methods for solving the local chemical equilibrium problem and the advection-diffusion equations for transport are the topic of sections 3.1 and 3.2 respectively. Section 4 deals with the formulation of the coupled problem, and its solution by a Newton-Krylov method. Preconditioners for the linear system and elimination methods for the non-linear system are detailed in section 5. Finally, in section 6 the methods are validated on two test cases, including the 1D MoMaS reactive transport benchmark, which is a fairly difficult test case for reactive transport [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF].

2. Reactive transport model. We consider a set of species subject to transport by advection and diffusion and to chemical reactions in a porous medium occupying a domain Ω ⊂ R d (d = 1, 2, 3). [-] The chemical phenomena involve both homogeneous and heterogeneous reactions. Homogeneous reactions, in the aqueous phase, include water dissociation, acid-base reactions and redox reactions, whereas heterogeneous reactions occur between the aqueous and solid phases, and include surface complexation, ion exchange and precipitation and dissolution of minerals (see [START_REF] Appelo | Geochemistry, Groundwater and Pollution[END_REF] for details on the modeling of specific chemical phenomena). Accordingly, we assume there are N s mobile species (X j ) j=1,...,N S in the aqueous phase and Ns immobile species in the solid phase ( Xj ) j=1,..., NS , and that there are N r homogeneous reactions, and Nr heterogeneous reactions.

In this work, we only consider equilibrium reactions, which means that the chemical phenomena occur on a much faster scale than transport phenomena, see [START_REF] Rubin | Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions[END_REF]. This assumption is justified for aqueous phase and ion-exchange reactions, but may not hold for reactions involving minerals. Such reactions should be modeled as kinetic reactions, with specific rate laws [START_REF] Bear | Theory and Applications of Transport in Porous Media[END_REF][START_REF] Fan | A fully-coupled flow-reactive-transport formulation based on element conservation, with application to CO 2 storage simulations[END_REF].

We can write the chemical system as

Ns j=1 (S cc ) ij X j 0 i = 1, . . . , N r homogeneous reactions, Ns j=1 (S cc ) ij X j + Ns j=1 (S cc ) ij Xj 0 i = 1, . . . Nr heterogeneous reactions, or in condensed form (2.1) S X X = S cc 0 S cc S cc X X 0 0 .
We denote by S the stoichiometric matrix, with the sub-matrices S cc ∈ R Nr×Ns , S cc ∈ R Nr×Ns and S cc ∈ R Nr× Ns . We assume that both the global stoichiometric matrix S and the "aqueous" stoichiometric matrix S cc are of full rank. As there will usually be more species than reactions, this just means that all reactions are "independent" (in the linear algebra sense !), and that this is also true of the reactions in the aqueous phase. Each reaction gives rise to a mass action law, linking the activities of the species. For simplicity, we assume all species follow an ideal model, so that their activity is equal to their concentration. We denote by c j (resp. cj ) the concentration of species X j (resp. Xj ). It will be convenient to write the mass action law in logarithmic form, so for a vector c with positive entries, we denote by log c the vector with entries log c j . We then have

(2.2) S cc 0 S cc S cc log c log c = log K log K ,
where K ∈ R Nr and K ∈ R Nr are the equilibrium constants for their respective reactions.

Next, we write the mass conservation for each species, considering both transport by advection and diffusion and chemical reaction terms:

(2.3) φ∂ t c +Lc = S T cc r + S T cc r, φ∂ t c = S T cc r, in Ω × [0, T f ]
where L denotes the advection-diffusion operator (written in 1D):

L(c) = ∂ x (-D∂ x c + uc) ,
φ is the porosity (fraction of void in a Representative Elementary Volume available for the flow), u is the Darcy velocity (we assume here permanent flow, so that u is considered as known), D is a diffusion-dispersion coefficient and r ∈ R Nr and r ∈ R Nr are vectors containing the reaction rates. We assume that the diffusion coefficient is independent of the species. This is a strong restriction on the model, but one that is commonly assumed to hold [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF][START_REF] Kräutle | A new numerical reduction scheme for fully coupled multicomponent transport-reaction problems in porous media[END_REF][START_REF] Saaltink | A mathematical formulation for reactive transport that eliminates mineral concentrations[END_REF][START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF]. The model is completed by appropriate initial and boundary conditions. Since we assume that all reactions are at equilibrium, the reaction rates are unknown and we now show how they can be eliminated.

2.1.

Elimination of the reaction rates -the coupled problem. We follow the approach of Saaltink et al. [START_REF] Saaltink | A mathematical formulation for reactive transport that eliminates mineral concentrations[END_REF] (see also a more general approach in [START_REF] Kräutle | A new numerical reduction scheme for fully coupled multicomponent transport-reaction problems in porous media[END_REF][START_REF]A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions[END_REF]) by introducing a kernel matrix U such that U S T = 0, i.e. such that columns of U T form of basis for the null-space of S. This can be done in several ways (see the above references, and section 3.1). Here we outline how one can compute such a matrix with the same structure as that of S.

Lemma 2.1. Assume that the matrix S is as defined in equation (2.1), that it has full rank, and that the submatrix S cc also has full rank. Then there exists a kernel matrix

(2.4) U = U cc U cc 0 U cc , U cc ∈ R (Ns-Nr)×Ns , U cc ∈ R (Ns-Nr)× Ns U cc ∈ R ( Ns-Nr)× Ns such that U S T = 0. Proof.
The existence of a kernel matrix U is well known (see any standard text on linear algebra, such as [START_REF] Strang | Introduction to Linear Algebra[END_REF] or [START_REF] Meyer | Matrix Analysis and Applied Linear Algebra[END_REF]). The content of the lemma is that the kernel matrix can be chosen with a block upper triangular structure.

The proof proceeds by computing each block of the product, and showing that the blocks in U can be chosen as specified. This is true because both C cc and C cc are of full rank (the later because the full matrix S is block triangular).

With the kernel matrix constructed in lemma 2.1, we can eliminate the reaction terms in equation (2.3). We start by defining the total analytic concentration for the mobile and immobile species respectively (these are the same as various total quantities defined in the classical survey by Yeh and Tripathi [START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF])

(2.5) T T = U cc U cc 0 U cc c c .
We also define the total mobile and immobile concentrations for the species in the aqueous phase

(2.6) C = U cc c, C = U cc c,
so that the total concentrations are given by (2.7)

T = C + C = U cc c + U cc c.
These new unknowns are all of dimension N c , where N c = N s -N r . We now multiply system (2.3) on the left by U . Because of our assumption that D is the same for all chemical species, multiplication by U commutes with the differential operator, and the system can be rewritten as (2.8) φ∂ t C + φ∂ t C + LC = 0, φ∂ t T = 0.

The coupled system consists of the (N s -N r ) + ( Ns -Nr ) conservation PDEs and ODEs (2.8), together with the N r + Nr mass action laws (2.2) and the relations connecting concentration and totals (2.6) and the second line of (2.5), for the N s + Ns concentrations and 2(N s -N r ) + Ns -Nr totals. Note that the ODEs for T are decoupled from the rest of the system. In section 4, we show how to eliminate the individual concentrations from the system, so that only the totals C and C remain as unknowns.

Remark 2.1. The formulation for the chemical system given in equations (2.2) and (2.5) generalizes the well known Morel formulation [START_REF] Morel | Principles and Applications of Aquatic Chemistry[END_REF], where a set of "principal" species is identified, and the remaining "secondary" species are written in terms of the principal ones.

The stoichiometric matrix S is split naturally in blocks, the mass action laws allow the elimination of the secondary unknowns, and the conservation laws lead to a non-linear system of equations for the principal species.

In the next section, we show that the "local chemical system" (2.2), (2.5) can be solved (for given (T, T )) in a similar way, without having first to explicitly identify the principal and secondary species.

3. Methods for solving chemistry and transport. We briefly recall in this section how the chemical equilibrium system, and the transport equation are solved.

3.1.

The chemical equilibrium problem. The subsystem formed by the mass action laws (2.2) and the definition of the totals (2.5) is a closed system that enables computation of the individual concentrations c and c given the totals T and T . This is actually the same system that would be obtained in a closed chemical system where now T and T would be known as the input total concentrations. This subsystem is what we call in the following the chemical equilibrium problem. It is a small (its size is the number of species) nonlinear system, that is notoriously difficult to solve numerically (see for example [START_REF] Carrayrou | New efficient algorithm for solving thermodynamic chemistry[END_REF][START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF][START_REF] Hoffmann | A parallel global-implicit 2-D solver for reactive transport problems in porous media based on a reduction scheme and its application to the momas benchmark problem[END_REF][START_REF] Marinoni | Thermodynamic equilibrium solutions through a modified Newton Raphson method[END_REF]).

It will be convenient in this subsection to temporarily ignore the distinction between mobile and immobile species. We thus define the vectors

ξ = c c ∈ R N ξ , κ = log K log K ∈ R Nκ τ = T T ∈ R (N ξ -Nκ)
(N ξ = N s + Ns is the total number of species and N κ = N r + Nr is the total number of reactions), and write the problem for the chemical equilibrium as

(3.1) S log ξ = κ, U ξ = τ,
with the stoichiometric matrix S ∈ R Nκ×N ξ , and the kernel matrix

U ∈ R (Nκ-N ξ )×N ξ is such that U S T = 0.
It has been shown by several authors [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF][START_REF] Carrayrou | New efficient algorithm for solving thermodynamic chemistry[END_REF][START_REF] Hoffmann | A parallel global-implicit 2-D solver for reactive transport problems in porous media based on a reduction scheme and its application to the momas benchmark problem[END_REF] that taking the logarithms of the concentrations as new unknowns in (3.1) is beneficial from the numerical point of view, as it automatically ensures that all concentrations will be positive (which might otherwise be difficult to enforce), and it makes all the unknowns to be of comparable size (in some cases, such as redox reactions, the concentrations have been seen to vary over several orders of magnitude). We take the same convention to define the exponential of a vector as for the logarithm: for a vector z the vector exp(z) is defined so that its ith entry is exp(z i ). By defining z = log(ξ) ∈ R N ξ , system (3.1) thus becomes

(3.2) Sz = κ, U exp(z) = τ,
In order to solve system (3.2), we take a cue from the solution of constrained least squares problem (see [8, chap. 5], and also [START_REF] Erhel | Characterizations of Solutions in Geochemistry: Existence, Uniqueness and Precipitation Diagram[END_REF] in the same context): we consider the first equation as a constraint, and determine its general solution as the sum of a particular solution, and an unknown of smaller dimension, that will be found by substituting it into the second equation. Following the references above, we first compute a QR decomposition of S T , as

Q T S T = R 0 , with Q = (Q 1 , Q 2 ) ∈ R N ξ ×N ξ
orthogonal and R ∈ R Nκ×Nκ upper triangular (and invertible, as S is assumed to be of full rank). The columns of Q 1 ∈ R N ξ ×Nκ form an orthonormal basis of range S T , while those of Q 2 ∈ R N ξ ×(N ξ -Nκ) form an orthonormal basis of Ker S. A first consequence is that one may choose U = Q T 2 in the second equation of (3.2). Now, any z ∈ R N ξ may be written uniquely as z = Q 1 y 1 + Q 2 y 2 , with y 1 ∈ R Nκ and y 2 ∈ R N ξ -Nκ . If we substitute this expression in the first equation of (3.2), we obtain

R T y 1 = κ, which determines y 1 uniquely. We set b 1 = Q 1 R -T κ, so that z = b 1 + Q 2 y 2 ,with b 1 known.
We now use the expression for z in the second equation of (3.2), and we obtain a non-linear system of size N ξ -N κ for y 2 ∈ R N ξ -Nκ :

(3.3) H(y 2 ) = def Q T 2 exp(b 1 + Q 2 y 2 ) -τ = 0.
This is the system that is to be solved numerically, and that is analogous to the nonlinear system for the principal species obtained via the Morel formalism. The number of unknowns has been reduced from the number of species to the number of species minus the number of reactions, that is the number of principal species. As a side remark, let us notice that the Jacobian matrix for this reduced chemical system has the form

J c = Q T 2 diag(c) Q 2 , c = b 1 + Q 2 y 2
, and this form is again the same as that obtained from the Morel formulation, with the difference that here matrix Q 2 is orthogonal. This may be important as the Jacobian matrices occuring in the solution of the chemical equilibrium problem may be very ill conditioned, as shown in [START_REF] Machat | Comparison of linear solvers for equilibrium geochemistry computations[END_REF]. Thus using orthogonal matrices limits the ill-conditioning to the diagonal matrix diag(c) where it is harmless.

To solve the reduced chemical problem (3.3), we used a variant of Newton's method. As is well known Newton's method does not always converge in practice and, especially for a code that is designed to be used in coupling applications, it is essential to ensure that the solver "always" works. We have found that using a globalized version of Newton's method (using a line search, cf. [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF]) was effective in making the algorithm converge from an arbitrary initial guess.

We now return to the context of solving the coupled problem, where it is important to distinguish between mobile and immobile species. Indeed, what is needed is the partition of the species between their mobile and immobile forms, rather than the individual concentrations (though they are still needed as intermediate quantities). The totals can be computed a posteriori using their definitions in equation (2.6).

It will be convenient to condense the chemical sub-problem by a function

(3.4) ψ C : R N c → R N c T → ψ C (T ) = C = U cc c
where c is obtained by solving the chemical problem (3.2), given T (and T , which we take as a constant) for ξ, and computing C as indicated above.

Transport model.

In this section we denote by c a generic unknown concentration. When the transport system is solved in the context of a coupled problem, the role of c will be played by C as defined in Section 2, cf. equation (2.8). The transport of a single species through a 1D porous medium Ω =]0, L[ is governed by the advection-diffusion equation (cf. (2.3)):

(3.5) φ∂ t c + ∂ x (-D∂ x c + uc) = φ∂ t c + ∂ x j = φq, 0 < x < L, 0 < t < T f ,
implicitly defining the flux j.

The initial condition is c(x, 0) = c 0 (x) and, in view of the applications, the boundary conditions are a Dirichlet condition (given concentration) c(0, t) = c d (t) at the left boundary (x = 0) and zero diffusive flux ∂c ∂x = 0 at the right boundary (x = L). We could easily take into account more general boundary conditions.

Discretization in space.

We treat the space and time discretization separately, as we will use different time discretizations for the different parts of the transport operator.

For space discretization we use a cell-centered finite volume scheme [START_REF] Eymard | Finite volume methods[END_REF].

The interval ]0, L[ is divided into N h intervals [x i-1 2 , x i+ 1 2 ] of length h i , where x 1 2 = 0, x N h + 1 2 = L. For i = 1, .
., N h we denote by x i the center and x i+1/2 the extremity of the element i. We denote by c i , i = 1, .., N h the approximate concentration in cell i.

We split the flux in equation (3.5) as the sum of a diffusive flux j d = -D ∂c ∂x and an advective flux j a = uc. We then integrate equation

(3.5) over a cell ]x i-1/2 , x i+1/2 [, to obtain (3.6) φ i h i dc i dt + j d,i+ 1 2 + j a,i+ 1 2 -j d,i-1 2 -j a,i-1 2 = h i φ i q i , i = 2, . . . , N h .
Our flux approximations come from finite differences. For the diffusive flux, we use harmonic averages for the diffusion coefficient (as used in mixed finite element methods) :

(3.7) j d,i+ 1 2 = -D i+ 1 2 c i+1 -c i h i+ 1 2 with D i+ 1 2 = 2D i D i+1 D i + D i+1 , D 1 2 = D 1 D N h + 1 2 = D N h and h i+ 1 2 = h i + h i+1 2
For the advective flux, we use an upwind approximation, so that (assuming for simplicity that u > 0), j a,i+ 1 2 = uc i These approximations are corrected to take into account the boundary conditions, both at x = 0 and at x = L. We give a matrix formulation, keeping time continuous for now. Since we will be using different discretizations for the diffusive and advective parts (see next section), we keep the matrices for advection and diffusion separate. With the notation:

α i = D i+ 1 2 h i+ 1 2 + D i-1 2 h i-1 2 , β i = - D i-1 2 h i-1 2 , γ i = - D i+ 1 2 h i+ 1 2 , i = 2, . . . , N h -1,
we define the matrices (with appropriate modifications for the boundary terms)

A d = tridiag(β i , α i , γ i ), A a = tridiag(-u, u, 0)
as well as M = diag(φ i h i ).

The semi-discrete system can be written as

(3.8) M dc dt + (A d + A a )c = M q, with the initial condition c i = 1 hi x i+1/2
x 1-1/2 c 0 (x i ).

Time discretization.

As the transport operator contains both advective and diffusive terms, it makes sense to use different time discretization methods for the different terms. Specifically, the diffusive terms should be treated implicitly, and the advective terms are better handled explicitly. Similarly to what was done above, we discretize the interval [0, T f ] with a time step ∆t, which we take as a constant for simplicity, and we denote by c n i the (approximate) value of c i (n∆t), and by c n the corresponding vector.

We compare several methods for discretization in time. Fully implicit both the diffusive and advective terms are treated implicitly: at each time step, we solve a linear system

(M + ∆t(A d + A a ))c n+1 = M c n + M ∆tq n+1 .
Explicit advection and implicit diffusion the diffusive terms are treated implicitly, and advective terms are handled explicitly under a CFL condition. With this scheme, the system to be solved at each time step is

(M + ∆tA d )c n+1 = (M -∆tA a )c n + M ∆tq n+1 .
As it has an explicit part, the scheme just defined is stable under the CFL condition φ i u∆t ≤ max i h i . As this may be too severe a restriction (some of our applications require integration over a very large time interval), we use an operator splitting scheme proposed by Siegel et al. [START_REF] Siegel | Solution of the advection-dispersion equation using a combination of discontinuous and mixed finite elements[END_REF] (see also [START_REF] Hoteit | Nuclear waste disposal simulations: Couplex test cases: Simulation of transport around a nuclear waste disposal site: The COUPLEX test cases[END_REF]) that is both unconditionally stable, and has a good behavior in advection dominated situations. Splitting (Explicit advection and implicit diffusion), with sub-time steps : this scheme works by taking several small time steps of advection, controlled by CFL condition within a large time step of diffusion. Thus CFL impacts advection, and larger time steps can be taken for diffusion. More precisely, the time step ∆t will be used as the diffusion time step, it is divided into N a time steps of advection ∆t a such that ∆t = N a ∆t a where N a ≥ 1, the advection time step will be controlled by CFL condition. Note that taking a single advection step amounts to using the implicit-explicit method seen previously. We solve equation (3.5) over the time step [t n , t n+1 ] by first solving the advection equation φ∂ t c + ∂ x (uc) = 0 over N a steps of size ∆t a each, and then solve the diffusion equation φ∂ t c + ∂ x (-D∂ x c) = φq starting from the value at the end of the advection step. Advection step. Denote the intermediate times by t n,m , m = 0, . . . , N a , with

t n,0 = t n , t n,Na = t n+1 . Each interval [t n , t n+1 ] is then divided into N a intervals [t n,m , t n,m+1 ], m = 0, ...N a -1. Let c n,
m be the approximate concentration c at time t n,m and c n,0 = c n . We discretize the advection equation in time, using the explicit Euler method, we obtain

(3.9) M c n,m+1 = (M -∆t a A a )c n,m , m = 0, . . . , N a -1.
Diffusion step. The diffusion part is discretized by an implicit Euler scheme, starting from c n,Na :

(3.10) (M + ∆tA d )c n+1 = M c n,Na + M ∆tq n+1 .
For further use, we note that all three discretizations methods can be written in a similar way, as (3.11)

Ac n+1 = Bc n + M ∆tq n+1 ,
where the matrices A and b are defined in each case as:

Fully implicit A = M + ∆t(A d + A a ) and B = M , Explicit-implicit A = M + ∆tA d and B = M -∆tA a ,
Splitting A = M + ∆tA d and B is defined implicitly by (3.9) Alternatively, we may want to follow the pattern begun in section 3.1, and abstract one transport step as a mapping from c n to c n+1 under the action of the source q n+1 . We define (3.12)

ψ T : R N h × R N h → R N h (c, q) → ψ T (c, q) = A -1 (Bc + M ∆tq)
4. Methods for solving the coupled system. In this section, we discuss methods for solving the coupled transport-chemistry system. It may be difficult to compute and store the Jacobian matrix (it may be very large, and contains a block that is the inverse of the chemical solution operator, that may be difficult to compute). Consequently we advocate a Newton-Krylov approach, as this requires only the multiplication of the Jacobian matrix by a given vector, the Jacobian matrix itself does not have to be stored. By exploiting the block structure of the Jacobian matrix, the matrix-vector product can be computed in terms of the individual Jacobians.

The efficiency of the Newton-Krylov method rests on the choice of an adequate preconditioner (see [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF] for several illustrations). In this paper, we show that block preconditioners for the Jacobian can be related to physics-based approximations, and also that a block Gauss-Seidel preconditioner is closely related to an elimination method at the non-linear level.

4.1. Formulation of the coupled system. We start with the coupled system that was obtained at the end of section 2. It consists of the (transformed) transport equations (2.8), together with the mass action laws (2.2). They are linked by the definition of the transformed variables T, C and C in equations (2.6) and (2.7).

Because chemistry is local, we can eliminate the individual concentrations at each point by using the "chemical solution" operator ψ C defined in equation (3.4). This only leaves C, C, T and T as unknowns that are solution of the following system (4.1)

φ∂ t C + φ∂ t C + LC = 0, φ∂ t T = 0, T = C + C, C = ψ C (T ).
As is clear from the second equation, T is constant in time (we will see that it is simply the cationic exchange capacity in the example below). To simplify the notation, the equation for T will be omitted in the sequel, and it has been dropped from the definition of ψ C . The next step is to discretize the transport equations in space and time, using any of the methods that were discussed in section 3.2. For ease of notation, we denote the discrete unknowns (vectors of size N h ) by the same letter in bold typeface as their continuous counterpart. Moreover, all quantities at time t n will be indicated by a superscript n.

We make use of a notational device, inspired from Matlab, that was introduced in [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF] to take into account the dependence of the unknowns both on the grid cell index and the chemical species index. For a block vector u ij , where i ∈ [1, N c ] represents the chemical index and j ∈ [1, N h ] represents the spatial index, we denote by • u :,j the column vector of concentrations of all chemical species in grid cell j.

• u i,: the row vector of concentrations of species i at all grid points; The unknowns are thus numbered first by chemical species, then by grid points, so that all the unknowns for a single grid point are numbered contiguously. As proposed by Erhel et al. [START_REF] Erhel | Solving partial differential algebraic equations and reactive transport models[END_REF], we make use of the Kronecker product and the vec notation (see for instance [START_REF] Golub | Matrix Computations[END_REF]). Given

V = V 1 , V 2 , . . . , V Nc ∈ R Nx×Nc (V i denotes the ith column of V ) we denote by vec(V ) the vector in R NxN c such that vec(V ) =      V 1 V 2 . . ., V Nc      .
We also extend the solution operator ψ C (defined in (3.4) to operate on the global vector:

(4.2) Ψ C : R Nc×N h → R Nc×N h T → Ψ C (T), with Ψ C (T) :,j = ψ C (T :,j ), for j = 1, . . . , N h
and define the right hand side vector by b n = BC n + M Cn . With these conventions, the discretized version of the coupled system (4.1) becomes

(4.3) A(C n+1 i,: ) T + M ( Cn+1 i,: ) T -(b n i,: ) T = 0, i = 1, . . . , N c , T n+1 ij -C n+1 ij -Cn+1 ij = 0, i = 1, . . . , N c , j = 1, . . . , N h , Cn+1 :,j -Ψ C (T n+1 :,j ) = 0, j = 1, . . . , N h .
This is a non-linear system of equations to be solved at each time step for the three unknowns (C n+1 , T n+1 , Cn+1 ) ∈ R 3NcN h , defined by the function f : R

3NcN h → R 3NcN h , such that (4.4) f   C T C  =   (A ⊗ I)C + (M ⊗ I) C -b n T -C - C C -Ψ C (T)   = 0 4.2.
Solution of the coupled system by Sequential Iterative Approach. The most classical method for solving the coupled problem (4.4) is the Sequential Iterative Approach (SIA) that consists of separately solving the transport equations and the chemical equations cf. [START_REF] Carrayrou | Operator-splitting procedures for reactive transport and comparison of mass balance errors[END_REF][START_REF] Lagneau | HYTEC results of the MoMas reactive transport benchmark[END_REF][START_REF] Saaltink | A mathematical formulation for reactive transport that eliminates mineral concentrations[END_REF][START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF].

It will be convenient to extend the solution operator for transport ψ T (defined in (3.12)) to operate on global vectors (4.5)

Ψ T : R Nc×N h × R Nc×N h → R Nc×N h (C, Q) → Ψ T (C, Q), with Ψ T (C, Q) i,: = ψ T (C i,: , Q i,: ), for i = 1, . . . , N c .
The solution operators Ψ C and Ψ T emphasize that both the transport and chemistry steps can be seen as black boxes. This is the basis for the approach followed in [START_REF] Heredia | TReacLab: an object oriented implementation of non intrusive operator splitting methods to couple independent transport and geochemical software[END_REF], where a flexible code that allow switching the transport and chemistry components is presented. We show in the next section that this black box approach is not restricted to the SIA.

At each time step, we iterate on the transport and chemistry problem. More precisely, for each iteration k, we first solve the transport equations: (4.6)

C n+1,k+1 = Ψ T C n , M ( Cn+1,k -Cn ) ∆t = -(A -1 ⊗ I) (M ⊗ I) Cn+1,k -b n
for C n+1,k+1 (there is one transport equation for each species). The new mobile total concentrations is added to the immobile concentrations Cn+1,k (from the previous iteration) to obtain T n+1,k+1 , and this last concentration provides the data for the local chemical problems

(4.7) Cn+1,k+1 = Ψ C (T n+1,k+1
).

The iterations are stopped when the difference between the solutions in the iterations k and k + 1 is small enough, relative to a specified tolerance

(4.8) C n+1,k+1 -C n+1,k C n+1,k+1 + Cn+1,k+1 -Cn+1,k Cn+1,k+1 <
It is known that in order both to control the errors due to the splitting, and to ensure convergence, this method requires a small time step [START_REF] Macquarrie | Reactive transport modeling in fractured rock: A state-of-the-science review[END_REF][START_REF] Saaltink | A comparison of two approaches for reactive transport modelling[END_REF][START_REF] Valocchi | Accuracy of operator splitting for advection-dispersionreaction problems[END_REF].

Solution of the non-linear system by a Newton-Krylov method.

In the geochemical literature, the SIA method is known as an operator splitting approach, but it is more properly a block Gauss-Seidel method on the coupled system. This suggests that other, potentially more efficient, methods could be used to solve system (4.4). A natural candidate is Newton's method. Since the system is large, and also because it involves implicitly defined functions, we turn to the Newton-Krylov variant.

Recall that at each step of the "pure" form of Newton's method for solving f (Z) = 0, one should compute the Jacobian matrix J = f (Z k ), solve the linear system (4.9)

J δZ = -f (Z k ),
usually by Gaussian elimination, and then set Z k+1 = Z k + δZ. In practice, one should use some form of globalization procedure in order to ensure convergence from an arbitrary starting point. If a line search is used, the last step should be replaced by Z k+1 = δZ + λZ k , where λ is determined by the line search procedure [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF].

The Newton-Krylov method (see [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF][START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF] and [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physics-based preconditioning to biogeochemical transport[END_REF], to which our work is closely related) is a variant of Newton's method where the linear system (4.9) that arises at each step of Newton's method is solved by an iterative method (of Krylov type), for instance GMRES [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. As the linear system is not solved exactly, the convergence theory for Newton's method does not apply directly. However, the theory has been extended to the class of Inexact Newton methods, of which the Newton-Krylov methods are representatives. The theory leads to a practical consequence, by giving specific strategies for the forcing term, that is the tolerance within which the linear system has to be solved (see [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF] or the short discussion in [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF]).

The main advantage of this type of method is that the full (potentially very large) Jacobian is not needed, one just needs to be able to compute the product of the Jacobian with a vector. As this Jacobian matrix vector product is a directional derivative, this leads to Jacobian free methods, where this product is approximated by finite differences. However, for the system considered here, this has several drawbacks: in addition to its inherent inaccuracy, computing the derivative by finite difference entails solving one more chemistry problem. It was shown in [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF] how this computation could be performed exactly. This is both cheaper and more accurate. Moreover, the above reference also shows how the natural block structure present in the coupled system (4.3) can be exploited to efficiently compute the residual. This is what is done in the present work, and we now give some details.

To compute the Jacobian matrix times vector product, start from equation (4.4). The Jacobian matrix of f has the following block form (4.10)

J f =   A ⊗ I 0 M ⊗ I -I I -I 0 -J C I   ,
where J C is a block diagonal matrix, whose jth block is ψ C (T :j ), for j = 1, . . . , N h , and the action of the Jacobian on a vector

v = v T C , v T T , v T C
T can be computed as

(4.11) J f   v C v T v C   =    (A ⊗ I)v C + (M ⊗ I)v C -v C + v T -v C v C -J C v T    .
Of course, the Kronecker product is just a notational device, and the matrices A ⊗ I or M ⊗ I) are never actually formed. Instead, we use a well know property of the Kronecker product (see [START_REF] Golub | Matrix Computations[END_REF]) :

(4.12) (A ⊗ I) vec(V ) = vec(V A T ),
and this just requires to multiply the matrix A by the concentration vector for each species.

Since the Krylov solvers can stagnate, resulting in slow convergence, possible strategies for preconditioning the linear system will be investigated in the next section.

Non linear and linear preconditioning.

Since the Jacobian matrix is not explicitly computed (it will just be used as a theoretical device in what follows), the only available options for preconditioning the system are those that respect the block structure of the matrix. We introduce two preconditioners derived from classical block-iterative methods, and we show that these methods have strong links to the SIA method from section 4.2, and also to a non-linear elimination method, to be introduced in section 5.2.

5.1.

Block preconditioning for the coupled system. Using a preconditioner means that instead of solving, at each Newton iteration, the system (5.1)

J f δZ = -f (Z k )
one solves one of the (mathematically equivalent) systems: right preconditioning J f P -1 δy = -f (Z k ), and then δZ k = P -1 δy, left preconditioning P -1 J f δZ = -P -1 f (Z k ). The matrix P is called a preconditioner, and it should be chosen so as to fulfill the often conflicting goals :

• P should be close to J f so that GMRES converges faster for the preconditioned system than for the original one,

• P is "easy" to invert, so that each iteration is not much more expensive than an iteration for J f alone. As already mentioned, in the context of the coupled problem (4.4), the entries of the Jacobian matrix J f are not known explicitly, so that methods that depend on the entries of the matrix (such as incomplete factorization methods) cannot be used. On the other hand, we know the block structure of J f (see (4.10)), and we can invert the upper left diagonal block (of course, when we write "invert" we really mean "solve a linear system with". We never actually form the inverse. See also remark 5.1 below). These two facts can be exploited by resorting to block preconditioners that respect the bock structure of the Jacobian. Specifically, we investigate block methods derived from classical iterative methods, namely Jacobi and Gauss-Seidel. Block Jacobi preconditioning The preconditioning matrix P for block Jacobi is:

P BJ =   A ⊗ I 0 0 0 I 0 0 0 I   ,
so that the action of the left-preconditioned matrix

P BJ J f on a vector v = (v T C , v T T , v T C ) T is : (5.2) P -1 BJ J f v =   v C + ((A -1 M ) ⊗ I)v C v C + v T -v C -J C v T   ,
Block Gauss-Seidel preconditioning Here, the preconditioning matrix P BGS for block Gauss-Seidel and its inverse are :

P BGS =   A ⊗ I 0 0 -I I 0 0 -J C I   , P -1 BGS =   A -1 ⊗ I 0 0 A -1 ⊗ I I 0 J C (A -1 ⊗ I) J C I   ,
so that the action of the left-preconditioned matrix P BGS J f on a vector

v = (v T C , v T T , v T C ) T is : (5.3) P -1 BGS J f v =   v C + (A -1 M ) ⊗ I v C v T -v C + (A -1 M ) ⊗ I v C v C -J C v C + J C (A -1 M ) ⊗ I v C .   ,
Here again, as in (4.11), neither A -1 nor the Kronecker products are computed. Rather, to compute w = (A -1 M ) ⊗ I v C , for a given vector v C ∈ R NcN h , one defines W ∈ R Nx×Nc as the solution of (5.4)

AW = M V T C where v C = vec(V C ), then let w = vec(W ).
This means that the action of the preconditioner can be computed by solving a transport step for each chemical species, and a multiplication by the Jacobian of the chemical operator (which in turns requires solving a linearized chemical problem for each grid cell). These are also the building blocks for the SIA formulation, and this shows that the fully coupled approach can be implemented at roughly the same cost per iteration as the SIA approach.

Remark 5.1. In both cases, notice that the preconditioning step involves the inverse of the transport block. In our 1D case, this is not a difficulty, as this block can be easily inverted. When we move to 2D or 3D problems, this step would become more problematic. However, the preconditioners could then be defined by replacing the matrix A in the definition of J f by the application of a (spectrally equivalent) matrix, such as several iterations of a multigrid solver. A multigrid method such as that proposed recently in [START_REF] Southworth | Reduction-based Algebraic Multigrid for Upwind Discretizations[END_REF] would be particularly appropriate.

Remark 5.2. In both cases (exact or inexact transport solver), one should expect mesh independent convergence, as the operator being solved by GMRES becomes bounded independently of the mesh. This expectation will be confirmed by the numerical results in section 6.1.5.

Elimination of unknown as non-linear preconditioning.

In this section, we consider alternative solution strategies, based on eliminating some of the unknowns from the system (4.4). We group these strategies under the heading of "non-linear preconditioning" because the elimination is done before the linearization, but it must be noted that we take advantage of the fact that the first (block) equation of (4.4) is linear and can be solved species by species. In this context, this strategy should be reminiscent of the "non-linear preconditioning" as introduced by Cai and Keyes [START_REF] Cai | Nonlinearly preconditioned inexact Newton algorithms[END_REF], where block of unknowns are eliminated locally. It also provides an extension to multi-component transport of previous work by Kern and Taakili [START_REF] Kern | Linear and nonlinear preconditioning for reactive transport[END_REF] that deals with preconditioning a model with one species undergoing sorption.

We first eliminate T from the original system (4.4), leading to a system with only C and C as unknowns

g C C = 0 with (5.5) g C C = (A ⊗ I)C + (M ⊗ I) C -b n C -Ψ C (C + C) .
The Jacobian of the new system is

J g = A ⊗ I M ⊗ I -J C I -J C .
and the Jacobian by vector product is :

(5.6)

J g v C v C = (A ⊗ I)v C + (M ⊗ I)v C -J C (v C + v C ) + v C ,
One can even go one step further, by eliminating both the unknowns T and C, to obtain a system with C as the single unknown:

(5.7) h( C) = C -Ψ C Ψ T C n , M C -Cn ∆t + C = 0.
This equation is presented in fixed point form, and solving it by fixed point method recovers the SIA method described in Section 4.2. But equation (5.7) can also be solved by Newton's method. We presented a Jacobian-free Newton-Krylov method in [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF]. Its main advantage is to require that one be able to evaluate the nonlinear residual function h, and the Jacobian times vector product is approximated by finite difference. Since it involves solving an additional chemical problem, this step is very expensive. In this work, we use the exact Jacobian of h, which requires that we "open the black box", to look in more detail at the structure of the Jacobian. To do this, we rewrite the equation above by making use of the precise definition of the transport operator Ψ T :

(5.8)

h( C) = C -Ψ C (A -1 ⊗ I) b n -(M ⊗ I) C + C = 0.
If we denote by J T = I -(A -1 M ) ⊗ I the Jacobian matrix of Ψ T , we easily see that the Jacobian by vector product for h is:

(5.9)

J h v = v -J C J T v = v -J C v + J C (A -1 M ) ⊗ I v.
Both the alternative formulation and its Jacobian involve the resolution of transport at each Newton step. One again expects that convergence of the linear solver will become independent of the mesh size [START_REF] Kern | Linear and nonlinear preconditioning for reactive transport[END_REF]. This will be confirmed in Section 6.

Links between block preconditioning and elimination.

The main advantage of the SIA approach is that it allows the reuse of existing software modules for solving transport and chemistry. This is an important practical issue, as these modules may have been developed independently, even by different groups. This advantage is offset by the possibly slower convergence when compared to GIA, see for example [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physics-based preconditioning to biogeochemical transport[END_REF]). The "h" method proposed in section 5.2 (keeping only C as unknown) belongs to the GIA family, in the sense that it solves both chemistry and transport as a coupled system. This method is in the same spirit as the one proposed by Knabner and his group [START_REF] Hoffmann | A parallel global-implicit 2-D solver for reactive transport problems in porous media based on a reduction scheme and its application to the momas benchmark problem[END_REF][START_REF] Kräutle | A new numerical reduction scheme for fully coupled multicomponent transport-reaction problems in porous media[END_REF][START_REF]A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions[END_REF]. The reduction method proposed in these papers is more general than the formulation in section 4.1, but their "resolution function" is identical to the chemical solution operator Ψ C .

However, the method aims at keeping some of the advantages of SIA: it allows to keep transport and chemistry as separate modules. Equation (5.8) shows that the evaluation of the residual when solving the system with h requires solving a transport step for each species, and then solving a local chemical equilibrium system at each grid cell. It was shown in [START_REF] Amir | A global method for coupling transport with chemistry in heterogeneous porous media[END_REF][START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF] that this is sufficient, provided one accepts to approximate the Jacobian matrix by vector product by a finite difference quotient. However, the present paper makes the point that it is both cheaper and more accurate to compute this matrix-vector product exactly. This requires more cooperation from the chemical solver, as one needs to access its Jacobian computation;

We also give some indications for comparing the cost of the "h" method with SIA. At each iteration of SIA, one has to solve one transport problem per component, and one chemical equilibrium problem per grid cell. For the "h" method, the same cost is incurred at each Newton iteration, and one has to add the cost of the Jacobian matrix vector product at each GMRES iteration. This translates to one additional transport problem for each component, and one linearized chemical problem for each cell. We neglect the overhead of constructing the Arnoldi basis in GMRES, as the number of iterations is expected to be small (as confirmed by the results in section 6.1.5). So the number of transport problem to be solved is the sum of the number of Newton and GMRES iterations, but the number of chemical equilibrium problems is only the number of Newton iterations. One can hope (and again this is confirmed by our numerical experiments) that the number of Newton iterations will remain small. Additionally, it has been observed that the SIA method needs a small time step to reduce the splitting errors (see [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF][START_REF] Lagneau | HYTEC results of the MoMas reactive transport benchmark[END_REF]).

It is more difficult to comment on the cost of GIA, as this will depend in a crucial way on the efficiency of the solver. With proper care and effort, DSA methods can be made very efficient [START_REF] Hammond | Evaluating the performance of parallel subsurface simulators: An illustrative example with PFLOTRAN[END_REF]. The results in [START_REF] Hoffmann | A parallel global-implicit 2-D solver for reactive transport problems in porous media based on a reduction scheme and its application to the momas benchmark problem[END_REF] also show that a method based on a structure similar to that of the "h" method was among the fastest on the MoMaS benchmark.

It is interesting to note that both the SIA method and the "h" function method from section 5.2 can be interpreted in terms of the block preconditioners introduced in section 5.1.

Indeed, the elimination method can be interpreted as a (linear) change of variables, given by (5.10)

       C = C + (A -1 M ) ⊗ I C T = T -(I -A -1 M ) ⊗ I C C = C , whose matrix B =   I 0 (A -1 M ) ⊗ I 0 I -(I -A -1 M ) ⊗ I 0 0 I  
is block triangular, so that the transformation is easily inverted. The transformed system

f    C T C   = f   C T C 
takes the simple form

(5.11) f    C T C   =     (A ⊗ I) C -b -C + T C -Ψ C T + (I -A -1 M ) ⊗ I C     = 0.
and the last equation is obviously identical to (5.8), after having eliminated the first two unknowns.

Because the transformation used is linear, and because Newton's method is invariant under a linear transformation, the iterates between the original and the transformed system will be related by the same transformation (provided the initial guesses are related similarly). Note that this will not necessarily be true for an inexact Newton's method, such as the Newton-Krylov method used in this work. However, Deuflhard points out [19, sec. 2.2.4] that because the Newton residual is invariant under affine transformations, GMRES is the natural choice for solving the linear system arising at each Newton iteration. Moreover, left preconditioning can be used with GMRES provided the preconditioning matrix

• is kept constant throughout the Newton iterations,

• and is incorporated in the convergence monitoring criteria.

Note that the Jacobian of the transformed system is the block triangular matrix (5.12)

J f =   A ⊗ I 0 0 -I I 0 0 -J C J h  
where J h = I -J C (I -A -1 M ) ⊗ I , defined in (5.9), is actually the Schur complement of J f with respect to its first two variables. This confirms the close link between the non-linear elimination method from section 5.2 and block Gauss-Seidel preconditioning in (5.3): as noted above, the Jacobian of h is exactly the Schur complement of the Jacobian of f . In both cases, one needs to compute (J C ((A -1 M ) ⊗ I))v. Of course, we are once again taking advantage of the fact that one of the operators is linear ! An interesting consequence of the change of variables is that, because J f and the change of variable matrix B are both block triangular, differentiating the identity f = f • B leads to a block triangular factorization of J f , (5.13)

J f = J f B =   A ⊗ I 0 0 -I I 0 0 -J C S     I 0 (A -1 M ) ⊗ I 0 I -(I -A -1 M ) ⊗ I 0 0 I   .
This factorization could be used as a basis for constructing efficient preconditioners for J f , much as in the spirit of [START_REF] Howle | Block preconditioners for coupled physics problems[END_REF], [START_REF] Ipsen | A note on preconditioning nonsymmetric matrices[END_REF] (see also [START_REF] Murphy | A note on preconditioning for indefinite linear systems[END_REF]). The matrix B is upper triangular with unit diagonal, so that all its eigenvalues are equal to 1. As noted in [START_REF] White | Block-preconditioned Newton-Krylov solvers for fully coupled flow and geomechanics[END_REF], under these conditions, Krylov subspace methods for the preconditioned system J -1 f J f converge in 1 iteration, giving J f as a perfect preconditioner.

Because it contains S as its (3, 3) block, J f cannot be formed, let alone inverted. A first solution is to replace S by an approximation S that is easier to invert. The simplest choice is to take S = I, which gives the block Gauss-Seidel preconditioner from section 5.1. But actually, systems with S can be solved by an iterative method, as the matrix vector product Sv can be computed by proceeding as for the Gauss-Seidel preconditioner at the end of section 5.1 (see equation (5.4)).

The numerical results in section 6.1.5 (cf. figure 6.9) will confirm that the performance of SIA and block Jacobi on the one hand and the elimination method and block Gauss-Seidel on the other hand are very similar.

6. Numerical results. 6.1. MoMaS Benchmark : 1D easy advective case. The MoMaS Benchmark has been designed to compare numerical methods for reactive transport models in 1D and 2D. Different methods for coupling have been used to solve this benchmark. The definition has been published in [START_REF] Carrayrou | Reactive transport benchmark of MoMaS, Computat[END_REF] and the results of participants are compared in the synthesis article [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF].

The geometry of the test case is shown in Figure 6.1. For the 1D test case, the domain is heterogeneous and composed of two porous media A and B. Medium A is highly permeable with low porosity and low reactivity in comparison with medium B. Their physical properties are given in Table 6.1. The Darcy velocity is constant over the domain and is equal to φu = 5.5 10 -3 LT -1 .

The chemical reactions are summarized in Table 6.2. The 7 reactions involve 9 mobile species (in the aqueous phase) and 3 immobile species (in the solid phase). The characteristic feature of this chemical system is that it contains large stoichiometric coefficients that range from -4 to 4 and a large variation of equilibrium constants between 10 -12 and 10 35 .

Boundary and initial conditions are presented in Table 6.3. The simulation involves an injection period followed by a leaching period, so that the system is returned to its initial state. Note that the total for the immobile species is given in Table 6.1. [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF] indicate that one needs to refine the mesh around medium B if one is to obtain accurate results. For all test cases, we use a mesh such that h A = 4h B . The computations were carried out with the various methods presented in this paper. When used with the appropriate numerical parameters they all gave comparable results. The figures in this section were obtained with the h-method. Figure 6.2 shows profiles of the concentrations of several species. The left and middle images show concentrations of (a part of) mobile species C1 and immobile species S at an early time t = 10, whereas the right image shows concentration of aqueous species C 2 on a smaller interval during the leaching period, at t = 5010. The middle image highlights the effect of the heterogeneity at x = 1. Figure 6.3 shows elution curves, that is evolution of the concentrations at the end of the domain (x = 2.1) as a function of time, for a mesh with 384 elements. Notice that the evolution of C5 follows a fairly complex pattern. It turns out that the accuracy for both species X3 and C5 is quite sensitive to the mesh size used. We come back to this point in detail in the next subsection. The results shown on these figures are in good agreement with those showed by various groups in the comparison paper [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF].

X 1 X 2 X 3 X 4 S K C 1 0 -1 0 0 0 10 -12 C 2 0 1 1 0 0 1 C 3 0 -1 0 1 0 1 C 4 0 -4 1 3 0 10 -1 C 5 0 4 3 1 0 10 +35 CS 1 0 3 1 0 1 10 +6 CS 2 0 -3 0 1 2 10 -1
As mentioned above, some of the species show a very sensitive dependence on the mesh and time step used. It has proven necessary to use very fine meshes, as well as small time steps, in order to resolve these species accurately. We now study in more detail how the accuracy depends on the time and space meshes 6.1.2. Influence of spatial discretization. The evolution of species X3 and C5, as a function of time, exhibits unphysical oscillations. The origin of these oscillations has been explained in [START_REF] Lagneau | HYTEC results of the MoMas reactive transport benchmark[END_REF], and is due to the interplay between the very stiff reactions and the spatial discretizations. They should decrease as the mesh is refined. Figures 6.4 shows that this is indeed what happens as we increase the number of discretization points. A mesh with 384 points gives qualitatively correct results, but we have also made use of a finer meshes in order to obtain more accurate results.

The effect of the mesh size on the accuracy of the results is also shown by looking in detail at two specific species: mobile species C 1 and immobile species S, both at time t = 10. They both exhibit a sharp peak, and we focus on the accuracy with which the location and amplitude of the peak can be determined. These elements were part of the comparison criteria for the benchmark, and were examined in detail in [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF].

We also discuss the influence of the discretization scheme, by comparing the three schemes introduced in section 3.2 from the point of view of accuracy. Their relative efficiencies are compared in section 6.1.4. In figures 6.5 and 6.6, the curves are labeled by the total number of discretization points in space. The coarsest mesh, corresponding to n = 96, has ∆x = 0.025 in medium A and ∆x = 0.00625 in medium B. In both cases the time step was chosen as ∆t = ∆t c (advective time step), with a CFL=0.1, as given by the smaller space mesh in medium B. Note that since the CFL is fixed, each curve also corresponds to a different time step. Figure 6.5 compares the concentrations of species C 1 on the interval [0, 0.2], at time t = 10, as computed by the three discretization schemes on increasingly refined meshes. As soon as the mesh has sufficiently many points to successfully resolve the solution, all three schemes give identical solutions. On the other hand, one needs at least 384 points (and preferably 768) to obtain a satisfactory solution. Here also, the three schemes give identical solutions when the mesh is fine enough. This time, one needs at least 768 mesh points to obtain a converged solution. One should compare Figure 6 in [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF], and also refer to Table 4 there, where the location and amplitude of the peak are tabulated for all the methods used in the benchmark. We have obtained values of x = 0.0175 for the location of the peak, and S = 0.985 for its height. These values are in the range reported by the other teams, but are different form the "mean" values as reported in [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF]. They are however very close to the "reference" values found in [START_REF] Carrayrou | Looking for some reference solutions for the reactive transport benchmark of MoMaS with SPECY[END_REF]. The fully implicit method here there is only one time step, and no stability restriction. We have compared three time steps, corresponding to CFL coefficients of 0.1, 1 and 3 respectively. The explicit-implicit method this method also imposes a single time step, and in addition it is subject to the stability condition of the explicit method, so the time step is restricted to CF L ≤ 1. We compare 2 time steps, corresponding to CFL coefficients of 0.1 and 1 respectively. Figure 6.7 compares the concentration of species S as computed by the three methods, for the time step sizes chosen as explained above, for two different mesh resolutions.

For all three methods the location of the peak was correctly determined even for larger values of the time step and the mesh size, but its amplitude was only correctly estimated for the finer mesh size.

For both splitting and explicit-implicit schemes, it is not necessary to use a CFL of 0.1 (a CFL of 1 is enough) for the peak amplitude to reach the value 1, but what is needed is to refine the mesh (up to 1536 nodes). However, the fully implicit scheme needs both a small time step corresponding to a CFL of 0.1 and a fine mesh with 1536 nodes to obtain the same results as the other schemes. a 1D model, so that solving a linear system is a not very time consuming. The conclusions reached below would have to be re-examined for a 2D (and even more for a 3D !) model. We have also compared the effect of using the fully implicit method with a variable time step (albeit not an adaptive choice, a sequence of pre-computed time steps is used for different phases of the solution evolution). Figure 6.8 compares the CPU times required by the h-method with different time-discretization schemes as the space mesh is refined. The time step was chosen as follows (in all cases, the smallest mesh size was used): For the splitting method The diffusion time step is 3 times larger than the advection time step that respects a CFL coefficient of 1. For the explicit-implicit method One time step is used controlled by a CFL condition of 1. For the fully implicit method One time step is used corresponding to a CFL coefficient of 3 (no stability restriction). For the fully implicit (variable) method Variables time steps are used, as described in table 6.4. The time steps are chosen so as to have a small time step when strong variations happen due to the reactions (especially during injection and leaching period) and a large time step is used for the intervals that represent the steady state. As expected, using a variable time step results in large savings (while maintaining the accuracy). Among the 3 schemes with fixed time step, the explicit-implicit method is the most expensive, with the fully implicit and the splitting methods leading to comparable costs. Fig. 6.8. CPU time required by the time discretization methods as the space mesh is refined 6.1.5. Influence of preconditioning strategy. In this section, we compare the various preconditioning strategies discussed in section 5. Our main criteria will be the number of linear and non-linear iterations. We have not tried to optimize the inexact Newton strategy, but have just relied on the default choices as provided in the Newton-Krylov code (we use the nsoli code from the book by C. T. Kelley [START_REF]Solving Nonlinear Equations with Newton's Method[END_REF]). In all experiments, GMRES was used without restart, and with a maximum number of allowed iterations fixed at 40. Figure 6.9 shows how these numbers change as the mesh is refined. The various linear preconditioning strategies (applied to the coupled system) are compared with the elimination, or nonlinear preconditioning, strategy. As predicted, the nonlinear elimination strategy has the smallest number both for non-linear and for linear iterations. It also shows a behavior that is independent of the mesh size. The unpreconditioned method is unsurprisingly non scalable, at least as far as the linear iterations are concerned. The number of non-linear iterations grows only weakly with the number of mesh points. The same is true for Jacobi preconditioning.

Gauss-Seidel preconditioning, on the other hand, show only a modest increase in the number of linear iterations, and a behavior for the non-linear iterations that is in between that of the unpreconditioned and the elimination strategies.

As explained in section 5, a limitation of the elimination strategy is that it requires an exact solution of the transport step. In this case, the Gauss-Seidel preconditioner might prove useful: replacing the transport solve step by an approximation, such as several iterations of a multi-grid solver, should lead to a more efficient solution method, with similar convergence behavior. We plan to explore this strategy in a forthcoming work.

Last, figure 6.10 shows the time required by the various methods. Since the cost of the methods is comparable, the ordering is the same as that in the previous figure. It confirms the good efficiency of the elimination strategy, with Gauss-Seidel preconditioning as a distant second. We can try and summarize the relative performances of the various methods as follows:

• the original f formulation is not numerically scalable, neither at the non-linear level, nor at the linear level; • the block Jacobi preconditioner applied to f does not bring any improvement;

• the block Gauss-Seidel preconditioner improves the linear performance, but does have a significant non-linear effect (nor was it expected); • the g formulation improves the linear performance, but not as much as Gauss-Seidel preconditioning; • the h formulation, after elimination of the C unknowns is the only method that gives a convergence independent of the mesh size. Of course, these conclusions are more or less natural: methods f and g keep the ill-conditioning from the second order operator (except for g at the linear level). The elimination method, leads to a bounded operator on L 2 (at least formally), and is expected to give mesh independent convergence. This good performance of the elimination method, at least on the linear level, can be confirmed by looking at the field of values of the matrix J h . The field of values of a matrix is the subset of the complex plane defined by

W (J h ) = x H J h x x H x , x ∈ C, x = 0 .
It includes the eigenvalues and is a convex set. For a non-symmetric matrix, the con-vergence of GMRES is better described by the field of values than by the eigenvalues (see [START_REF] Benzi | Field-of-values convergence analysis of augmented Lagrangian preconditioners for the linearized Navier-Stokes problem[END_REF] or [START_REF] Klawonn | Block triangular preconditioners for nonsymmetric saddle point problems: field-of-values analysis[END_REF]). We conjecture that the field of values of the Schur complement J h can be bounded away from 0, independently of N x . Though we have no proof at the moment, we have the following numerical confirmation on figure 6.11, which shows the field of values and isolines of the -pseudo-spectra of the Jacobian matrix J h for 2 different mesh sizes. One can see that the convex hull of the field of values is indeed approximately independent of the mesh size. The figure was obtained thanks to the Eigtool software [START_REF] Wright | Eigtool[END_REF][START_REF] Wright | Large-scale computation of pseudospectra using ARPACK and Eigs[END_REF]. 6.2. Test 2. Ions exchange in a natural system. Our second test case comes for a field study that includes experimental results (see Valocchi et al [START_REF] Valocchi | Transport of ion-exchanging solutes in groundwater: Chromatographic theory and field simulation[END_REF]). We follow the setup given in Fahs et al [START_REF] Fahs | On the efficiency of the direct substitution approach for reactive transport problems in porous media[END_REF], as this reference includes more details for the numerical simulation. In this test case, four aqueous components (Na + , Ca 2+ , Mg 2+ and Cl -are injected into an homogeneous landfill. During the transport, the aqueous components react with the ion exchange sites of the soil (S). Three reactions of ion exchange occur and lead to three adsorbed species S-Na, S2-Ca, S2-Mg. Chemical reactions, constants of equilibrium, initial and boundary conditions are summarized in Table 6 equal to 0.08m , the time step is equal to 0.11089 h corresponding to CFL coefficients of 1. Figure 6.12 shows the evolution of the concentrations in Ca 2+ and Mg 2+ at the end of the column, as a function of time (we use the same units as Fahs et al [START_REF] Fahs | On the efficiency of the direct substitution approach for reactive transport problems in porous media[END_REF] which are different than those originally used by Valocchi et al [START_REF] Valocchi | Transport of ion-exchanging solutes in groundwater: Chromatographic theory and field simulation[END_REF]). The results are in 6.12. Variation of Ca 2+ and Mg 2+ concentrations as function of time in the output of the domain qualitative agreement with those in both references (Figure 1 in Fahs et al [START_REF] Fahs | On the efficiency of the direct substitution approach for reactive transport problems in porous media[END_REF], Figure 11 in Valocchi et al [START_REF] Valocchi | Transport of ion-exchanging solutes in groundwater: Chromatographic theory and field simulation[END_REF]), though it appears difficult to make a precise comparison as the curves are in logarithmic scale on both axes.

Finally, figure 6.13 compares the number of linear and non-linear iterations for the various methods presented. In this case, the number of non-linear iterations was very close for all methods (with Block-Jacobi preconditioning a close winner), and both the Gauss-Seidel preconditioning and the h-method again giving a convergence independent of the mesh size An alternative formulation and block preconditioners for linear system were used to accelerate the convergence of the Krylov method and to reduce CPU time. The results show that the alternative formulation requires less CPU time than other preconditioners, and the number of linear and non linear iterations becomes almost independent of the mesh.

The reactive transport benchmark 1D problem proposed by GNR MoMaS was used to demonstrate the efficiency of the method.

Natural extensions of this work to multidimensional situations are under way, as well as extensions to handle kinetic reactions.

Fig. 6 . 1 .

 61 Fig. 6.1. Geometry of the domain Medium A Medium B Porosity φ (-) 0.25 0.5 Total immobile concentration T S 1 10

Fig. 6 . 2 .

 62 Fig. 6.2. Concentrations profiles. Left figure: C 1 at t = 10, middle figure S at t = 10, right figure C 2 at t = 5010.

Fig. 6 . 3 .

 63 Fig. 6.3. Elution curves (concentrations at right end of the domain as a function of time). Left image: species X3, middle image: species C5, right image: species C2

Fig. 6 . 4 .

 64 Fig. 6.4. Elution curves of X 3 (left) and C 5 (right) concentration, for various mesh resolutions

Fig. 6 . 5 .

 65 Fig. 6.5. Concentration of species C 1 at t = 10 for the different discretization schemes and various mesh resolutions

Figure 6 .

 6 Figure 6.6 compares the concentration of species S (a sorbed species) on the interval bracketing the location of the peak [0, 0.15], at time t = 10, as computed by the three discretization schemes on increasingly refined meshes. Here also, the three schemes give identical solutions when the mesh is fine enough. This time, one needs at least 768 mesh points to obtain a converged solution. One should compare Figure6in[START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF], and also refer to Table4there, where the location and amplitude of the peak are tabulated for all the methods used in the benchmark. We have obtained values of x = 0.0175 for the location of the peak, and S = 0.985 for its height. These values are in the range reported by the other teams, but are different form the "mean" values as reported in[START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF]. They are however very close to the "reference" values found in[START_REF] Carrayrou | Looking for some reference solutions for the reactive transport benchmark of MoMaS with SPECY[END_REF].

Fig. 6 . 6 .

 66 Fig. 6.6. Concentration of species S at t = 10 for the different discretization schemes and various mesh resolutions

6. 1 . 4 .Fig. 6 . 7 .

 1467 Fig. 6.7. Effect of the CFL condition on the concentrations in S (top figure 768 nodes, bottom figure 1536 nodes)

Fig. 6 . 9 .

 69 Fig. 6.9. Non-linear (left) and linear (right) iterations for different preconditioning strategies

Fig. 6 . 10 .

 610 Fig. 6.10. CPU time as a function of the mesh size for the various solution methods

Fig. 6 . 11 .

 611 Fig. 6.11. Field of values (dashed line) and pseudo-spectra (color) for matrix J h , with Nx = 768 and Nx = 1536

Fig. 6 . 13 .

 613 Fig. 6.13. Number of non-linear (left) and linear (right) iterations for ion exchange example

Table 6 . 1

 61 Physical properties of the materials

Table 6

 6

	Total	T 1 T 2 T 3 T 4
	Conc.	0	-2	0	2
	Injection 0.3 0.3 0.3 0
	t∈[0,5000]
	Leaching	0	-2	0	2
	t∈[5000,..]
	.2
	Morel tableau for Chemical equilibrium

Table 6 . 3

 63 Boundary conditions 6.1.1. Sample results. The results obtained in

Table 6 .

 6 4 CFL values for variable time step simulation

	Start time	0	20	100 2500 3200 5000 5100
	End time 20 100 2500 3200 5000 5100 6000
	CFL value	1	5	10	5	40	1	40

Table 6 . 5

 65 .5. Na + Ca 2+ Mg 2+ Cl - Chemical reactions, initial and boundary conditionsIn this test, we consider a column of length L=16 m. The values of the transport parameters are given in Table6.[START_REF] Benzi | Field-of-values convergence analysis of augmented Lagrangian preconditioners for the linearized Navier-Stokes problem[END_REF] The simulated time is T=5000 h. The mesh size is

	S	log K

Table 6 . 6

 66 Transport parameters

* This work was partially supported by GNR MoMaS, CNRS-2439 (PACEN/CNRS, ANDRA, BRGM, CEA, EDF, IRSN), and by the Hydrinv(EuroMediterrannean 3+3) project