Newton-Krylov methods for coupling transport with chemistry in porous media
Creators: name: Kern, Michel
name: Amir, Laila
name: Kern, Michel
Publisher: name: XVI International Conference on Computational Methods in Water Resources
Year: 2006

description: 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 bio-remediation, and in the evaluation of underground water quality.

When looking at a coupled system, the role of chemistry is simply to separate the
species into mobile and immobile species (immobile species may come from either
sorption or precipitation reactions, and are not subject to transport). In this work,
we assume that the medium is saturated, and that surface reactions do not change the
porosity. 

Multi-species chemistry involves the solution of ordinary differential equations (if
the reactions are kinetic) or nonlinear algebraic equations (if we assume local
equilibrium). When simulating a coupled system, these equations have to be solved at
each (grid) point, and at every time (step), leading to a huge coupled non-linear
system. As has been observed several times, it is essential to use efficient
numerical methods so as to be able to handle the size of systems occurring in the
applications.

If we assume local equilibrium, than the coupled system may be written as a DAE.
Since methods and software for solving DAEs have reached a high level of maturity (at
least for small index systems, which is the case here), it is natural to try and use
this technology. One important issue is the size of the system to be solved, as all
chemical species at all grid points are coupled. For any realistic configuration, it
will not be possible to form, let alone factor, the Jacobian matrix. A better
solution is to use Newton--Krylov methods, where the linear system at each Newton
iteration is solved by an iterative method. We can thus keep the fast convergence of
Newton's method, while only requiring Jacobian matrix--vector products, and
these can be approximated by finite differences.

We will discuss two methods, based on different choices for the main unknowns: 

In the ``one level'' method, we take as unknowns both the total   dissolved and   
fixed concentrations, but also the concentrations for the component species that
are the main unknowns of the chemical problem. This is closely related to the Direct
Substitution Approach  used by the geochemists. By carefully formulating the problem,
we can put it into the DAE framework. At each time step, a single non-linear system
has to be solved, and the Jacobian matrix couples both the transport matrix, and the
chemistry Jacobian matrix. The method is called one-level as we solve transport and
chemistry at the same time. Provided that we can access this Jacobian matrices (and
this may not always be possible in practice), this may be less expensive than the
two-level methods below.

The second, ``two-level'' method takes as unknowns only the   total dissolved and
fixed concentrations. The system to be solved at each time step requires the
solution of the chemical sub-problem at each Newton iteration, when evaluating the
residual. On the orher hand, the Jacobian matrix again involves the transport matrix,
as  well as the Jacobain of the ``chemical solution operator'' mapping. This is only
known implicitly, and the Jacobian has again to be evaluated by finite diffferences.
Even though this method may be more expensive than the one-level method, 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. 

A major issue for both methods will be the choice of a suitable preconditioner.