Amir Laila
email: laila.amir@inria.fr

Michel Kern
email: michel.kern@inria.fr

Laila Amir

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

Keywords: Geochemistry, transport in porous media, Newton-Krylov methods, advection-diffusionreaction equations Mathematics Subject Classification (2000) 76V05, 65M99

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

Introduction

The simulation of multi-species reacting systems in porous media is of importance in several different fields: for computing the near field in nuclear waste simulations, in the treatment of bio-remediation, in CO2 sequestration simulations and in the evaluation of underground water quality. This work deals with numerical methods for solving coupled transport and chemistry problems. The transport of solutes in porous media is described by partial differential equations of advection-diffusion type, wheres multi-species chemistry involves the solution of ordinary differential equations (if the reactions are kinetic) or nonlinear algebraic equations (if local equilibrium is assumed). After discretization, one is led to a system of nonlinear equations, coupled the unknowns for all chemical species at all grid points.

After the influential paper by Yeh and Tripathi [START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF], operator splitting methods, where transport and chemistry are solved for separately at each time step (possibly iterating to convergence), became the methods of choice. Some representative papers where operator splitting methods are used are [START_REF] Saaltink | On the behavior of approaches to simulate reactive transport[END_REF], [START_REF] Saaltink | A comparison of two approaches for reactive transport modelling[END_REF], [START_REF] Carrayrou | Operator-splitting procedures for reactive transport and comparison of mass balance errors[END_REF], [START_REF] Kanney | Convergence of iterative split operator approaches for approximating nonlinear reactive transport problems[END_REF], [START_REF] Lucille | Chemtrap: a hydrogeochemical model for reactive transport in porous media[END_REF], [41]. Operator splitting methods are easy to implement, and the splitting errors can be controlled by carefully restricting the time step. On the other hand, the time-step restriction can become their main drawback, as it can be difficult to get the fixed point iteration to converge for more difficult problems.

More recently, global methods have become popular, due to the increase in computing power now available. In this approach, the full non-linear system is solved in one step, usually by some form of Newton's method. Most papers use the Direct Substitution Approach (see [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physicsbased preconditioning to biogeochemical transport[END_REF], [START_REF] Fahs | On the efficiency of the direct substitution approach for reactive transport problems in porous media[END_REF]), where one substitutes the chemical equations in the transport equations. On the other hand, the problem can also be put in the form of a Differential Algebraic Equations (DAE), enabling the use of powerful software (see [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF]). Finally, the chemical equations can be eliminated locally, and a system involving transport equations, with a source term coming from the reactions has to be solved. This approach is taken in [START_REF] Kraeutle | A new numerical reduction scheme for fully coupled multicomponent transport-reaction problems in porous media[END_REF][START_REF] Kraeutle | A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions[END_REF], where additionally a reduction method leads to a smaller system. Most of the papers quoted above employ a Newton method for solving the nonlinear system at each time step, with the difficulty that the Jacobian matrix has to be computed, stored and factored. This can become problematic for large problems, and Hammond et al. [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physicsbased preconditioning to biogeochemical transport[END_REF] have used the Jacobian-Free Newton-Krylov method, where the Newton correction is solved for by an iterative method. The Jacobian is only needed through the computation of a directional derivative. The method keeps the fast convergence of Newton's method, while only requiring Jacobian matrix-vector products, and these can be approximated by finite differences.

The method presented in this paper is a global method where the chemical equations are eliminated locally, leading to a nonlinear system where the transport and chemistry subsystems remain separated. Thus the residual can be evaluated by calling separately written transport and chemistry modules. The system is then solved by a Newton-Krylov method, and it will be shown how the Jacobian matrix-vector product can also be computed by the same module. Thus the main contribution of this paper is to show that a global method can be implemented while still keeping transport and chemistry modules separated. This property will be referred to as using "blackbox solvers". As the chemical equilibrium equations are not substituted in the transport equations, the transport and chemistry parts of the nonlinear residual are easily identified, and can each be computed by calling on standard solution modules.

An outline of the paper is as follows. In section 2 the chosen model is explained, and the methods used for solving the (non-reactive) transport part, and the chemical equilibrium system are detailed. Section 2.3 shows how we obtain the coupled model. Couple formulations and coupling algorithms are the subject of section 3, beginning with a review of existing methods, while our approach is presented in section 3.2. Numerical results, in particular experience with the MoMaS benchmark, are shown in section 4.

Reactive transport equations

In this work, the transport of several reacting species in a single phase flow through a porous medium is considered. The species can react both between themselves and with the porous matrix. In this section the numerical methods used to solve the individual subsystems of the coupled problem will be described.

Transport model

The transport of a single species through a porous medium (a domain Ω ⊂ R d , with d = 1, 2 or 3), with porosity φ, in a known Darcy field u, subject to dispersion and molecular diffusion, follows the linear advectiondispersion equation

φ ∂c ∂t + L(c) = q, in Ω (1)
where

L(c) = ∇ • (u c) -∇ • (D∇c),
is the transport operator, and q is a source term. The diffusion-dispersion tensor D is given by

D = d e I + |u| (α L E(u) + α t (I -E(u))) , E i j (u) = u i u j |u| 2 ,
where d e is the molecular diffusion coefficient, and α L (resp. α t ) is the longitudinal (resp. transverse) dispersivity coefficient.

In this work, we restrict to a one dimensional problem, so that the transport equation over a bounded interval Ω =]0, L[ can be written as

φ ∂c ∂t + ∂ ∂x -D ∂c ∂x + uc = q, 0 < x < L, 0 < t < T, (2)
where the porosity φ and the diffusion-dispersion coefficient D can both depend on space. Because the flow is assumed compressible, the velocity u is taken to be a constant. 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). More general boundary conditions could easily be accommodated.

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 a cell-centered finite volume scheme will be used, see for instance [START_REF] Eymard | Finite volume methods[END_REF]. The interval

[0, L] is divided into N g intervals [x i-1 2 , x i+ 1 2 ] of length h i , where x 1 2 = 0, x N g + 1 2
= L. For i = 1, .., N g , denote by x i the center and x i+1/2 the right end of element i. Finally, denote by c i , i = 1, .., N g the approximate solution in cell i.

Equation ( 2) is written in the form

φ ∂c ∂t + ∂ϕ ∂x = q, ( 3
)
where the flux ϕ(x,t) = -D ∂c ∂x + uc has been split as the sum of a diffusive flux ϕ d = -D ∂c ∂x and an advective flux ϕ a = uc.

Equation ( 3) is integrated over a cell ]x i-1/2 , x i+1/2 [ giving

φ i h i dc i dt + ϕ d,i+ 1 2 + ϕ a,i+ 1 2 -ϕ d,i-1 2 -ϕ a,i-1 2 = h i q i , i = 2, . . ., N g . (4)
The flux approximations required to close the system are provided by finite differences. The diffusive flux needs a value for the diffusion coefficient, which is taken as the harmonic average (as done in mixed finite element methods):

ϕ d,i+ 1 2 = -D i+ 1 2 c i+1 -c i h i+ 1 2 (5)
with

D i+ 1 2 = 2D i D i+1 D i + D i+1 , D 1 2 = D 1 , D N g + 1 2 = D N g and h i+ 1 2 = h i + h i+1 2
For the advective flux, an upwind approximation is used, so that (assuming u > 0), ϕ 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. The semi-discrete system can be summarized by the finite dimensional system

M dc dt + Lc = q + g, ( 6
)
where c ∈ R N g now represents the vector of cell concentrations, L ∈ R N g ,N g is the matrix form of the transport operator, M ∈ R N g ,N g is a mass matrix accounting for variable porosity and mesh size, q ∈ R N g is a give source term and g ∈ R N g represents the effects of the boundary conditions.

Time discretization

Let denote by ∆t the time step (taken constant for simplicity) used to discretize the time interval [0, T ] and denote by c n i the (approximate) value of c i (n∆t). The first and most straightforward alternative is to discretize equation ( 6) by the backward Euler method, see for instance [START_REF] Ascher | Numerical Methods for Evolutionary Differential Equations[END_REF]. This is the method that is used in section 3 to keep the description simple, but is not the recommended method, as it leads to an overly diffusive scheme.

Better alternatives are obtained by exploiting the structure of the transport operator, and by using different time discretizations for the advective and for the diffusive parts. Specifically, the diffusive terms should be treated implicitly, and the advective terms are better handled explicitly.

If this idea is applied directly to equation ( 6), the resulting fully discrete scheme is only stable under a CFL (Courant-Friedrichs-Lewy) condition 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), an alternative is to use an operator splitting scheme, as 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 (editors: Alain Bourgeat and Michel Kern)[END_REF][START_REF] Mazzia | A time-splitting technique for the advection-dispersion equation in groundwater[END_REF]). In this work, splitting is used only within the (linear) transport step, but recent papers by Kačur et al. [START_REF] Frolkovič | Semi-analytica solution of a contaminant transport equation with nonlinear sorption in 1D[END_REF][START_REF] Kačur | Solution of contaminant transport with equilibrium and nonequilibrium adsorption[END_REF] apply splitting directly to a transport with sorption model by solving (analytically) a nonlinear advection step, followed by a nonlinear diffusion step. This is different from operator splitting as used in geochemical models, as the chemistry terms are solved for together with the transport terms.

