Newton-Krylov methods for coupling transport with chemistry in porous media

Abstract

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.