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 effective 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 multirate interaction by exchange with the mean (MRIEM). We derive and analyze the macroscopic descriptionof 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. Key Points We develop a new approach to reactive transport modeling that accounts for the dynamics of small-scale concentration fluctuations The temporal evolution of integratedmixing metrics agrees with the characteristic trends of fully resolved systems Experimental observations of mixing-limited reactive transport are successfully reproduced