The splitting scheme works by taking several small time steps of advection, controlled by a CFL condition, within a large time step of diffusion. The scheme has been shown to be unconditionally stable, and has a good behavior in advection dominated situations.

More precisely, the time step ∆t will be used as the diffusion time step, it is divided into M time steps of advection ∆t c such that ∆t = M∆t c where M ¿1, the advection time step will be controlled by CFL condition.

Advection step

The interval [t n ,t n+1 ] is divided into M intervals [t n,m ,t n,m+1 ], m = 0, ...M -1, where t n,0 = t n , t n,M = t n+1 . Denote c n,m i the approximate concentration c at time t n,m and c n,0 = c n . The advection equation is discretized in time using the explicit Euler method to obtain

     φ i ( c n,m+1 i -c n,m i ∆t c ) + u c n,m i -c n,m i-1 h i-1/2 = 0, i = 2, . . . , N g , c n,m+1 1 = c g (t n,m+1 ) m = 0, . . . , M -1. ( 7
)
Diffusion step The diffusion part is discretized by an implicit Euler scheme, starting from c n,M i :

- D i-1 2 h i-1 2 ∆tc n+1 i-1 + φ i h i + D i+ 1 2 h i+ 1 2 ∆t + D i-1 2 h i-1 2 ∆t c n+1 i - D i+ 1 2 h i+ 1 2 ∆tc n+1 i+1 = φ i h i c n,M i + q i h i ∆t, i = 2, . . . , N g -1 (8)
As above, 2 equations accounting for the boundary conditions must be added.

Chemical equations

