A comparison of Eulerian and Lagrangian transport and non-linear reaction algorithms

Abstract When laboratory-measured chemical reaction rates are used in simulations at the field-scale, the models typically overpredict the apparent reaction rates. The discrepancy is primarily due to poorer mixing of chemically distinct waters at the larger scale. As a result, realistic field-scale predictions require accurate simulation of the degree of mixing between fluids. The Lagrangian particle-tracking (PT) method is a now-standard way to simulate the transport of conservative or sorbing solutes. The method’s main advantage is the absence of numerical dispersion (and its artificial mixing) when simulating advection. New algorithms allow particles of different species to interact in nonlinear (e.g., bimolecular) reactions. Therefore, the PT methods hold a promise of more accurate field-scale simulation of reactive transport because they eliminate the masking effects of spurious mixing due to advection errors inherent in grid-based methods. A hypothetical field-scale reaction scenario is constructed and run in PT and Eulerian (finite-volume/finite-difference) simulators. Grid-based advection schemes considered here include 1st- to 3rd-order spatially accurate total-variation-diminishing flux-limiting schemes, both of which are widely used in current transport/reaction codes. A homogeneous velocity field in which the Courant number is everywhere unity, so that the chosen Eulerian methods incur no error when simulating advection, shows that both the Eulerian and PT methods can achieve convergence in the L1 (integrated concentration) norm, but neither shows stricter pointwise convergence. In this specific case with a constant dispersion coefficient and bimolecular reaction A + B → P , the correct total amount of product is 0.221MA0, where MA0 is the original mass of reactant A. When the Courant number drops, the grid-based simulations can show remarkable errors due to spurious over- and under-mixing. In a heterogeneous velocity field (keeping the same constant and isotropic dispersion), the PT simulations show an increased reaction total from 0.221MA0 to 0.372MA0 due to fluid deformation, while the 1st-order Eulerian simulations using ≈ 106 cells (with a classical grid Peclet number Δx/αL of 10) have total product of 0.53MA0, or approximately twice as much additional reaction due to advection error. The 3rd-order TVD algorithm fares better, with total product of 0.394MA0, or about 1.14 times the increased reaction total. A very strict requirement on grid Peclet numbers for Eulerian simulations will be required for realistic reactions because of their nonlinear nature. We analytically estimate the magnitude of the effect for the end-member cases of very fast and very slow reactions and show that in either case, the mass produced is proportional to 1 / P e , where Pe is the Peclet number. Therefore, extra mass is produced according to D , where the dispersion includes any numerical dispersion error. We test two PT methods, one that kills particles upon reaction and another that decrements a particle’s mass. For the bimolecular reaction studied here, the computational demands of the particle-killing methods are much smaller than, and the particle-number-preserving algorithm are on par with, the fastest Eulerian methods.