A global method for coupling transport with chemistry in heterogeneous porous media
Abstract
Modeling reactive transport in porous media, using a local chemical equilibrium assumption, leads to a system of advectiondiffusion 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 nonlinear 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 NewtonKrylov 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 NewtonKrylov 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.
Keywords:
Geochemistry transport in porous media NewtonKrylov methods advection–diffusion–reaction equationsMsc:
76V05 65M99∎
1 Introduction
The simulation of multispecies 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 bioremediation, 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 multispecies 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 Yeh and Tripathi (1989), 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 Saaltink et al (2001), Saaltink et al (2000), Carrayrou et al (2004), Kanney et al (2003), Lucille et al (2000), Samper et al (2009). 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 timestep 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 nonlinear system is solved in one step, usually by some form of Newton’s method. Most papers use the Direct Substitution Approach (see Hammond et al (2005), Fahs et al (2008)), 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 de Dieuleveult et al (2009)). 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 Kraeutle and Knabner (2005, 2007), 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. Hammond et al (2005) have used the JacobianFree 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 NewtonKrylov 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 “black–box 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 (nonreactive) 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.
2 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.
2.1 Transport model
The transport of a single species through a porous medium (a domain , with or 3), with porosity , in a known Darcy field , subject to dispersion and molecular diffusion, follows the linear advection–dispersion equation
(1) 
where
is the transport operator, and is a source term. The diffusion–dispersion tensor is given by
where is the molecular diffusion coefficient, and (resp. ) 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 can be written as
(2) 
where the porosity and the diffusion–dispersion coefficient can both depend on space. Because the flow is assumed compressible, the velocity is taken to be a constant.
The initial condition is and, in view of the applications, the boundary conditions are a Dirichlet condition (given concentration) at the left boundary () and zero diffusive flux at the right boundary (). More general boundary conditions could easily be accommodated.
2.1.1 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 cellcentered finite volume scheme will be used, see for instance Eymard et al (2000). The interval is divided into intervals of length , where . For , denote by the center and the right end of element . Finally, denote by , the approximate solution in cell .
Equation (2) is written in the form
(3) 
where the flux has been split as the sum of a diffusive flux and an advective flux .
Equation (3) is integrated over a cell giving
(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):
(5) 
with
For the advective flux, an upwind approximation is used, so that (assuming ),
These approximations are corrected to take into account the boundary conditions, both at and at . The semidiscrete system can be summarized by the finite dimensional system
(6) 
where now represents the vector of cell concentrations, is the matrix form of the transport operator, is a mass matrix accounting for variable porosity and mesh size, is a give source term and represents the effects of the boundary conditions.
2.1.2 Time discretization
Let denote by the time step (taken constant for simplicity) used to discretize the time interval and denote by the (approximate) value of . The first and most straightforward alternative is to discretize equation (6) by the backward Euler method, see for instance Ascher (2008). 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 . 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. Siegel et al (1997) (see also Hoteit et al (2004); Mazzia et al (2000)). In this work, splitting is used only within the (linear) transport step, but recent papers by Kačur et al. Frolkovič and Kačur (2006); Kačur et al (2005) 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 will be used as the diffusion time step, it is divided into M time steps of advection such that where M >1, the advection time step will be controlled by CFL condition. Equation (3) will be solved over the time step by first solving the advection equation over steps of size each, and then solving the diffusion equation starting from the value at the end of the advection step.
Advection step
The interval is divided into intervals , , where . Denote the approximate concentration c at time and . The advection equation is discretized in time using the explicit Euler method to obtain
(7) 
Diffusion step
The diffusion part is discretized by an implicit Euler scheme, starting from :
(8) 
As above, 2 equations accounting for the boundary conditions must be added.
2.2 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 Rubin (1983)). This would not be the case if mineral dissolution was taken into account, as these reactions typically need kinetic models.
Consider a set of chemical species linked by reactions
where is the stoichiometric matrix. Following Morel Morel and Hering (1993), 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 )

mobile components ,

fixed components ,

mobile secondary species ,