The chemical model is described in this section. In this study, we assume a local chemical equilibrium at every point, which means that the chemical phenomena occur on much faster scale than transport phenomena. This is a common modeling assumption for reactive transport in porous media, at least when the only reactions considered are aqueous phase and sorption reactions (these are "sufficiently fast" reactions according to Rubin [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 would not be the case if mineral dissolution was taken into account, as these reactions typically need kinetic models.

Consider a set of N e chemical species (X j ) j=1,...,N e linked by N r reactions

N e ∑ j=1 ν i j X j ⇆ 0, i = 1, . . ., N r
where ν is the stoichiometric matrix. Following Morel [START_REF] Morel | Principles and Applications of Aquatic Chemistry[END_REF], we distinguish between component and secondary species by extracting a full rank matrix from ν. Component species are a minimal subset of the species such that the other secondary species can be written in terms of them (in a unique way). Each secondary species gives rise to a reaction that expresses how it is formed in terms of the components, and to a mass action law that gives the value of its activity in terms of the component activities. Similarly, each component gives rise to a conservation equation, expressing how the given total concentration of such a component is distributed among the component itself and the secondary species. Additionally, in the context of reactive transport, it is required to know how the species are split between those that are in solution, and those that have been adsorbed on the solid matrix (in this paper we do not take precipitation into account). We thus introduce (with obviously N e = N c + N s + N x + N y ) -mobile components c j , j = 1, . . . , N c , -fixed components s j , j = 1, . . . , N s , -mobile secondary species x j , i = 1, . . ., N x , -fixed secondary species y j . i = 1, . . . , N y .

We have identified the name of the species with their concentrations, and we assume an ideal solution (activities and concentrations are identified). Mobile secondary species x can be expressed as linear combinations of mobile components while secondary fixed species depend on both mobile and fixed components. Therefore the mass action laws are written as

x i = K xi N c ∏ j=1 c S i j j , i = 1, . . . , N x , y i = K y i N c ∏ j=1 c A i j j N s ∏ j=1 s B i j j , i = 1, . . ., N y , (9)
where K xi and K y i are the equilibrium constants, and S i j , A i j and B i j are the entries of the stoichiometric matrices

S ∈ R N c ×N x , A ∈ R N c ×N y and B ∈ R N s ×N y .
Mass conservation for each component is expressed in the form

c + S T x + A T y = T, s + B T y = W, ( 10
)
where T j is the total concentration of the mobile component j, and W j is the total concentration of the fixed component j (T and W are vectors of size N c and N s respectively). In the case of ion exchange, the second mass conservation equation is simply B T y = W , and W is the Cationic Exchange Capacity of the porous matrix (see Appelo and Postma [START_REF] Appelo | Geochemistry, Groundwater and Pollution[END_REF]). As will be seen later, in the context of coupled transport and chemistry, T j is given by the transport model and W is constant. In a closed chemical system, T j would be part of the data (total concentration of the components). Due to the wildly different orders of magnitude of the concentrations that are commonly encountered, the chemical problem is reformulated by using as main unknowns the logarithms of the concentrations. This has the added advantage that concentrations are automatically positive, and has become the standard way to solve the problem [28]. An additional advantage has been pointed out by Samper et al. [41]: by taking the logarithms of the concentrations as unknowns, the Jacobian of the nonlinear system is symmetric, and with a proper choice of the component species, it can be shown to be diagonally dominant, and thus nonsingular. The symmetry can also be seen on equation ( 16) below. Let log u be the vector with entries log u i , where u i are the entries of vector u. Equations ( 9) can then be rewritten as a linear system

log x = log K x + S log c log y = log K y + A log c + B log s (11)
The nonlinear system of equations ( 10) and ( 11) forms what will be called the chemical problem. In the sequel, it will be assumed that this problem always has a (positive) solution (c, s), for all feasible values of the data T and W . This is true in our simplified settings because the chemical equilibrium problem is a consequence of the minimization of the Gibbs free energy, which can be shown to be convex in the absence of minerals (see [START_REF] Shapiro | Mass action laws and the Gibbs free energy function[END_REF]).

To solve the chemical problem, a variant of Newton's method is used. As is well known, Newton's method is not always convergent, unless the initial point is sufficiently close to the solution. However, and this is especially true in the context of a coupled code where the chemical problem will be solved repeatedly, it is essential to ensure that the solver "never" fails. 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. In order to get a smaller system, the secondary concentrations are eliminated, and the system to be solved involves only

lc = log c ∈ R N c and ls = log s ∈ R N s . Define the function H : R N c +N s → R N c +N s by H lc ls = exp(lc) + S T exp(log K x + Slc) + A T exp(log K y + Alc + Bls) exp(ls) + B T exp(log K y + Alc + Bls), ( 12
)
where the notation exp(v) for a vector v means the vector with elements exp(v j ), then equations ( 10) and ( 11) are equivalent to:

H lc ls = T W . ( 13
)
This is the nonlinear system that to be solved for lc and ls, given T and W . The secondary concentrations can then be computed from equation [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF].

When solving the coupled problem, the distribution of the species between their mobile form and their fixed form will be needed. The individual concentrations must still be solved for, but they are intermediate quantities. Once the component concentrations have been computed as described in the previous paragraph, one can compute for each species its mobile part C j and its fixed part F j by

C = c + S T x, F = A T y. ( 14
)
Note that, by definition, the relationship T = C + F holds.

In the formulation to be presented below, it will be convenient to represent the mapping from the vector of total concentrations to the vector of fixed concentrations. This mapping, denoted by Ψ, is defined by first solving the chemical problem [START_REF] Eymard | Finite volume methods[END_REF], then computing F by [START_REF] Fahs | On the efficiency of the direct substitution approach for reactive transport problems in porous media[END_REF]. More precisely

ψ : R Nc → R Nc T → ψ(T ) = A T y, ( 15
)
where equation ( 13) is first solved for lc and ls, then y is computed by [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF].

It is important to keep in mind that computing Ψ(T ) means solving the chemical system (plus some simple computations), as this will be the most expensive part when evaluating the residual of the coupled system (see eq.( 26) in section 3.2). As this will be useful later on, the computation of the Jacobian of Ψ is outlined here. Assume Ψ(T ) itself has been computed, so that the nonlinear system (13) has been solved. First, the Jacobian matrix of H should also be computed as part of the solution process. This is almost certainly needed for solving the chemical problem, if Newton's method is used. Differentiating equation ( 12) leads to:

H ′ lc ls = diag(exp(lc)) 0 0 diag(exp(ls)) + S T A T 0 B T diag(x) 0 0 diag(y) S 0 A B , (16)
where diag(v) is the diagonal matrix with vector v along the diagonal. Then, by an application of the implicit function theorem (see for instance [START_REF] Rudin | Principles of mathematical analysis, 3rd edn[END_REF]), and by differentiating equation [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF], there comes

ψ ′ (T ) = A T diag(y) A B H ′ lc ls -1 I 0 . ( 17
)
It should be stressed that the Jacobian of H is needed to computed the Jacobian of Ψ (inverting it is straightforward, as this will usually be a small matrix). This may prove problematic in practice for several reasons. First, the chemical solver may not give access to the Jacobian, even if it is used internally. This is a limitation to the "black-box" approach. Second, for more realistic chemical models, including non-ideal chemistry, and taking minerals into account, computing the Jacobian may be much more difficult than the fairly simple computation outlined above. As a last resort, one could compute the Jacobian by finite differences, but it will be argued in section 3.2 that, for this particular problem, the analytical computation is more efficient.

Coupled transport and chemistry

The starting point for the coupled model is the following set of equations for the total, mobile and fixed concentrations of each component

       φ ∂C j ∂t + φ ∂F j ∂t + L(C j ) = 0 j = 1, . . . , N c ∂W j ∂t = 0, j = 1, . . . , N s (18)
These equations can be derived from the individual conservation equations by standard algebraic manipulations, see for instance Yeh and Tripathi [START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF]. It is the formulation given in the benchmark definition [START_REF] Carrayrou | Presentation of the MoMaS reactive transport benchmark[END_REF], see also [START_REF] Saaltink | A mathematical formulation for reactive transport that eliminates mineral concentrations[END_REF], [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF]. The second equation is obvious, as W j was taken as a constant (at each point in space).

Taking into account the relation T j = C j + F j , j = 1, . . . , N c noted above, the first equation of the system is equivalent to

φ ∂T j ∂t + L(C j ) = 0 j = 1, . . . , N c , (19)
where T j is the total concentration, C j the total mobile concentration, and F j the total fixed concentration for component j.

From now on, L will denote the discretized transport operator, as defined in equation [START_REF] Carrayrou | Presentation of the MoMaS reactive transport benchmark[END_REF]. Each unknown concentration depends on both the grid point index, and the chemical species index. We will use a notation inspired from Matlab. For a concentration u i j , where i ∈ [1, N g ] represents the spatial index and j ∈ [1, N c ] represents the chemical index, we shall denote by -u :, j the column vector of concentrations of species j at all grid points; -u i,: the row vector of concentrations of all chemical species in grid cell x i . The unknowns will be numbered first by chemical species, then by grid points. Thus all the unknowns for a single grid point are numbered contiguously.

The coupled problem is obtained by putting together equation ( 19) above with the definition of the chemical solution operator ψ, defined in eq. ( 15) (the subscript T denotes transposition):

         M ∂C :, j ∂t + M ∂F :, j ∂t + L(C :, j ) = g :, j j = 1, . . . , N c T i j = C i j + F i j , i = 1, . . . , N g , j = 1, . . . , N c F i,: = ψ(T T i,: ) T , i = 1, . . . , N g ( 20
)
This system is then discretized in time to obtain the fully discrete coupled nonlinear system. In this work we restrict to a simple backward Euler scheme with constant step-size, noting that other more sophisticated, strategies are obviously possible (in particular, an adaptive step-size is essential for efficiency). Denoting time indexes by a superscript, the following system is obtained

           M C n+1 :, j -C n :, j ∆t + M F n+1 :, j -F n :, j ∆t + L(C n+1 :, j ) = g :, j j = 1, . . . , N c T n+1 i j = C n+1 i j + F n+1 i j i = 1, . . . , N g , j = 1, . . . , N c F n+1 i,: = ψ((T n+1 i,: ) T ) T i = 1, . . . , N g (21)
This is the system to be solved at each time step.

Formulation and coupling algorithms

The formulation of reactive transport seen above gives rise to a large system of nonlinear equations. For complex problems, its solution will require a large amount of computer time, which makes it important to choose an appropriate method. In this section, several formulations and approaches that have appeared in the literature will be reviewed. Thanks to the relationship T = C +F, it is easy to eliminate one of the 3 variables, and this leads to different formulations for the coupled problem, depending on which variables are kept in the transport equation. We keep the system continuous in time, as it makes the notation somewhat lighter, but the same manipulations can obviously be done at the discrete level too.

According to Saaltink et al. [START_REF] Saaltink | A mathematical formulation for reactive transport that eliminates mineral concentrations[END_REF], see also Salignac [START_REF] Salignac | Transport multi-espèces et réactions géochimiques en aquifère : développement et validation du modèle couplé HYTEC 2D[END_REF], one can derive three main formulations from the system given in [START_REF] Hoteit | Nuclear waste disposal simulations: Couplex test cases: Simulation of transport around a nuclear waste disposal site: The COUPLEX test cases (editors: Alain Bourgeat and Michel Kern)[END_REF]:

-formulation (TC) where T is the principal variable , C the transported variable

M ∂T : j ∂t + L(C : j ) = g :, j (22)
This is the formulation used by Erhel et al. in [START_REF] De Dieuleveult | A numerical model for coupling chemistry and transport[END_REF][START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF], as it lends itself best to a DAE type algorithm. It is not convenient for our purpose, as the transport equation then involves both T and C, and is thus not easily used with an existing transport solver. -formulation (TT) where T is the principal variable , T transported variable

M ∂T : j ∂t + L(T : j ) + L(F : j ) = g :, j (23)
This seems to be the least satisfactory formulation, as the transport operates on the fixed species, and for this reason it will not be considered further. -formulation (CC) where C is the principal variable , C transported variable

M ∂C : j ∂t + M ∂F : j ∂t + L(C : j ) = g :, j (24)
This is formulation 4 in Saaltink et al. [START_REF] Saaltink | A mathematical formulation for reactive transport that eliminates mineral concentrations[END_REF], and is the formulation chosen below. It has been reported that this formulation is the least suitable for use in an operator split algorithm, because C and F are used at different time levels (to compute the data for the chemical problem). When this formulation is used in a global method this should not matter as much, as the iterations are ran to convergence, and both values should eventually get close to their limits.

Formulation (CC) will be used in the rest of the paper, because it takes the form of a standard transport operator, with a source term coming from the chemical part. Its structure is closely related to the system describing single species transport with sorption, as seen for instance in [START_REF] Frolkovič | Semi-analytica solution of a contaminant transport equation with nonlinear sorption in 1D[END_REF], or [START_REF] Kanney | Convergence of iterative split operator approaches for approximating nonlinear reactive transport problems[END_REF], with the main differences that the unknown is a vector of concentration, and mostly that what plays the role of the sorption isotherm is the implicitly defined function Ψ introduced in (15).

Review of former approaches

At each time step, the system given by (21) (one transport equation for each component and one chemical system for each grid point) forms a large nonlinear system, whose size is the number of components times the number of grid points. This system has traditionally been solved by a sequential two-steps approach, as reviewed below (cf. [START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF]). However, this method suffers from several defects: it may severely restrict the step size to ensure convergence, and if used non-iteratively it is only first order in time, which may introduce additional errors (cf. [START_REF] Carrayrou | Operator-splitting procedures for reactive transport and comparison of mass balance errors[END_REF]). Due to its quadratic convergence rate, Newton's method would be an ideal candidate for solving the system. On the other hand, a practical difficulty has to be reckoned with: Newton's method requires the solution of a linear system with the Jacobian matrix at each iteration step. In realistic situations, it will not be possible to store, much less factor, the Jacobian matrix. As will be seen in section 3.2, this difficulty can be overcome by resorting to an iterative method for solving the linear system.

Sequential approach

The sequential approach consists of separately solving the chemical equations and the transport equations. The method has been used in numerous papers: see for instance [START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF], and also [START_REF] Kanney | Convergence of iterative split operator approaches for approximating nonlinear reactive transport problems[END_REF], [START_REF] Lucille | Chemtrap: a hydrogeochemical model for reactive transport in porous media[END_REF], [START_REF] Saaltink | A mathematical formulation for reactive transport that eliminates mineral concentrations[END_REF], [START_REF] Carrayrou | Operator-splitting procedures for reactive transport and comparison of mass balance errors[END_REF] or [START_REF] Mayer | Formulation of the multicomponent reactive transport code MIN3P and implementation of MoMaS benchmark problems[END_REF]¡. At each iteration, a transport equation for each component is solved first, with a source term given by the (change in) fixed concentration at the previous iteration. This total mobile concentrations will be added to a total fixed concentration computed in the previous iteration, to obtain the total used as data for solving a chemical problem at each grid point. These steps are then iterated until convergence.

In the geochemical literature, this is known as an operator splitting approach (usually called Standard Iterative Approach, or SIA), but it is more properly a block Gauss-Seidel methods on the coupled system, as each subsystem is solved alternatively. The method is quite appealing, as it is easy to implement starting from separate transport and chemistry codes, and can provide good accuracy if implemented carefully, as shown in the references above). As will be seen below, theses advantages can be retained in the Newton-Krylov framework.

The Standard Non-Iterative Approach (SNIA) is the case where only one iteration of the method is carried out at each time step. In that case, splitting errors can become important, and the method is not really suitable for difficult problems.

The SIA approach does not suffer from splitting errors if the tolerance is small enough, but it may require a small time step to obtain convergence in the case of stiff problems. The main drawback of the method is thus that the size of the time step is used to control convergence, and not based on the physical character of the solution.

Direct Substitution Approach

As computing power increased, it was recognized that the operator splitting methods of the previous sections could not satisfactorily handle difficult problems, and more tightly coupled method came to more widespread use.

The Direct Substitution Approach method consists in solving for the individual concentrations of the components, that is substituting equations ( 10)- [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF] in equation (1) (this can be done explicitly, as in Hammond et al. [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physicsbased preconditioning to biogeochemical transport[END_REF], or implicitly, as in Kräutle et al. [START_REF] Kraeutle | A new numerical reduction scheme for fully coupled multicomponent transport-reaction problems in porous media[END_REF][START_REF] Kraeutle | A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions[END_REF], or Saaltink et al. [START_REF] Saaltink | A mathematical formulation for reactive transport that eliminates mineral concentrations[END_REF]). It is also possible to reformulate the problem as a differential algebraic system (DAE), and to take advantage of the high quality software available for such problems, as in Erhel et al. [START_REF] De Dieuleveult | A global strategy for solving reactive transport equations[END_REF], [START_REF] De Dieuleveult | A global approach for reactive transport: application to the benchmark easy test case of MoMaS[END_REF] or [START_REF] De Dieuleveult | Un modèle numérique global et performant pour le couplage géochimie-transport[END_REF] . A high performance parallel implementation is described by Hammond et al. [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physicsbased preconditioning to biogeochemical transport[END_REF], using a Jacobian-Free Newton-Krylov method (see section 3.2).

The main advantages of this approach are to avoid the errors caused by the separation of operators, and to allow fast convergence independently of the time step, but its principal drawback is the need to form and to store the Jacobian matrix especially for a large problem. Moreover, sometimes it may be difficult to calculate the exact derivatives for geochemical processes especially when precipitation phenomena or kinetic reactions are taken into account.

The size of the system can be made smaller by means of a reduction method, cf Kräutle et al. [START_REF] Kraeutle | A new numerical reduction scheme for fully coupled multicomponent transport-reaction problems in porous media[END_REF][START_REF] Kraeutle | A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions[END_REF], and [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]. The reduction method makes a change of variables in the chemical system, so that a set of decoupled transport equation is first solved, leaving a smaller nonlinear system, that is still solved with Newton's method.

A Newton-Krylov based fully coupled method

As was already mentioned in the previous section, Hammond et al. [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physicsbased preconditioning to biogeochemical transport[END_REF] have used a Newton-Krylov method for solving the system obtained from the DSA approach. Substituting the chemical equations in the transport operator is the most straightforward way of formulating the coupled problem, but leads to a system where chemistry and transport terms are mixed, and makes it virtually impossible to separate the transport and chemistry modules. However, this separation is seen as one of the important advantages of the operator splitting approaches.

By coupling the formulation given in section 2.3 with the Newton-Krylov framework, a strongly coupled method that can be implemented by keeping transport and chemistry separate is obtained. Thus, the chemical equations are not directly substituted in the transport equation, but the function Ψ introduced previously in ( 15) is used to represent the effect of chemistry. Different formulations could be adopted depending on the choice of unknowns (refer back to section 3). In this work, both the total mobile and fixed concentrations, and also the total concentrations (though they could easily be eliminated) are chosen as main unknowns.

Even though this method may be more expensive than the methods based on DSA, its main advantage is to make it possible to treat chemistry as a black-box, even in the Newton-Krylov context. This may be important, as chemical simulators are becoming increasingly sophisticated.

Recall (equation ( 21)) that the nonlinear system to be solved at each time step is

         (M + ∆tL)C n+1 :, j + MF n+1 :, j -b n :, j = 0, j = 1, . . . , N c T n+1 -C n+1 -F n+1 = 0, F n+1 i,: -ψ (T n+1 i,: ) T T = 0, i = 1, . . . , N g ( 25
)
where b n :, j = MC n :, j + ∆t g n+1 :, j + M∆tF n :, j is known. Denoting by

G : R 3N c N g → R 3N C N g the function G   C T F   =     (M + ∆tL)C :, j + MF :, j -b n :, j j∈[1,N c ] T -C -F F i,: -ψ (T i,: ) T T i∈[1,N g ]     , (26)
the nonlinear problem to be solved at each time step is G(Z) = 0, where Z denotes the vector

  C n+1 T n+1 F n+1   .
Recall that at each step of the "pure" form of Newton's method for solving G(Z) = 0, one should compute the Jacobian matrix J = G ′ (Z k ), solve the linear system (usually by Gaussian elimination)

J δZ = -G(Z k ) (27)
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.

The main drawback of the method for large scale problems is again the need to form, and then factor, the Jacobian matrix. For coupled problem such as the one studied in this paper, there is the additional difficulty of simply computing the Jacobian: the numerical methods for transport and chemistry are quite different, and it is even possible that the simulation codes have been written by different groups.

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 physicsbased preconditioning to biogeochemical transport[END_REF], to which our work is closely related), is a variant of Newton's method where the linear system that arises at each step of Newton's method is solved by an iterative method (of Krylov type). The main advantage of this type of method is that the full Jacobian is not needed, one just needs to be able to compute the product of the Jacobian with a vector. As this is a directional derivative, this leads to the Jacobian free methods, where this product is approximated by finite differences. However, for some problems, it may be possible to compute the needed directional derivative exactly. As will be seen below, this is the case for our coupled problem, provided the Jacobian of the chemical problem can be computed. This is both cheaper and more accurate.

The main contribution of this paper is to show that the formulation given above lends itself to an implementation of Newton's method that allows to keep the two codes separate. This is in keeping with thte philosophy set forth in the review paper by Keyes and Knoll [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF] that a Newton-Krylov solver can often be made by wrapping a classical split-step solver. This is what is being done here, as the formulation to which the Newton-Krylov method is applied is the one used for operator splitting. Additionally, it will be shown below that the Jacobian may even be formed in block form, provided the individual codes provide their Jacobians (this is obviously easier for transport than for chemistry), and this obviously carries over to the Jacobian-vector product.

At this point, it is appropriate to add a few comments on the size of the problems envisioned. The examples used in this work are small scale, one dimensional, problems. They can hardly be called large. On the other hand, we believe they are representative of the problems that will be encountered in more realistic applications. For such problems, in 2 or 3 space dimensions, involving tens or hundreds of thousands of grid points and several tens of chemical species, the nonlinear system will indeed be very large, and a method like that of Hammond et al. [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physicsbased preconditioning to biogeochemical transport[END_REF], or like the method presented in this section will be necessary.

A Krylov subspace method (see for instance [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF]) is used to approximately solve the linear system in equation [START_REF] Lagneau | HYTEC results of the MoMas reactive transport benchmark[END_REF]. The linear iterates are drawn from the Krylov subspace, K j = span{r 0 , Jr 0 , J 2 r 0 , . . ., J j-1 r 0 }. In the GMRES method (see [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]), the iterates are defined to minimize the residual JδZ j + G(Z) 2 over K j . Other methods, such as Bi-CGSTAB [45] or QMR [START_REF] Freund | QMR: a quasi-minimal residual method for non-Hermitian linear systems[END_REF] could be used as well.

As the linear system is not solved exactly, the convergence theory for Newron's method does not apply directly. However, the theory has been extended by Dembo et al. [START_REF] Dembo | Inexact newton methods[END_REF] to the class of Inexact Newton methods, of which the Newton-Krylov methods are representatives. The main consequences of this analysis are summarized below.

An important issue in such methods is the stopping criterion for the inner linear iteration. A stopping criterion of the form

JδZ + G(Z k ) ≤ η k G(Z k ) (28)
in this context, as the initial iterate is usually 0. The choice of the forcing term η k should strike a balance between two conflicting goals:

-Keep the (local) convergence of Newton's methods; -Avoid over-solving, that is taking too many linear iterations when still far away from the nonlinear solution.

The first goal will tend to require a small value for η k , while the second one obviously tends to make η k larger. It has been shown (see theorem 6.1.4 in [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF]) that provided η k is bounded away from 1, the inexact Newton's method will converge, and that superlinear convergence obtains if η k goes to zero faster than G(Z k ) . Based on this result, the strategy proposed by Kelley in [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF] (after the choice in [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF]) computes η k as

η k = γ G(Z k ) 2 / G(Z k-1 ) 2 , (29)
where γ ∈ [0, 1] is a parameter (the value suggested in [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF] is γ = 0.9). Safeguards are added to this choice in order to prevent η k to become too close to 1, or too small. It is also necessary to globalize the algorithm, and this can be done using a line search, just as in the "classical" Newton's method. The other main practical advantage of the Newton-Krylov methods is that they do not require forming the Jacobian matrix. All that is needed is the ability to compute the product of the Jacobian matrix by an arbitrary vector, in order to enlarge the Krylov subspace. This matrix-vector product can be interpreted as a directional derivative. This means that, for complex functions G it may not be necessary to compute the Jacobian, at the cost of one extra evaluation of the function itself. It turns out, however, that in our case, this trade-off is not advantageous. Indeed, it is well known that the most expensive part of the evaluation of G is the solution of the chemical problem at each grid point. On the other hand, it was shown above that computing the Jacobian of ψ is actually cheaper than computing ψ itself (once ψ has already been computed), as it only involves the solution of a linear system (see equation ( 17)), whereas computing ψ itself requires the solution of a nonlinear system.

It will now be shown how the method can be implemented, given modules for transport and chemistry. The first ingredient needed is the computation of the residual, that is evaluating the function G defined in [START_REF] Kraeutle | A reduction scheme for coupled multicomponent transport-reaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions[END_REF]. Given a vector Z = C T F , Z is first split into its three components, and each sub-vector is regarded as a N g × N c matrix, as in section 2.3. Then G(Z) is computed by block:

-For the transport block, the transport operator is applied to each species C :, j , with a source term given by -M F :, j -F n :, j ∆t , for j = 1, . . . , N c (F n denotes F at the previous time step); -The second block is the trivial computation T -C -F; -The third block is the solution of the chemical problem at each grid point: F :,i -Ψ(T i,: ), for i = 1, . . . , Ng.

This shows that the first block will only need transport related quantities, whereas the third block will only call chemistry related ones. Actually, these are the same computations that would be needed for implementing a operator splitting method.

As far as the Jacobian matrix-vector product is concerned, and using the computation in section 2.2, the action of the Jacobian on a vector v =

v c v T v F
(that is the directional derivative of G in the direction of the vector v) can be computed as

J   v C v T v F   =     (M + ∆tL)v C :, j + Mv F :, j j∈[1,N c ] -v C + v T -v F v F i,: -v T i,: (ψ ′ (T T i,: )) T i∈[1,n x ]     . (30)
Even though it is not used as such in this work, it is valuable to examine the structure of the Jacobian. As the previous computation shows, the Jacobian also has a natural block structure. Recall that the unknowns are numbered by species at each point in space. Then the block corresponding to the action of L can be written using the Kronecker product (see for instance [19]) as A = (M + ∆tL) ⊗ I. Then the Jacobian matrix is

J =   A 0 M -I I -I 0 -ψ ′ (T T ) I   , (31)
where ψ ′ (T ) = diag(ψ ′ (T T 1,: ), . . ., ψ ′ (T T N g ,: ) is the Jacobian of ψ, and for each i = 1, . . . , N g , ψ ′ (T T i,: ) is a small N c by N c block. The structure of the Jacobian is illustrated on figure 1, for the case N g = 10, N c = 3. It is a 3 × 3 block matrix, each bock being of size N g × N c . We can clearly see the different parts of the Jacobian: the transport part in the upper left corner has 3 diagonals corresponding to the Kronecker product structure (remember that L is tridiagonal), and the chemistry part at the bottom has 10 3 × 3 blocks.

It would in principle possible to compute and store the Jacobian matrix according to equation [START_REF] Mazzia | A time-splitting technique for the advection-dispersion equation in groundwater[END_REF] as a sparse matrix, and to compute the matrix-vector product using a general purpose routine. The advantage of the method given in equation [START_REF] Mayer | Formulation of the multicomponent reactive transport code MIN3P and implementation of MoMaS benchmark problems[END_REF] is that the structure of the Jacobian is fully exploited, which leads to a much more economical computation.

Numerical results

Ion exchange

The following example of advective transport in the presence of cation exchanger is adopted as a first test case comparison of both approaches. The example is used in the documentation of PHREEQC-2 [START_REF] Parkhurst | User's guide to PHREEQC (version 2)-A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations[END_REF] as Example 11. The one-dimensional simulation problem describes a column experiment where the chemical composition of the effluent from a column containing a cation exchanger is simulated. Initially, the column contains a sodiumpotassium-nitrate solution in equilibrium with the cation exchanger. The column is then flushed with three Fig. 1: The block structure of the Jacobian matrix pore volumes of calcium chloride solution, so that an equilibrium state with calcium and chloride is reached. Calcium, potassium, and sodium react to equilibrium with the exchanger at all times. The flow and transport parameters used for this example are presented in Table 1, and the initial and injected concentrations are listed in Table 2. The Cationic Exchange Capacity for the exchanger is 1.1 mmol/l.

Na + + X -⇋ NaX K + + X -⇋ KX 1 2 Ca 2+ + X -⇋ 1 2 CaX 2
with NaX, KX and CaX 2 are (sorbed) complexes, and X indicates exchange site with charge -1

Comparison with Phreeqc

Figure 2 shows elution curves, that is the evolution of the concentration of the various species at the end of the column, as a function of time. The sorbed potassium and sodium ions are successively replaced by calcium. Because potassium exchanges more strongly than sodium (as indicated by a larger value of log K in the exchange reaction), sodium is released first, followed by potassium. Finally when all of concentration has been released, the concentration of calcium increases to its steady-state value, the potassium is displaced from the exchanger and the concentration in solution increases to balance the Cl -concentration. Both the sequential method and the global method described in section 3.2 have been applied to the test case described in section 4.1. Both the computational demands and the accuracy of the solutions will be compared. As can be seen on figure 2, the results obtained are close to those computed by PhreeqC. One can still see differences both in the location and amplitude of the peak in potassium concnetration, and in the region where the three curves cross. These results are also comparable to those obtained by Xu et al. [START_REF] Xu | Modeling of non-isothermal multi-component reactive transport in field scale porous media flow systems[END_REF].

Performance of the method

The CPU times for the iterative splitting, non iterative splitting and global approaches are compared on figure 3. The CPU time required for each method is plotted versus the number of the nodes of the grid. As expected, it can be seen that the non-iterative method requires much less CPU time than the iterative methods. On the other hand, the global approach described in the paper requires less time than the iterative splitting, at least for the simple chemical system considered here. For a single time step, the iterative splitting approach requires between 20 and 27 iterations on the average. The number of fixed point iterations increase with the number of the nodes in the grid. On the other hand, the number of Newton iterations for the Newton-Krylov method is less than 6, independently of the number of nodes. The number of Krylov iteration for each Newton step, however, does increase with the number of nodes. We go back to this issue in subsection 4.2.1.

The 1D "easy" MoMaS Benchmark

The global and the splitting approaches will now be applied to the 1D easy GDR Momas Benchmark, as described in the introductory paper to this special issue [START_REF] Carrayrou | Presentation of the MoMaS reactive transport benchmark[END_REF], see also the original description in [START_REF] Carrayrou | GDR MoMaS benchmark -reactive transport[END_REF]. Let us just recall that the model is a one-dimensional column, made of 2 different media: the part in the middle is less conductive but more reactive than the surrounding medium. The chemical system has 5 components (4 mobile components and a fixed component), and 7 secondary species. The equilibrium constants vary over 50 orders of magnitude, and the stoichiometric coefficients can be as large as 4, making the problem highly non-linear.

First, results showing the evolution of the component species at various times, and using several spatial and temporal resolutions are shown on figure 4a. The left figure is at time t = 10, the right one at t = 50. As expected, the concentrations remain almost constant in the middle (reactive) region. Meshes with 220, 440, 660 and 880 points have been used, and in each case the time step is chosen as 0.9 times the limit fixed by the CFL condition. For these early times, the dependence on the mesh is not very strong. Elution curves (concentrations 5, first for t going from 0 to 400 (figure 5a), then for t going from 4900 to 5300 (figure 5b). The elution curves show that the correct limiting behavior is reached before the leaching phase begins.

The output results required in the benchmark definition are included. Most were obtained with a 220 points mesh, which may not be sufficient, as will be seen below. It has not yet been possible to obtain results with a finer mesh resolution for significantly longer times.

Figures 6a, figures 6b and 6c (elution curve for the total dissolved concentration of component X3,a nd species C1) show an oscillations pattern that has been observed by other groups working on the benchmark. These oscillations have been convincingly explained by V. Lagneau [START_REF] Lagneau | HYTEC results of the MoMas reactive transport benchmark[END_REF] as being due to the interaction of the very rapid chemistry and the discrete nature of the grid. They are a discretization artifact, but appear independently of the method. They can be reduced by using a more refined grid.

Figures 7a and7b show the influence on the mesh, by showing the concentration over a small spatial region, for time t = 10 . The concentrations are computed with 4 meshes of increasing resolution. The peaks in the solution are not resolved satisfactorily for the coarser mesh, with 220 points, but 660 (and better 880) points give the correct location and amplitude. Even if the method as it is currently implemented cannot yet be considered as robust, its ability to locate these solution features with reasonably coarse meshes was seen as one of its strong points. Unfortunately, this may still not be enough to eliminate the oscillations shown on figure 6. This issue is currently being worked on, part of the difficulty being that increasing the mesh resolution may not be sufficient. As the nonlinear problem becomes more difficult, it may be necessary to increase the maximum number of iterations allowed to make sure the Newton-Krylov method has converged.

Performance of the method

The benchmark was intended to be a difficult test for numerical methods, and this is indeed the case. On the average, more than 20 Newton iterations are required at each time step, and between 15 and 40 conjugate gradient steps are needed at each nonlinear iterations. Figure 8 shows a typical time step: the solid curve shows the cumulative number of conjugate gradient (alternatively, the number of matrix vector products), and the dots represent the nonlinear iterations.

Statistics for a single time step are gathered in table 3, for three different mesh resolutions (220, 440 and 660 points). They give the number of non-linear iterations (NNI) for a (typical) time step, and the total number of linear iterations (NLI) accummulated over the whole Newton iteration. The number of nonlinear iterations depends only weakly on the mesh resolution, whereas the number if linear iterations increases with the mesh resolution. Table 3 shows that the solver spends a large proportion of its time in the linear solver, despite the adaptive choice of the forcing parameter (equation ( 29)). Moreover, the number of linear iterations for each nonlinear iteration also increases with the mesh resolution. Actually, this is expected, as the solution of the linearized problem includes the solution of the transport operator, wich has an elliptic-like structure, so that its condition number grows like the square of the number of grid points. This problem could be alleviated by using a suitable preconditioner that would make the number of iterations independent of the mesh resolution (a domain decomposition preconditioner could be used as in [START_REF] Achdou | A domain decomposition preconditioner for an advection-diffusion problem[END_REF]). As noticed by Hammond et al. [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physicsbased preconditioning to biogeochemical transport[END_REF], designed a matrix-free preconditioner (so as to be compatible with the Newton-Krylov framework) is a challenge. Natural choices would exploit the block structure of the Jacobian, the simpler ones being based on block-Jacobi, or block Gauss-Seidel. Operator-splitting as a preconditioner has also been proposed in [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physicsbased preconditioning to biogeochemical transport[END_REF]. These possibilities are currently being explored, exploiting the block structure of the Jacobian, and the results will be reported in a forthcoming paper [START_REF] Taakili | Linear and nonlinear preconditioning for a model of transport in porous media with sorption, in preparation 45. van der Vorst HA[END_REF].

Conclusions-Perspectives

In this paper, it was shown that a global method for coupling transport with chemistry based on the Newton-Krylov technology can be implemented while keeping the transport and chemical solvers separated. The results shown are promising: it is possible to solve efficiently geochemical problems using the method, although there remains several issues that need to be addressed.

-The first is to run test cases on more demanding configurations, where the method can be expected to show its full potential. This includes the other MoMaS test cases, with a more complex chemistry model, and also an implementation of the method in 2 and 3 dimensions. -It will then certainly be necessary to explore the question of how to precondition the Jacobian, in order to reduce the number of Krylov iterations. An natural avenue is to reuse the operator splitting methods, as proposed by [START_REF] Hammond | Application of Jacobian-free Newton-Krylov with physicsbased preconditioning to biogeochemical transport[END_REF]. A similar study is being carried out for a related, but simpler model, see [START_REF] Taakili | Linear and nonlinear preconditioning for a model of transport in porous media with sorption, in preparation 45. van der Vorst HA[END_REF]. -The results reported above used a fixed time step, which was clearly insufficient for the large interval of integration. To successfully solve difficult problems like the benchmark above, it will clearly be necessary to use adaptive time stepping. -A more difficult problem will be to take into account precipitation-dissolution phenomena in the chemical model. As the models are non-differentiable, this makes it more difficult to employ Newton's method.

As was apparent from the numerical experiments, the method also shows some limitations. The most serious is its high cost, as each evaluation of the residual involves the solution of a chemical problem at each grid point. The fact that the method has two levels of nonlinear iterations means that it may not be as robust as other global methods based on a single level of iterations. Finding a good preconditioner may not be a limitation, but most strategies will involve solving more transport problems, which will also incur a high cost. by ANDRA, BRGM, CEA, EDF and IRSN. Both authors thank the referees for their detailed comments, which led to significant improvements in the contents of the paper.

Equation ( 3 )

 3 will be solved over the time step [t n ,t n+1 ] by first solving the advection equation φ ∂c ∂t + ∂ ∂x (uc) = 0 over M steps of size ∆t c each, and then solving the diffusion equation φ starting from the value at the end of the advection step.

Fig. 2 :

 2 Fig. 2: Elution curves (concentrations at the end of the column) versus time, for the problem of section 4.1. Left: global method, right: PhreeqC reference.

Fig. 3 :

 3 Fig. 3: Computing time for 3 methods applied to the ion exchange of section 4.1

(a) t = 10 (b) t = 50 Fig. 4 :

 10504 Fig. 4: Concentration of all components at times t = 10 and t = 50, for various mesh resolutions

Fig. 5 :Fig. 6 :

 56 Fig. 5: Concentration of the components X1 and X4 at the end of the column (x = 2.1) as a function of time

10 FigTable 3 :

 103 Fig. 7: Concentration profiles

Table 1 :

 1 Physical parameters

	Darcy velocity Diffusion coefficient Length of column Mesh size Duration of experiment Time step	2.78 10 -6 m/s 5.56 10 -9 m 2 /s 0.08 m 0.0002 m 1 day 720 s	Component Ca Cl K Na	C init 0 0 2.0 10 -4 1.0 10 -3	C inflow 0.6 10 -3 1.2 10 -3 0 0

Table 2 :

 2 Initial and injected concentrationsThe chemical reactions for this example are:

Acknowledgments

The first author would like to express her sincere thanks to ITASCA Consultants, France for providing the necessary Ph.D fellowship to carry out this research in INRIA-Rocquencourt, France. The second author's work was supported by Groupement MoMaS CNRS-2439. We gratefully acknowledge sponsorship of GDR MoMAS Fig. 8: Iterations


Marie-Christing Christine Lafaye
email: mclafaye@univ-lr.fr

Georges Louis
email: glouis@univ-lr.fr

Antoine Wiedemann

De UML au DDL, en passant par Java et un outil d'ORM

Keywords: ORM, base de données, intégrité des données, modèle du domaine, IDM, UML ORM, database, data integrity, domain model, MDE, UML

teaching and research institutions in France or abroad, or from public or private research centers.

Qualité des données gérées avec un ORM.

Introduction

Le modèle du domaine [START_REF] Larman | UML 2 et les Design Patterns, 2 ème edn[END_REF], noté ci-après MD est un modèle central au cours du développement des systèmes d'information, car il décrit les concepts du métier et pas une application particulière [START_REF] Ambler | The Object Primer : Agile Model-Driven Development with UML 2[END_REF]. Dans un modèle du domaine, le type d'un attribut ne doit pas être un concept complexe du domaine, implanté par une référence dans un langage à objets. Il est possible d'y utiliser non seulement les types primitifs d'UML (Integer, Boolean, String, UnlimitedNatural) mais aussi d'autres types courants comme Date, CodePostal, Adresse, types énumérés . . .. Il est comparable à un diagramme E/R [START_REF] Chen | The Entity-Relationship Model-Toward a Unified View of Data[END_REF]) a contrario des modèles objets dont les structures de données sont plus complexes [START_REF] Cauvet | Processus de conception orientée objet : transformation d'un schéma conceptuel en un schéma logique[END_REF]. Il s'agit d'un modèle multi-usages souvent utilisé comme support dans la conception des bases de données relationnelles. Nous voulons aider le concepteur à dériver d'un MD une implémentation Java des classes persistantes, un mapping ORM (Object Relational Mapping), et surtout le schéma relationnel qui en résulte. Nous faisons l'hypothèse que la cohérence des données stockées est garantie car d'une part nous déduisons des multiplicités des associations et des hiérarchies de spécialisation/généralisation du MD les contraintes référentielles, et d'autre part nous prenons en compte toutes les contraintes d'unicité collectées par le concepteur.

De nombreuses approches produisent le schéma relationnel à partir d'un diagramme E/R. Cependant, elles ne permettent pas la génération de code SQL acceptant les valeurs nulles [START_REF] Codd | Missing information (applicable and inapplicable) in relational databases[END_REF]. De plus UML est aujourd'hui largement répandu et nous souhaitons permettre au concepteur de travailler dans ce langage. On trouve dans la littérature de nombreux patrons de conception du schéma relationnel à partir d'un MD UML, mais la transformation n'est pas automatisée et la prise en compte exhaustive des contraintes d'unicité non assurée [START_REF] Dorsey | Oracle8 Database Design Using UML Object Modeling[END_REF][START_REF] Soutou | De UML à SQL[END_REF][START_REF] Urban | An advanced course in database systems, beyond relational databases[END_REF]. Les outils d'ORM, exploitent ces patrons (Hibernate, 2009), (Sun 2009), [START_REF] Fowler | Patterns of Enterprise Application Architecture[END_REF]. Ils cachent la base de données sous-jacente. La qualité du schéma relationnel dépend de la spécification du mapping. Nous proposons d'utiliser les techniques d'ingénierie dirigée par les modèles (IDM) pour obtenir automatiquement ce mapping.

Dans la section 2 nous présentons les principes de fonctionnement des outils ORM. Dans la section 3 nous décrivons notre proposition d'utilisation de l'IDM pour la transformation d'un MD en un schéma relationnel. Nous avons expérimenté l'outil ORM Hibernate, et nous détaillons dans la section 4 le mode d'utilisation des annotations disponibles pour que cet outil génère automatiquement le DDL cible souhaité. La section 5 est la conclusion.

Object Relational Mapping (ORM)

Le but d'un outil d'ORM est de fournir au développeur d'applications en langage à objet un service de persistance dans une base de données relationnelle. Faire fonctionner un application Java avec un support persistant relationnel suppose de résoudre les problèmes posés par les différences entre le modèle objet de la programmation et le modèle relationnel. Il faut par exemple rendre cohérentes les façons dont l'un et l'autre des paradigmes mettent en oeuvre l'identification des individus. Les objets voient leur identification gérée de manière transparente par le système, tandis que les tuples relationnels sont identifiés explicitement par la valeur des attributs composant leur clé primaire. L'association dans UML et la référence de la programmation objet s'appuient sur ces identifiants implicites, au contraire de la dépendance de référence du modèle relationnel qui repose sur des identifiants explicites.

Pour exécuter les ordres SQL insert, select, update et delete (CRUD) requis par les opérations CRUD sur les objets, les solutions ORM utilisent l'API Java Database Connectivity (JDBC). [START_REF] Fowler | Patterns of Enterprise Application Architecture[END_REF] propose d'inclure le code SQL soit directement dans la classe du domaine (patron Active Record), soit dans une classe annexe à la classe du domaine (patron Data Mapper ou patron Data Access Object (DAO) de Sun (Sun 2002)). Ce mapping fixe la correspondance entre les attributs des classes Java et les colonnes du schéma relationnel. Il peut être utilisé pour générer le DDL de la base de données, puis il pilote l'instanciation des objets et l'enregistrement des tuples par l'outil ORM (matérialisation/dématérialisation des objets de l'application [START_REF] Larman | UML 2 et les Design Patterns, 2 ème edn[END_REF]).

La plateforme Java EE 5 défini l'API Java Persistence (JPA) comme standard pour l'ORM (EJB 2006). JPA exploite des solutions présentes dans les plateformes de persistance répandues sur le marché comme Hibernate, Orable TopLink, Java Data Objects (JDO), et la persistance gérée par conteneur des Enterprise JavaBeans (EJB). L'API définit un modèle de persistance standard pour les applications Java. JPA admet que les métadonnées de mapping puissent être fournies de deux manières : soit par un fichier XML, soit par des annotations directement sur le code Java. Hibernate, développé depuis 2001, implante complètement JPA. Mais il possède un jeu d'annotations supplémentaire sur le code Java : l'API Hibernate Annotations (Hibernate, 2009). Ces annotations restent moins puissantes qu'un document XML pour spécifier un mapping.

Ingénierie dirigée par les modèles pour la conception du schéma relationnel à partir du modèle du domaine

Nous détaillons dans la figure 1 les différentes étapes du processus de conception assistée que nous proposons. Trois espaces techniques [START_REF] Bézivin | Model Driven Engineering : An Emerging Technical Space[END_REF] se côtoient. Dans l'espace UML le concepteur construit manuellement un MD. L'utilisation de stéréotypes UML permet de sélectionner les classes dont les objets seront persistants. Nous fournissons des annotations spécifiques pour exprimer les contraintes d'unicité sur le MD (un stéréotype « UC » sur une contrainte UML associée à une classe persistante). Nous tentons de garder ces annotations aussi légères et ergonomiques que possible. Ce MD annoté, dont un exemple est fourni figure 3, est la variable d'entrée des différentes transformations de modèle opérées dans l'espace IDM. Les règles de transformation que nous utilisons sont les suivantes : 2) Chaque attribut de type primitif devient une colonne sur laquelle une contrainte NOT NULL est définie par défaut. Nous acceptons les types primitifs UML et SQL (notamment Date).

3) Chaque association binaire est transformée soit en une clé étrangère, soit en une table de jointure avec deux clés étrangères (cf table 1). En fonction des multiplicités de l'association, il est nécessaire de créer des contraintes NOT NULL, UNIQUE et des triggers.

4) Chaque association n-aire devient une table de jointure composée de n clés étrangères.

5) La transformation des classes association s'effectue par l'application des règles 1 et 3 ou 1 et 4 ci-dessus. 6) Chaque généralisation est transformée en une clé étrangère (qui est aussi clé primaire) dans chaque spécialisation (ceci s'applique à l'héritage multiple). 7) Chaque contrainte « UC » (UML) est transformée en une contrainte UNIQUE (SQL). 8) Par défaut, l'identification par valeur des instances est réalisée en adoptant le patron de conception « Identity Field » [START_REF] Fowler | Patterns of Enterprise Application Architecture[END_REF] Nous obtenons deux types de résultats : une représentation CWM d'un schéma relationnel [START_REF] Lafaye | L'approche MDA pour la conception des bases de données relationnelles, une illusoire simplicité[END_REF][START_REF] Wiedemann | Relational Database Modeling[END_REF], ou un mapping ORM spécifié au choix sous forme d'un fichier XML ou d'annotations Java. Les règles de transformation que nous appliquons sont identiques dans les deux cas.

Traitement de l'identification et des contraintes d'unicité

Pour transformer les contraintes UML « UC » en contraintes SQL UNIQUE conformément à la règle 7 de la section 3, le mapping Hibernate contient des annotations @UniqueConstraint (lignes 19 et 34 de la figure 4).

Conformément à la règle 8 de la section 3, l'attribut artificiel sur lequel est définie la clé primaire est aussi ajouté à la classe Java, avec les annotations JPA @Id et @GeneratedValue (exemple id_TRAI lignes 38 et 39 figure 4). Pour définir une clé primaire composée, il est nécessaire de définir une classe interne Id et d'indiquer son rôle d'identifiant grâce aux annotations Embeddable (exemple Eval lignes 55 à 80).

Traitement des associations

Le traitement des associations comprend la définition de références Java adéquates dans le code des classes métier, et le mapping Hibernate. Nous distinguons trois familles d'associations binaires : 1-1, 1-N et N-M, comme c'est indiqué dans la première colonne de la table 1 Remarque 3 : Associations [1..1-1..1] Pour nous la meilleure implantation relationnelle de ces associations est une table résultante de la fusion des deux classes, et de la construction sur cette table de deux vues, recréant par projection les deux entités initiales. La BCNF [START_REF] Codd | Recent investigations in relational database systems[END_REF]

Généralisations/Spécialisations

Plusieurs stratégies sont proposées dans la littérature [START_REF] Ambler | Mapping Objects to Relationnal Databases : O/R Mapping in detail[END_REF][START_REF] Soutou | De UML à SQL[END_REF][START_REF] Urban | An advanced course in database systems, beyond relational databases[END_REF]

Conclusion

Lorsque plusieurs applications vont partager un même modèle du domaine, l'effort nécessaire pour définir ce modèle et pour lui donner une implémentation de qualité est tout-à-fait justifié. Nous avons montré ici une technique pour obtenir ce résultat en utilisant les patrons de conception et les outils habituellement appliqués pour le passage d'un modèle UML à des classes Java dont la persistance est assurée par une base de données relationnelle (ORM). Un des aspects sur lesquels nous insistons est la qualité du schéma relationnel, et l'obtention des contraintes qui assureront la cohérence des données enregistrées.

Notre proposition peut se résumer ainsi : 1) Annoter le modèle UML des données persistantes, en définissant des contraintes d'unicité, puis 2) Transformer systématiquement le modèle UML en classes Java, annotées de manière à 3) Obtenir un schéma relationnel de qualité, avec un maximum de contraintes d'intégrité.

