Accuracy of operator splitting for advection‐dispersion‐reaction problems

An operator‐splitting approach is often used for the numerical solution of advection‐dispersion‐reaction problems. Operationally, this approach advances the solution over a single time step in two stages, one involving the solution of the nonreactive advection‐dispersion equation and the other the solution of the reaction equations. The first stage is usually solved with a finite difference, finite element, or related technique, while the second stage is normally solved with an ordinary differential equation integrator. The only generally published guidelines on numerical accuracy suggest that the discretization errors associated with each stage must be small in order to achieve high accuracy of the overall solution. However, in this note we demonstrate that there is an inherent mass balance error present in the operator‐splitting algorithm for problems involving continuous mass influx boundary conditions. The mass balance error does not exist for instantaneous mass input problems. These conclusions are based upon analysis of a simple first‐order decay problem for which each stage of the calculation can be performed analytically (i.e., without discretization error). For this linear decay problem we find that the product of the first‐order decay coefficient (k) times Δt must be less than approximately 0.1 in order for the mass balance error to be less than 5%. We also present a variant of the normal operator‐splitting algorithm in which the order of solving the advection‐dispersion and reaction operators is reversed at each time step. This modification reduces the mass balance error by more than a factor of 10 for a wide range of kΔt values.