fixed secondary species .
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 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
(9) 
where and are the equilibrium constants, and , and are the entries of the stoichiometric matrices , and .
Mass conservation for each component is expressed in the form
(10) 
where is the total concentration of the mobile component , and is the total concentration of the fixed component ( and are vectors of size and respectively). In the case of ion exchange, the second mass conservation equation is simply , and is the Cationic Exchange Capacity of the porous matrix (see Appelo and Postma Appelo and Postma (2005)). As will be seen later, in the context of coupled transport and chemistry, is given by the transport model and is constant. In a closed chemical system, 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 van der Lee (1993). An additional advantage has been pointed out by Samper et al. Samper et al (2009): 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 be the vector with entries , where are the entries of vector . Equations (9) can then be rewritten as a linear system
(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 , for all feasible values of the data and . 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 Shapiro and Shapley (1965)).
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. Kelley (1995)) 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 and . Define the function by
(12) 
where the notation for a vector means the vector with elements , then equations (10) and (11) are equivalent to:
(13) 
This is the nonlinear system that to be solved for and , given and . The secondary concentrations can then be computed from equation (11).
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 and its fixed part by
(14) 
Note that, by definition, the relationship 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 (13), then computing by (14). More precisely
(15)  
where equation (13) is first solved for and , then is computed by (11).
It is important to keep in mind that computing 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 itself has been computed, so that the nonlinear system (13) has been solved. First, the Jacobian matrix of 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:
(16) 
where is the diagonal matrix with vector along the diagonal. Then, by an application of the implicit function theorem (see for instance Rudin (1976)), and by differentiating equation (11), there comes
(17) 
It should be stressed that the Jacobian of 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 “blackbox” approach. Second, for more realistic chemical models, including nonideal 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.
2.3 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
(18) 
These equations can be derived from the individual conservation equations by standard algebraic manipulations, see for instance Yeh and Tripathi Yeh and Tripathi (1989). It is the formulation given in the benchmark definition Carrayrou et al (2009), see also Saaltink et al (1998), de Dieuleveult et al (2009). The second equation is obvious, as was taken as a constant (at each point in space).
Taking into account the relation noted above, the first equation of the system is equivalent to
(19) 
where is the total concentration, the total mobile concentration, and the total fixed concentration for component .
From now on, will denote the discretized transport operator, as defined in equation (6). 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 , where represents the spatial index and represents the chemical index, we shall denote by

the column vector of concentrations of species at all grid points;

the row vector of concentrations of all chemical species in grid cell .
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):
(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 stepsize is essential for efficiency). Denoting time indexes by a superscript, the following system is obtained
(21) 
This is the system to be solved at each time step.
3 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 , 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. Saaltink et al (1998), see also Salignac Salignac (1998), one can derive three main formulations from the system given in (20):

formulation (TC) where is the principal variable , the transported variable
(22) This is the formulation used by Erhel et al. in de Dieuleveult and Erhel (2007); de Dieuleveult et al (2009), as it lends itself best to a DAE type algorithm. It is not convenient for our purpose, as the transport equation then involves both and , and is thus not easily used with an existing transport solver.