La qualité du modèle relationnel obtenu se caractérise comme suit : 1) Les valeurs nulles, dans leur variété inapplicable au sens de [START_REF] Codd | Missing information (applicable and inapplicable) in relational databases[END_REF], sont utilisées pour limiter le nombre de tables de jointure (les seules associations binaires qui donnent lieu à une table sont les associations N-N). Ceci limite le nombre de jointures à l'exécution.

2) Les contraintes d'unicité collectées sur le modèle UML sont traduites en contraintes d'unicité (UNIQUE) dans le schéma relationnel.

3) On évite la multiplication de clés artificielles inutiles. C'est notamment le cas pour le traitement des généralisations/spécialisations d'UML.

Le traitement des valeurs nulles « inapplicables » est lié à la dynamique des objets. Spécifier cette dynamique au niveau du modèle UML devrait permettre de compléter notre outillage afin d'introduire la vérification automatique de nouvelles contraintes sur la base de données, qui assureraient que les valeurs nulles dans les tables restent bien en permanence cohérentes avec l'état des instances du modèle.

Figure 1 -

 1 Figure 1 -Diagramme d'activité du processus de conception, depuis le modèle du domaine jusqu'aux artefacts relationnels

  Dans cet article, nous détaillons les mappings qui serviront d'entrée dans l'espace ORM à la transformation permettant d'obtenir in fine le DDL cible. Ce DDL, équivalent à celui représenté dans le modèle CWM, contient la description des tables et des contraintes d'intégrité : NOT NULL, PRIMARY KEY, UNIQUE, FOREIGN KEY. . .et minimise le nombre d'identifiants systèmes (surrogate key).4. Spécification du mapping Hibernate Nous illustrons notre travail à l'aide du modèle du domaine présenté figure 3, dans lequel un stage est suivi par un étudiant, et donne lieu à une soutenance. Plusieurs critères sont utilisés pour évaluer un stage, chacun d'entre eux fait l'objet d'une note attribuée par un enseignant. Le concepteur a rajouté sur ce modèle des contraintes d'unicité qui seront exploitées dans le modèle relationnel. Tous les attributs de la base de données sont NOT NULL, sauf les attributs Training.end et Training.mark à cause de la multiplicité [0..1] dans le MD UML. Dans la figure 4, le modèle du domaine UML est implémenté en Java. Le mapping Hibernate y est ajouté par des annotations sur le code Java, et par un complément en XML. Le DDL obtenu est présenté figure 2. 4.1. Traitement d'une classe Le traitement d'une classe du modèle du domaine comprend toujours celui de ses attributs et de son identification. L'annotation JPA @Entity indique qu'une classe doit être persistante ( classe Training des lignes 34 à 54 de la figure 4). Dans le MD, les attributs possèdent par défaut une multiplicité [1..1]. Notre mapping ajoute donc par défaut l'annotation JPA @Column{nullable=false}.

