The presence of solute concentration fluctuations at spatial scales much below the scale of resolution is a major challenge for modeling reactive transport in porous media. Overlooking small-scale fluctuations, which is the usual procedure, often results in strong disagreements between field observations and model predictions, including, but not limited to, the overestimation of e˙ective reaction rates. Existing innovative approaches that account for local reactant segregation do not provide a general mathematical formulation for the generation, transport and decay of these fluctuations and their impact on chemical reactions. We propose a Lagrangian formulation based on the random motion of fluid particles carrying solute concentrations whose departure from the local mean is relaxed through multi-rate interaction by exchange with the mean (MRIEM). We derive and analyze the macroscopic description of the local concentration covariance that emerges from the model, showing its potential to simulate the dynamics of mixing-limited processes. The action of hydrodynamic dispersion on coarse-scale concentration gradients is responsible for the production of local concentration covariance, whereas covariance destruction stems from the local mixing process represented by the MRIEM formulation. The temporal evolution of integrated mixing metrics in two simple scenarios shows the trends that characterize fully-resolved physical systems, such as a late-time power-law decay of the relative importance of incomplete mixing with respect to the total mixing. Experimental observations of mixing-limited reactive transport are successfully reproduced by the model.

In this work we employ a reduced-order basis of conservative chemical components to model reactive transport using a Lagrangian (particle tracking) method. While this practice is well-understood in the Eulerian (grid-based) context, its adaptation to a Lagrangian context requires a novel reformulation of particle transport properties. Because the number of conservative-species particles need not change during simulation, spatial resolution stays constant in time, and there is no increase in computational expense due to increasing numbers of product particles. Additionally, this treatment simplifies the interaction between equilibrium and kinetic reactions and allows the use of species-dependent transport operators at the same time. We apply this method to model a suite of simple test problems that include equilibrium and kinetic reactions, and results exhibit excellent match with base-case Eulerian results. Finally, we apply the new method to model a 2D problem concerning the mobilisation of cadmium by a CO$_2$ leak, showing the potential applicability of the proposed methodology.