formulation (TT) where is the principal variable , transported variable
(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 is the principal variable , transported variable
(24) This is formulation 4 in Saaltink et al. Saaltink et al (1998), 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 and 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 Frolkovič and Kačur (2006), or Kanney et al (2003), 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).
3.1 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 twosteps approach, as reviewed below (cf. Yeh and Tripathi (1989)). However, this method suffers from several defects: it may severely restrict the step size to ensure convergence, and if used noniteratively it is only first order in time, which may introduce additional errors (cf. Carrayrou et al (2004)). 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.
3.1.1 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 Yeh and Tripathi (1989), and also Kanney et al (2003), Lucille et al (2000), Saaltink et al (1998), Carrayrou et al (2004) or Mayer and MacQuarrie (2009)<. 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 GaussSeidel 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 NonIterative 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.
3.1.2 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)–(11) in equation (1) (this can be done explicitly, as in Hammond et al. Hammond et al (2005), or implicitly, as in Kräutle et al. Kraeutle and Knabner (2005, 2007), or Saaltink et al. Saaltink et al (1998)). 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. de Dieuleveult et al (2009), de Dieuleveult and Erhel (2009) or de Dieuleveult (2008) . A high performance parallel implementation is described by Hammond et al. Hammond et al (2005), 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. Kraeutle and Knabner (2005, 2007), and Hoffmann et al (2009). 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.
3.2 A Newton–Krylov based fully coupled method
As was already mentioned in the previous section, Hammond et al. Hammond et al (2005) 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
(25) 
where is known.
Denoting by the function
(26) 
the nonlinear problem to be solved at each time step is , where denotes the vector .
Recall that at each step of the “pure” form of Newton’s method for solving , one should compute the Jacobian matrix , solve the linear system (usually by Gaussian elimination)
(27) 
and then set . 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 , 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 Kelley (1995), Knoll and Keyes (2004) and Hammond et al (2005), 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 Knoll and Keyes (2004) that a NewtonKrylov solver can often be made by wrapping a classical splitstep solver. This is what is being done here, as the formulation to which the NewtonKrylov 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. Hammond et al (2005), or like the method presented in this section will be necessary.
A Krylov subspace method (see for instance Kelley (1995)) is used to approximately solve the linear system in equation (27). The linear iterates are drawn from the Krylov subspace, . In the GMRES method (see Saad and Schultz (1986)), the iterates are defined to minimize the residual over . Other methods, such as BiCGSTAB van der Vorst (1992) or QMR Freund and Nachtigal (1991) 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. Dembo et al (1982) 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
(28) 
in this context, as the initial iterate is usually 0. The choice of the forcing term should strike a balance between two conflicting goals:

Keep the (local) convergence of Newton’s methods;

Avoid oversolving, 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 , while the second one obviously tends to make larger. It has been shown (see theorem 6.1.4 in Kelley (1995)) that provided is bounded away from 1, the inexact Newton’s method will converge, and that superlinear convergence obtains if goes to zero faster than . Based on this result, the strategy proposed by Kelley in Kelley (1995) (after the choice in Eisenstat and Walker (1996)) computes as
(29) 
where is a parameter (the value suggested in Kelley (1995) is ). Safeguards are added to this choice in order to prevent 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 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 tradeoff is not advantageous. Indeed, it is well known that the most expensive part of the evaluation of 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 defined in (26). Given a vector , is first split into its three components, and each subvector is regarded as a matrix, as in section 2.3. Then is computed by block:

For the transport block, the transport operator is applied to each species , with a source term given by , for ( denotes at the previous time step);

The second block is the trivial computation ;

The third block is the solution of the chemical problem at each grid point: , for .
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 (that is the directional derivative of in the direction of the vector ) can be computed as
(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 can be written using the Kronecker product (see for instance Horn and Johnson (1990)) as . Then the Jacobian matrix is
(31) 
where is the Jacobian of , and for each , is a small by block.
The structure of the Jacobian is illustrated on figure 1, for the case , . It is a block matrix, each bock being of size . 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 is tridiagonal), and the chemistry part at the bottom has 10 blocks.
It would in principle possible to compute and store the Jacobian matrix according to equation (31) as a sparse matrix, and to compute the matrix–vector product using a general purpose routine. The advantage of the method given in equation (30) is that the structure of the Jacobian is fully exploited, which leads to a much more economical computation.
4 Numerical results
4.1 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 PHREEQC2 Parkhurst and Appelo (1999) as Example 11. The onedimensional 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 sodiumpotassiumnitrate solution in equilibrium with the cation exchanger. The column is then flushed with three 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 .
Darcy velocity  

Diffusion coefficient  
Length of column  
Mesh size  
Duration of experiment  1 day 
Time step 
Component  