Figure 3 -

 3 Figure 3 -Exemple de modèle du domaine annoté

Remarque 1 :

 1 navigabilité des associations UML L'association UML de type 1-N (Training-Student par exemple) est l'association la plus communément rencontrée en modélisation. En UML, elle est implicitement bidirectionnelle, de même, dans le modèle relationnel, puisque la jointure est commutative. La référence Java n'autorise cependant la navigation que dans un sens, de la classe qui possède la référence vers la classe référencée. Représenter le double sens de navigation en Java oblige à définir deux références mutuelles dans Training et Student. Pour tenir compte de la cardinalité maximum * du côté de Training dans le modèle métier UML, la référence Student.trainings est déclarée comme une collection. Mais contrairement à UML, où les deux extrémités sont liées par l'association à laquelle elles appartiennent, il n'en est rien en Java : les deux références définissent à elle deux non pas une association bidirectionnelle Training ↔ Student, mais deux associations unidirectionnelles, de sens opposé, Training → Student et Training ← Student. L'implémentation Java de la figure 4 ne contraint pas une instance t de Training appartenant à la collection s.trainings d'une instance de Student d'avoir sa référence t.trainee égale à s. Au mapping de la référence Training.trainee par l'annotation ManyToOne, lignes 47 à 51, et au mapping de Student.training par l'annotation OneToMany, ligne 25, il est nécessaire d'ajouter la propriété mappedBy="trainee" pour préciser que les deux références Java représentent les deux extrémités d'une même association.Remarque 2 : cardinalités minimum Nous proposons de traiter les cardinalités minimum à 1 des multiplicités 1-1 en ajoutant une contrainte supplémentaire, qu'Hibernate prend en charge avec la propriété nullable = false. La classe Training avec l'annotation JPA @JoinColumn, et sa propriété nullable = false, ligne 49 précise qu'une instance de Training doit être reliée à une instance de Student. Pour les cardinalités minimums à 1 des multiplicités 1-N il faudrait ajouter un trigger dans le code SQL (Hibernate n'en génère pas).

