Multicomponent reactive transport with chemical equilibrium reactions, involving 
advective and dispersive solute transport coupled with the nonlinear reaction, is 
fundamental feature of subsurface environments. At the equilibrium, the governing 
system of equations is formed by the partial differential equations (PDEs) of the 
transport operator and the nonlinear algebraic equations (AEs) describing the 
chemical reactions. These two sets of equations are coupled and needed to be solved 
using the operator splitting (OS) approach or the global approach. 
With (OS), the transport and reaction equations are separated and solved 
sequentially for each time step. The OS includes the Sequential non-iterative 
approach (SNIA) and the Sequential iterative approach (SIA), which iterates between 
transport and chemistry until convergence for each time step.

With the global approach, the governing equations of transport and reaction are 
solved simultaneously. We distinguish two methods in the global approach. (i) The 
Differential Algebraic Equations Approach (DAE) solves a large system formed by both 
transport and chemistry equations (Miller, 1983). (ii) The direct substitution 
approach (DSA) solves a reduced system obtained after substituting the chemistry 
algebraic equations in the transport partial differential equations (Shen and 
Nikolaidis, 1997; Jenning et al., 1982).

Since the reference work of Yeh and Tripathi (1989), the (OS) approach is widely 
used to simulate reactive transport problems since it allows different numerical 
methods to be used for the reactive and transport components. However, the (OS) 
approach can introduce operator-splitting errors which are avoided with the global 
approach.

In this work, we show that DSA and DAE have the same numerical behavior. For a fine 
discretization and/or for a large number of chemical species, DSA is shown to be 
more efficient then DAE. 
Both DSA and SIA give accurate results. Contrarily to DSA, SIA solves sequentially 
two small systems (transport and chemistry). It was shown in Saaltink (2001) that 
DSA runs faster than SIA in chemically difficult cases and the SIA may become faster 
than the DSA for very large, chemically simple problems. In this work, we combine 
DSA with a very efficient linear solver UMFPACK. Our numerical experiments suggest 
that for all cases, DSA is shown to be more efficient than SIA. DSA requires less 
iteration to reach the convergence and allows large time steps contrarily to SIA.

Comparisons between DSA and SNIA show that SNIA spends less CPU time than DSA. 
However OS errors introduced with SNIA are proportional to the time step size and 
can therefore be significant. When SNIA is combined with an adaptative time stepping 
procedure based on a posteriori control error, DSA can be more efficient then SNIA.