Ca  0  
Cl  0  
K  0  
Na  0 
The chemical reactions for this example are:
NaX  
KX  
with NaX, KX and are (sorbed) complexes, and indicates exchange site with charge 1
4.1.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 steadystate value, the potassium is displaced from the exchanger and the concentration in solution increases to balance the 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. Xu et al (1999).
4.1.2 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 noniterative 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.
4.2 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 Carrayrou et al (2009), see also the original description in Carrayrou et al (2006). Let us just recall that the model is a onedimensional 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 nonlinear.
First, results showing the evolution of the component species at various times, and using several spatial and temporal resolutions are shown on figure (a)a. The left figure is at time , the right one at . As expected, the concentrations remain almost constant in the middle (reactive) region. Meshes with , , and 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 at the end of the column as functions of time) are shown on figure 5, first for going from 0 to 400 (figure (a)a), then for going from 4900 to 5300 (figure (b)b).
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 (a)a, figures (b)b and (c)c (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 Lagneau and van der Lee (2009) 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 (a)a and (b)b show the influence on the mesh, by showing the concentration over a small spatial region, for time . 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.
4.2.1 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 nonlinear 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.
Mesh 220  Mesh 440  Mesh 660  
NNI  NLI  NNI  NLI  NNI  NLI 
25  494  18  551  25  636 
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 ellipticlike 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 Achdou et al (2000)). As noticed by Hammond et al. Hammond et al (2005), designed a matrixfree preconditioner (so as to be compatible with the NewtonKrylov framework) is a challenge. Natural choices would exploit the block structure of the Jacobian, the simpler ones being based on blockJacobi, or block GaussSeidel. Operatorsplitting as a preconditioner has also been proposed in Hammond et al (2005). These possibilities are currently being explored, exploiting the block structure of the Jacobian, and the results will be reported in a forthcoming paper Taakili and Kern (2009).
5 Conclusions– Perspectives
In this paper, it was shown that a global method for coupling transport with chemistry based on the NewtonKrylov 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 Hammond et al (2005). A similar study is being carried out for a related, but simpler model, see Taakili and Kern (2009).

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 nondifferentiable, 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.
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 INRIARocquencourt, France. The second author’s work was supported by Groupement MoMaS CNRS2439. We gratefully acknowledge sponsorship of GDR MoMAS 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.
References
 Achdou et al (2000) Achdou Y, Tallec PL, Nataf F, Vidrascu M (2000) A domain decomposition preconditioner for an advectiondiffusion problem. Computer Methods in Applied Mechanics and Engineering 184(24):145 – 170, DOI DOI:10.1016/S00457825(99)002273
 Appelo and Postma (2005) Appelo CAJ, Postma D (2005) Geochemistry, Groundwater and Pollution, 2nd edn. CRC Press
 Ascher (2008) Ascher UM (2008) Numerical Methods for Evolutionary Differential Equations. Society for Industrial & Applied Mathematics
 Carrayrou et al (2004) Carrayrou J, Mosé R, Behra P (2004) Operatorsplitting procedures for reactive transport and comparison of mass balance errors. Journal of Contaminant Hydrology 68(34):239–268