Remarque 4 :

 4 Association [0..1-1..1] Pour pouvoir faire en sorte que la clé étrangère soit aussi clé primaire (exemple : Defence), nous devons utiliser un mapping XML comme celui des lignes 82 à 89 de la figure 4.

  DDL cible et annotations Hibernate pour les associations binaires

  . Un attribut artificiel « id_<alias classe> » utilisé comme clé primaire est automatiquement géré par le SGBD (surrogate key). Dans certains cas, la clé primaire d'une table est construite autrement. Dans l'exemple de la figure 3, les spécialisations Teacher et Student héritent de l'identifiant de la généralisation Person. La classe association ternaire Eval possède un identifiant construit par la concaténation des clés primaires des trois classes associées. Enfin, compte tenu de la multiplicité [0..1-1..1] de l'association Training-Defence, la classe Defence utilise la clé étrangère référençant Training comme clé primaire. Quand une contrainte d'unicité définie par le concepteur est bâtie sur plusieurs attributs, nous pensons qu'elle ne constitue pas un bon candidat à l'identification. Sur la classe Training, la contrainte d'unicité {trainee, start} est composée. Ici, trainee est un rôle et a pu être utilisé par le concepteur pour la construction de sa contrainte en raison des multiplicités de l'association Training-Student.

  Pour garder la possibilité d'avoir dans l'application des instances Java dont les attributs NULL n'ont pas de valeur, nous utilisons des types objet et non des types primitifs. Par exemple, Training.mark est déclaré du type objet Java Integer, et non pas du type primitif int (ligne 46) et Eval.mark (ligne 80).

	Tables :		Dépendances de référence :
	create table Criterion ( ... ); create table Eval ( criterion int8 not null , teacher int8 not null , training int8 not null , mark int8 not null , primary key ( criterion , teacher , training ) ); create table Person ( ... ); create table Student ( no int4 not null , id_PERS int8 not null , primary key ( id_PERS ), unique ( no ) ); create table Teacher ( ... ); create table Training ( id_TRAI int8 not null , " end " timestamp , mark int4 , start timestamp not null , topic varchar (255) not null , trainee int8 not null , primary key ( id_TRAI ), unique ( trainee , start ) ); create table Defence ( training int8 not null , date timestamp not null , room varchar (255) not null , mark float4 , primary key ( training ) );	1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30	alter table Eval add constraint FK_TRAI foreign key ( training ) references Training ; alter table Eval add constraint FK_CRIT foreign key ( criterion ) references Criterion ; alter table Eval add constraint FK_TEAC foreign key ( teacher ) references Teacher ; alter table Student add constraint FKF3371A1B9D5CEDAD foreign key ( id_PERS ) references Person ; alter table Teacher add constraint FKD6A63C29D5CEDAD foreign key ( id_PERS ) references Person ; alter table Training add constraint FK_STUD foreign key ( trainee ) references Student ; alter table Defence add constraint FK_TRAI foreign key ( training ) references Training ;	31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
	Figure 2 -DDL PostgreSQL cible du schéma relationnel
	l'ajoutons pas dans les cas de multiplicité d'attribut [0..1] (exemple Training.end
	ligne 45 figure 4).
				Nous ne

  . La troisième colonne indique dans quelle table la clé étrangère représentant l'association est ajoutée. Par exemple, pour la multiplicité [0..1-0..1] le concepteur peut choisir d'implanter la clé étrangère soit dans la table obtenue par transformation pour la classe A soit dans la table pour la classe B. Pour les associations de type N-M, il est nécessaire de rajouter un table supplémentaire (table de jointure), avec deux clés étrangères. La quatrième colonne précise le cas où la clé étrangère admet des valeurs nulles, et la cinquième le cas où il nécessaire de créer une contrainte UNIQUE sur la clé étrangère. Les septième, et huitième colonnes détaillent le code et les annotations Hibernate qui doivent être ajoutés aux classes Java pour obtenir le DDL souhaité.

  sera ainsi respectée. JPA et Hibernate permettent l'obtention d'un table unique pour deux classes Java, mais pas la définition des vues.

  pour implanter une hiérarchie de spécialisations/généralisations dans un schéma relationnel. Hibernate permet trois d'entre elles : une table par hiérarchie, une table par classe et une table par spécialisation. Nous choisissons la stratégie une table par classe par l'attribut strategy de l'annotation Inheritance, ligne 13 de la figure 4.

	@Entity public class Criterion {	1
	...	2
	}	3
	// mapping fourni en XML dans un	4
	// fichier séparé	5
	public class Defence {	6
	long id_DEFE ;	7
	Date date ;	8
	Float mark ;	9
	String room ;	10
	Training training ;	11
	}	12
	@Entity @Inheritance ( strategy =	13
	InheritanceType . JOINED )	14
	public class Person {	15
	...	16
	}	17
	18 @Entity @Table ( uniqueConstraints =