Carrayrou et al (2006)
Carrayrou J, Dimier A, Kern M, Knabner P, Leterrier N (2006) GDR MoMaS
benchmark – reactive transport. published electronically at
http://www.gdrmomas.org/ex_qualifications.html
 Carrayrou et al (2009) Carrayrou J, Kern M, Knabner P (2009) Presentation of the MoMaS reactive transport benchmark. Computational Geosciences This issue
 Dembo et al (1982) Dembo RS, Eisenstat SC, Steihaug T (1982) Inexact newton methods. SIAM Journal on Numerical Analysis 19(2):400–408, DOI 10.1137/0719025, URL http://link.aip.org/link/?SNA/19/400/1
 de Dieuleveult (2008) de Dieuleveult C (2008) Un modèle numérique global et performant pour le couplage géochimietransport. Thèse de doctorat, Université de Rennes 1
 de Dieuleveult and Erhel (2007) de Dieuleveult C, Erhel J (2007) A numerical model for coupling chemistry and transport. In: International Conference on SCIentific Computation And Differential Equations, SciCADE 2007
 de Dieuleveult and Erhel (2009) de Dieuleveult C, Erhel J (2009) A global approach for reactive transport: application to the benchmark easy test case of MoMaS. Comput Geosci This issue
 de Dieuleveult et al (2009) de Dieuleveult C, Erhel J, Kern M (2009) A global strategy for solving reactive transport equations. Journal of Computational Physics 228(17):6395–6410, DOI doi:10.1016/j.jcp.2009.05.044
 Eisenstat and Walker (1996) Eisenstat SC, Walker HF (1996) Choosing the forcing terms in an inexact Newton method. SIAM Journal on Scientific Computing 17(1):16–32, URL citeseer.ist.psu.edu/article/eisenstat94choosing.html
 Eymard et al (2000) Eymard R, Gallouët T, Herbin R (2000) Finite volume methods. In: Ciarlet PG, Lions JL (eds) Handbook of Numerical Analysis, vol VII, North–Holland, pp 713–1020
 Fahs et al (2008) Fahs M, Carrayrou J, Younes A, Ackerer P (2008) On the efficiency of the direct substitution approach for reactive transport problems in porous media. Water, Air, & Soil Pollution pp 299–208, doi 10.1007/s1127000896912
 Freund and Nachtigal (1991) Freund RW, Nachtigal NM (1991) QMR: a quasiminimal residual method for nonHermitian linear systems. Numerische Mathematik 60:315–339
 Frolkovič and Kačur (2006) Frolkovič P, Kačur J (2006) Semianalytica solution of a contaminant transport equation with nonlinear sorption in 1D. Compuational Geosciences 10(3):279–290, DOI 10.1007/s1059600690239
 Hammond et al (2005) Hammond GE, Valocchi A, Lichtner P (2005) Application of Jacobianfree Newton–Krylov with physicsbased preconditioning to biogeochemical transport. Advances in Water Resources 28:359–376
 Hoffmann et al (2009) Hoffmann J, Kräutle S, Knabner P (2009) A parallel globalimplicit 2D solver for reactive transport problems in porous media based on a reduction scheme and its application to the MoMaS benchmark problem. Comput Geosci This issue
 Horn and Johnson (1990) Horn RA, Johnson CR (1990) Matrix Analysis. Cambridge University Press
 Hoteit et al (2004) Hoteit H, Ackerer P, Mosé R (2004) 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). Computational Geosciences 8:99–124, DOI doi:10.1023/B:COMG.0000035074.37722.71, URL http://www.ingentaconnect.com/content/klu/comg/2004/00000008/%00000002/05379190
 Kanney et al (2003) Kanney JF, Miller CT, Kelley CT (2003) Convergence of iterative split operator approaches for approximating nonlinear reactive transport problems. Advances in Water Resources 26(247261)
 Kačur et al (2005) Kačur J, Malengier B, Remešikova M (2005) Solution of contaminant transport with equilibrium and nonequilibrium adsorption. Comput Methods Appl Mech Engrg 194:479–489, DOI http://dx.doi.org/10.1016/j.cma.2004.05.017
 Kelley (1995) Kelley CT (1995) Iterative methods for linear and nonlinear equations, Frontiers in Applied Mathematics, vol 16. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, with separately available software
 Knoll and Keyes (2004) Knoll DA, Keyes DE (2004) Jacobianfree NewtonKrylov methods: a survey of approaches and applications. J Comput Phys 193(2):357–397, DOI http://dx.doi.org/10.1016/j.jcp.2003.08.010
 Kraeutle and Knabner (2005) Kraeutle S, Knabner P (2005) A new numerical reduction scheme for fully coupled multicomponent transportreaction problems in porous media. Water Resources Research 41(W09414), DOI doi:10.1029/2004WR003624
 Kraeutle and Knabner (2007) Kraeutle S, Knabner P (2007) A reduction scheme for coupled multicomponent transportreaction problems in porous media: Generalization to problems with heterogeneous equilibrium reactions. Water Resources Research 43(W03429), DOI doi:10.1029/2005WR004465
 Lagneau and van der Lee (2009) Lagneau V, van der Lee J (2009) HYTEC results of the MoMas reactive transport benchmark. Comput Geosci This issue
 van der Lee (1993) van der Lee J (1993) CHESS, another speciation and surface complexation computer code. Tech. Rep. LHM/RD/93/39, CIG École des Mines de Paris, Fontainebleau
 Lucille et al (2000) Lucille PL, Burnol A, Ollar P (2000) Chemtrap: a hydrogeochemical model for reactive transport in porous media. Hydrological processes 14:2261–2277
 Mayer and MacQuarrie (2009) Mayer K, MacQuarrie K (2009) Formulation of the multicomponent reactive transport code MIN3P and implementation of MoMaS benchmark problems. Comput Geosci This issue
 Mazzia et al (2000) Mazzia A, Bergamaschi L, Putti M (2000) A timesplitting technique for the advectiondispersion equation in groundwater. J Comput Phys 157(1):181–198, DOI http://dx.doi.org/10.1006/jcph.1999.6370
 Morel and Hering (1993) Morel FMM, Hering JG (1993) Principles and Applications of Aquatic Chemistry. Wiley, NewYork
 Parkhurst and Appelo (1999) Parkhurst DL, Appelo C (1999) User’s guide to PHREEQC (version 2) A computer program for speciation, batchreaction, onedimensional transport, and inverse geochemical calculations. Tech. Rep. 994259, USGS
 Rubin (1983) Rubin J (1983) Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions. Water Resources Research 19:1231–1252
 Rudin (1976) Rudin W (1976) Principles of mathematical analysis, 3rd edn. McGrawHill, New York :
 Saad and Schultz (1986) Saad Y, Schultz MH (1986) GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci Statist Comput 7(3):856–869
 Saaltink et al (1998) Saaltink M, Ayora C, JCarrera (1998) A mathematical formulation for reactive transport that eliminates mineral concentrations. Water Resources Research 34(7):1649–1656
 Saaltink et al (2000) Saaltink M, Carrera J, Ayora C (2000) A comparison of two approaches for reactive transport modelling. Journal of Geochemical Exploration 6970:97–101
 Saaltink et al (2001) Saaltink M, Carrera J, Ayora C (2001) On the behavior of approaches to simulate reactive transport. Journal of Contaminant Hydrology 48:213–235
 Salignac (1998) Salignac AL (1998) Transport multiespèces et réactions géochimiques en aquifère : développement et validation du modèle couplé HYTEC 2D. PhD thesis, École des Mines de Paris
 Samper et al (2009) Samper J, Xu T, Yang C (2009) A sequential partly iterative approach for multicomponent reactive transport with CORE2D. Computational Geosciences 13:301–316, DOI 10.1007/s1059600891195, URL http://dx.doi.org/10.1007/s1059600891195
 Shapiro and Shapley (1965) Shapiro NZ, Shapley LS (1965) Mass action laws and the Gibbs free energy function. J SOC Indust Appl Math 13(2):353–375
 Siegel et al (1997) Siegel P, Mosé R, Ackerer P, Jaffré J (1997) Solution of the advectiondispersion equation using a combination of discontinuous and mixed finite elements. Journal for Numerical Methods in Fluids 24:595–613
 Taakili and Kern (2009) Taakili A, Kern M (2009) Linear and nonlinear preconditioning for a model of transport in porous media with sorption, in preparation
 van der Vorst (1992) van der Vorst HA (1992) BiCGSTAB: A fast and smoothly converging variant of BiCG for the solution of nonsymmetric linear systems. SIAM Journal Sci Stat Comput 13(2):631–644
 Xu et al (1999) Xu T, Samper J, Ayora C, Manzano M, Custodio E (1999) Modeling of nonisothermal multicomponent reactive transport in field scale porous media flow systems. Journal of Hydrology 214:144–164
 Yeh and Tripathi (1989) Yeh GT, Tripathi VS (1989) A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components. Water Resources Research 25:93–108