Lagrangian Modeling of Mixing‐Limited Reactive Transport in Porous Media: Multirate Interaction by Exchange With the Mean

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.


Introduction
The inherent difficulty of properly representing the interaction of reactive chemicals occurring over multiple spatio-temporal scales in complex hydrodynamic settings renders reactive transport modeling in porous media a major challenge in subsurface hydrology [Dentz et al., 2011;Sanchez-Vila and Fernàndez-Garcia, 2016;Benson et al., 2017;Valocchi et al., 2019].With the exception of highly idealized settings or incredibly small samples, generally in porous media it is unfeasible to obtain a completely resolved flow field within real porous media geometries based on the complete microscopic equations (e.g., Navier-Stokes and Advection-Diffusion). Instead one typically describes the system with macroscopic equations in an equivalent continuum [Icardi et al., 2019, and references therein].By doing so, one essentially ignores detailed resolution of local velocity and concentration fluctuations occurring at the pore-scale, below the scale of the equivalent continuum.The system is represented by macroscopic variables and properties, which aim to represent subscale fluctuations in an effective manner, obtained for instance by volume averaging [Quintard and Whitaker, 1994;Whitaker, 1999;Wood et al., 2003].However, these effective parameters really only aim to capture mean behaviors, and processes that depend nonlinearly on subscale fluctuations may often not be well described.Similarly, since macroscopic properties such as the hydraulic conductivity can vary substantially in space within real aquifers, one may further upscale flow and transport in heterogeneous media with a new set of effective parameters [Dagan, 1989;Gelhar, 1993;Rubin, 2003].This step further reduces the apparent complexity of the system, but again does not contain potentially important information below the scale of the effective model.Among available macroscopic models, the upscaled advection-dispersion-reaction equation (ADRE) is the most widely used for modeling reactive transport at all practical spatial scales.It is embedded as the standard in most popular reactive transport codes [e.g., Cederberg et al., 1985;Mangold and Tsang, 1991;Yeh and Tripathi, 1991;Steefel and Lasaga, 1994;Walter et al., 1994;Saaltink et al., 2004;De Simoni et al., 2005;Bea et al., 2009;Steefel et al., 2015, and references therein].However, field and laboratory observations, numerical simulations and theoretical developments have demonstrated time and time again that the upscaled ADRE fails to adequately represent mixing and chemical reactions at all scales [Rashidi et al., 1996;Cao and Kitanidis, 1998;Gramling et al., 2002;Palanichamy et al., 2009;Tartakovsky et al., 2008;Fernàndez-Garcia et al., 2008;Edery et al., 2009;Sanchez-Vila et al., 2010;de Anna et al., 2014a,b;Porta et al., 2016], because of its disregard for the local concentration fluctuations and the use of scale-averaged concentrations to compute reactions.In fact, there is the general consensus that reaction rates observed in the field tend to be much lower than those measured in laboratory experiments because of the presence of anti-correlated local fluctuations of reactant concentrations [e.g., Chiogna and Bellin, 2013;Ding et al., 2017].This issue is often encapsulated as the inability of the classical advection-dispersion model to distinguish mixing from spreading.
In order to obtain better predictions, effective transport models should somehow incorporate the sub-scale mixing limitation effects.Several such approaches have been proposed in recent years, both from the Eulerian and from the Lagrangian perspective (see Porta et al. [2016] and references therein).The Eulerian approaches typically focus on very specific initial and boundary conditions, corresponding to the mixing of two reactants moving across a column-shaped porous medium, forming a sharp interface at = 0, as in the famous laboratory experiments of Gramling et al. [2002].As such, existing effective solutions typically contain a time-decaying term controlling either an apparent kinetic reaction rate [Sanchez-Vila et al., 2010], a pre-defined concentration covariance function [Chiogna and Bellin, 2013], or a mobile-mobile mass exchange rate coefficient [Ginn, 2018].Hochstetler and Kitanidis [2013] consider a constant, Damköhler-dependent efficiency term multiplying the reaction rate, which accounts for the effect of reactant segregation.While all the above-mentioned approaches provide interesting simplified interpretations of the physical process, they do not provide general differential equations governing the transport of local concentration fluctuations, and hence they might not be readily applicable to broader sets of initial and boundary conditions.
An exception to this is the pioneering work by Kapoor and Gelhar [1994], who proposed a concentration variance conservation equation.However, the authors did not attempt to export their theories to multi-component reactive transport.While seminal, the work's main weakness is arguably its reliance on the assumption of a stationary concentration micro-scale, which should in reality be subjected to aging as acknowledged later by Kapoor and Kitanidis [1998].Oates [2007] applied some of these concepts to model reactive transport in porous media, by imposing a global time-dependent concentration micro-scale that stabilizes at late times.
On the other hand, Lagrangian approaches have been proposed to reproduce mixing-limited reactive transport, motivated by hydrologic settings [Benson and Meerschaert, 2008;Edery et al., 2009;Ding et al., 2013;Paster et al., 2013Paster et al., , 2014;;Schmidt et al., 2017Schmidt et al., , 2018Schmidt et al., , 2019;;Benson et al., 2019].These emulate reactant segregation by specifying the covariance structure of the initial condition, which usually involves employing a determined finite number of particles (initial positions) in the simulation.However, it is unclear whether these approaches are capable of reproducing the transient evolution of concentration fluctuations in realistic settings, where such fluctuations may be actively promoted by local-scale heterogeneous advection, even if the initial or boundary conditions are perfectly smooth [e.g., de Dreuzy et al., 2012].
Here we present a novel Lagrangian approach to simultaneously account for ( ) coarse-scale advective-dispersive behavior as well as ( ) the generation, transport and decay of local concentration fluctuations.The model aims to offer not just a solution specific to one setting, but rather a mathematical framework to potentially represent a broad array of settings and transport problems.Unlike the aforementioned Lagrangian approaches, the proposed model does not rely on a noisy initial condition to represent reactant segregation, but in fact converges to the desired solution with sufficient particles (i.e., the particle number is not a model parameter, but rather a numerical discretization).Eulerian (two-moment) implementations of the proposed model are possible, but Lagrangian implementation is more natural and straightforward.
The paper is structured as follows.In §2 we develop the conceptual and mathematical model leading to the differential equation describing the local concentration fluctuations perceived by a random-walking Lagrangian particle.In §3 we derive the resulting Eulerian differential equation describing the transport, generation and decay of concentration point-covariance; we also provide the temporal evolution of the spatial integral of the former (or mixing state) in two simple cases of initial and boundary conditions with pseudo-analytical solution.In §4 we implement the proposed model to reproduce the reaction product concentration data corresponding to two benchmark laboratory experiments.Finally, in §5 we summarize our main conclusions.

Conservative transport and mixing
By definition, all continuum-scale models of transport in porous media assume or solve a flow field with some degree of coarse-graining; that is, the velocity variability below a threshold resolution (corresponding to some representative elementary volume) is removed and replaced by an upscaled dispersion tensor.We distinguish two spatial scales, above and below the aforementioned threshold, which hereafter we refer to as coarse scale and local scale, respectively.Coarse-scale concentrations of species A at position x and time , A (x, ), are often assumed to obey the upscaled advection-dispersion equation, where v is the coarse-scale velocity, and D is the dispersion tensor, which represents the combined effect of velocity fluctuations at the local scale (around v) and molecular diffusion.(1) assumes that the porosity (volume of fluid per unit volume of medium) is constant.Hereafter, we also assume that v and D are spatially and temporally constant.These assumptions simplify the presentation and analysis of the model, but generalization should be readily possible.One manner for solving equation ( 1) is via Random Walk Particle Tracking (RWPT) [e.g., Salamon et al., 2006], a Lagrangian approach in which particles = 1, . . ., carry solute mass of one or several chemical species, and their trajectory over small time intervals [ , + Δ ] is defined as a combination of deterministic advective displacements and a Wiener random process emulating dispersion, where X is the position of particle , B is a matrix such that BB T = D, and is a vector of random numbers drawn independently from a standard normal distribution.Here, similar to Benson and Bolster [2016] and Engdahl et al. [2017], each particle is assigned a static mass of solvent, , and a variable concentration of solute A, A, ; therefore the mass of A carried by is A, .
Given any particle attribute , one may define its interpolation onto the Eulerian space (x) [Monaghan, 2005] at any point x in the model domain Ω , being the number of spatial dimensions, as where summation is implied over all particles = 1, . . ., ; (x) is a smoothing function (kernel) whose size is scaled by a bandwidth ℎ [see, e.g., Wand and Jones, 1994]; and ( is the particle density.Note that operation (3) computes a weighted average of over particles located at X ≈ x.We refer to (x) as local average in the ideal limiting case where → ∞ and ℎ → 0. In this case, the kernel (x) approaches the -dimensional Dirac delta function (x).
In numerical applications (i.e., for finite ), o ne m ust u se a n averaging b andwidth ℎ large enough to avoid issues related to subsampling [e.g., Sole-Mari and Fernàndez-Garcia, 2018;Sole-Mari et al., 2019a].Hence, in this case we refer to (x) resulting from any numerical approximation of operation (3) as small-volume average, which is assumed to stabilize and approach the ideal local average given a sufficiently large and correspondingly small ℎ.More details are given in Appendix B.
given particle motion equation (2) and in the limit where → ∞, Δ → 0, the averaged concentrations A (x, ) converge to being governed by the Fokker-Planck equation [Risken, 1989], which is equivalent to the ADE (1) when D is spatially constant; otherwise a correction can be applied to the drift term, see LaBolle et al. [1996].
To summarize, in the proposed Lagrangian model, each numerical particle represents a discrete amount of a chemical solution traveling through a porous medium, moving by displacements representing the scale-averaged advection (deterministic) and upscaled dispersion (normal random).Consequently, at the coarse scale, the concentration field o beys t he a dvection-dispersion equation (ADE).This type of Lagrangian model is widely used by researchers and practitioners in hydrology to simulate nonreactive transport of solutes.
While coarse-scale concentrations of non-reactive chemicals may agree reasonably well with the ADE under certain conditions [Dagan, 1984], concentration fluctuations m ay s till o ccur at the local scale.These local-scale fluctuations, w hich a re n ot explicitly a ccounted for i n classical formulations, may drive the outcome of nonlinear processes, such as chemical reactions, far from what would be predicted by the ADE [Kang et al., 2019].Note that by using equation (2) (or similar stochastic formulations), where each particle follows its own unique random path, it is implied that at any given time each particle is sampling the local-scale fluid velocity fi eld independently.Analogously, in the proposed model, particle concentrations A, ( ) are assumed to represent the local-scale concentrations, and may therefore be at disequilibrium with the averaged A (x, ).Hence, hereafter we refer to A, ( ) as local concentrations.Figure 1 is a schematic representation that illustrates our proposed conceptual model.The local-scale structured spatial variability of concentrations in the physical system is emulated by the stochastic variability of local concentrations experienced by overlapping Lagrangian particles in the model.Because it is defined by interpolation (see Appendix B), the coarse-scale concentration is a smooth function in the Eulerian space, whereas local departures from the well-mixed equilibrium or local fluctuations are only defined on the Lagrangian particles.In order to represent the evolution of these local fluctuations we need to define a mixing model.
A simple representation of the local mixing as seen by a particle could be to assume a diffusion process driven by a coefficient within a fluctuation s tructure o f d imension and characteristic mixing length ℓ , with the subscript designating local (micro) scale variables,  [Villermaux, 1972;Pope, 2000].Its practical numerical implementation requires some careful consideration, related to features such as mass conservation.Implementation aspects, including small-volume approximations of the averaging operator, are discussed in Appendix B. An appealing advantage of this kind of mixing model is that a local value's variation in time depends only on its current degree of departure from the mean, potentially avoiding direct particle-particle interaction.Importantly in our context, the process is Markovian, i.e., the time-derivative (6) depends only on the current state.We note, however, that equation ( 6) is overly simplistic.Previous attempts to apply the IEM model (from an Eulerian perspective) to laminar flow and transport in heterogeneous porous media have concluded that, at the beginning of new contact between solutions with different chemical composition, one should account for a growing stage of ℓ before it reaches a stable asymptotic value [Kapoor and Kitanidis, 1998;de Dreuzy et al., 2012].
That is, a single constant value of cannot reproduce the distinct stages of the mixing process, and one should consider not only a slow linear mixing, but also a fast stretching-enhanced mixing stage.Moreover, fluctuations may occur across multiple overlapping length s cales.As acknowledged by Villermaux [1983] in the context of IEM applied to turbulent mixing, "several stages for mixing, each with their own time constants should be considered, possibly in series or in parallel".Here, in order to account for the impact of age on the mixing process, we propose a parallel multi-rate interaction by exchange with the mean (MRIEM), based on representing local mixing as occurr ing within different virtual zones = 1, . . ., Z , each being sampled by a fraction of the particle ( = 1).Within each zone , each particle sees a virtual local concentration A, , ( ) of each species A, possibly at disequilibrium with A (X ( ), ), such that The values of and are assumed to depend on local-scale flow and transport conditions.An idealized local mixing process representation that would resemble such a multi-rate behavior is given in Appendix A. The zone-concentration values A, , ( ) do not necessarily have a physical meaning individually, but are instead intended to emulate the complex transient nature of the mixing process.In principle, parameter sets and can be different for each species to account, for instance, for different values of the local-scale diffusion coefficient.
Given any Lagrangian-defined attribute and its local average across the particle space , it can be shown (see Appendix C) that, if particles move according to (2), the following relation holds between the Eulerian and the Lagrangian time-derivatives of : Then, by combining ( 10) with ( 8) and ( 9) we see that That is, the local mixing process described by (9) does not modify the coarse-scale description of non-reactive transport, driven by the particle displacements in (2).
One of the simplest implementations would comprise only two zones, one of them with an instantaneous mixing rate (i.e., very fast in relation to the time scale of interest), Hereafter, we refer to this particular case as dual-rate, to as the slow mixing fraction, and to as the slow mixing rate coefficient.In this case the local concentration in zone 1 is always at equilibrium with the coarse-scale concentration.Combining equations ( 8), ( 9) and ( 12) we may write: where for conciseness from here on, inside the derivative we use the notation A, ( ) ≡ A (X ( ), ).
In the dual-rate model (12), equation ( 13) can be interpreted as such: While following its random path described by (2), fluid-particle may experience variations in the perceived coarse-scale concentration of A. Conceptually, these changes correspond to the particle seeing itself involved in new mixing events, that is, in the formation of new fluctuation structures.Only a portion Figure 2 illustrates the conceptual decoupling of transport as a combination of spreading and mixing, in the true physical system as well as in the proposed Lagrangian model in one coarse-scale dimension.In the physical system, spreading represents the growth of the width (variance) of a solute plume due to velocity variability, which, alone, does not generate new contact between otherwise segregated solute molecules.Mixing, on the other hand, is precisely the generation of new contact between formerly segregated solutes, and is the result of local dispersion applied to the structure generated by spreading (see upper-right part of Figure 2).This decoupled picture is a simplification, because there is a continuous interplay between the two processes.The rate of spreading is influenced non-linearly by local dispersion [e.g., van Milligen and Bons, 2012].Similarly, mixing is influenced by the growth of contact surfaces, which is controlled by local advection [e.g., Villermaux, 2012].In the proposed Lagrangian model, spreading is represented by particle motion (2), which controls the coarse-scale behavior of concentrations; mixing is represented by the relaxation equation ( 9), which mitigates local departures from equilibrium experienced by individual particles, arising due to the aforementioned random motion.

Reactive transport
The proposed Lagrangian model can be extended to reactive transport applications by following the premise that chemical reactions occur at the local scale and thus are controlled exclusively by local concentrations defined on Lagrangian particles.We provide a brief summary on the incorporation of kinetic transformation ( §2.2.1) and equilibrium speciation ( §2.2.2).Naturally, the two may be integrated together, such as in Molins et al. [2004].

Kinetic reactions
Consider multiple kinetic reactions labeled = 1, . . ., R , with reaction rate laws that model reactions as a function of solute concentrations, (C), where C ≡ [ A , B , . . .] T ; and stoichiometric coefficients A, , B, , . . ., which indicate the generation/consumption of concentration per unit extent of reaction.Equation ( 9) is then extended to: where By combining ( 14) with ( 10) we obtain the coarse-scale Eulerian description Equation ( 16) elucidates that, in the occurrence of local concentration fluctuations and for nonlinear reaction systems, i.e., R A (C) ≠ R A (c), the Eulerian description of reactive transport does not obey the classical form of the advection-dispersion-reaction equation (ADRE) where reactions are computed directly from averaged concentrations.This is in contrast with the conservative transport case, where as shown by equation ( 11), coarse-scale concentrations do follow the classical ADE.

Equilibrium reactions
In the case of equilibrium reactions, a common approach is to compute the transport of chemically conservative components [e.g., Saaltink et al., 1998;Molins et al., 2004;De Simoni et al., 2005], and then speciation is provided by solving equilibrium system, A, ( ) = E A (U ( )). ( where E A (U) combines the law of mass action and the different stoichiometries to find the equilibrium concentrations, from components U ≡ { A , B , . . .}.In ( 17), U follow the conservative transport and mixing model presented in §2.1.By taking the particle average on both sides of ( 17), we obtain the coarse-scale description of the equilibrium reaction system We note from (18) that, similar to the kinetic reaction example, A ≠ E A (u), as opposed to classical well-mixed reactive transport approaches.

Concentration fluctuation dynamics
The dynamics of local concentration fluctuations are the key feature of the proposed reactive transport model.These local fluctuations, which emulate the sub-scale concentration variability of the porous medium, are what drives reactive transport away from the classical ADRE, as discussed in §2.2.Therefore, the analysis of their behavior is crucial in order to understand how the proposed model may be able to reproduce the mixing dynamics of physical systems, including those modeled in §4.
We study the behavior of concentration fluctuations through the local concentration covariance of two chemically conservative compounds A and B (where the particular single-compound case is implicitly included as A ≡ B).First, in §3.1, we present and discuss the partial differential equation describing covariance generation, transport and destruction.Then, in §3.2, we study integrated mixing metrics for two specific cases with closed form solutions, facilitating the comprehension of the model's behavior.These analytical solutions may also be used in future work for quantitative comparison to fully-resolved simulations, a task which is out of this work's scope.

Local concentration covariance PDE
Let us define the local fluctuation as the departure from well-mixed equilibrium on particles, A, ( ) A, , ( ) By definition, A (x, ) = 0. We study the concentration covariance of species A and B, which we denote as AB (x, ), As noted in §2.1, the assumption that the mixing dynamics of A and B can be described by the same sets of parameters { 1 , . . ., Z }, { 1 , . . ., Z }, is made here only for the sake of simplicity, and this assumption could be relaxed.It can be shown (see Appendix D.1) that, in the proposed model, if A and B are chemically conservative, each " " entry of the local concentration covariance evolves according to: The total local concentration covariance can be then obtained as the sum of all entries, as indicated by ( 21).Let us consider, once again, the specific case represented by ( 12).Then, one may write: It is worth remarking that, for B = A, expression ( 23) is mathematically equivalent to the concentration variance conservation equation introduced by Kapoor and Gelhar [1994, equation 56], as long as one assumes a scalar proportionality in their proposed dual-dispersivity system such that fulfills A = 2 ( + A), where and A are the microdispersivity and macrodispersivity tensors, respectively.Hence, the two conceptual models have clear similarities since in our case is the fraction of non-instantaneous mixing, and in Kapoor and Gelhar [1994], by analogy, it is the square root of the fraction of total dispersivity that is attributed to the macrodispersivity, i.e., to the nonmixed spreading.Nevertheless, there are important nuances that distinguish the models, as will be discussed in §4.1.
In any case, determining concentration variance/covariance is not the main, or at least not the only, purpose of our model.As outlined in §2.2, the distribution of local concentrations represented by the particles affects local processes such as chemical reactions.The concentration covariance is a powerful tool to assess contact between solutions, and therefore a good proxy for the potential magnitude of incomplete mixing effects on chemical reactions.

Analytical solutions
A metric that is commonly used to characterize spatial fluctuations of solute concentrations is the so-called mixing state [Bolster et al., 2011;de Dreuzy et al., 2012], which is defined as the spatial integral of the squared concentrations.Here we extend this definition, for any two solutes, A and B, as the spatial integral of the product of the two concentrations.In the present two-scale context, this may be written as AB ( ) where Ω is the model domain, AB is the mixing state, c AB is the ideal mixing when sub-scale fluctuations are not considered, and Σ AB is the contribution of the local fluctuations to the mixing state (which may be either positive or negative): c AB ( ) In the particular case A = B one recovers the classical definition.Additionally, we also quantify the relative deviation from the ideal well-mixed behavior: AB ( ) Here, quantity AB ( ) is analogous to the ( ) from de Dreuzy et al. [2012] for a single species.
In some cases with simple boundary and initial conditions, closed-form solutions exist for the integrals in (25).Below, we provide and discuss two such simple but representative cases.

Continuous injection
Consider a mean-uniform stationary flow in an infinitely long domain, which at the coarse scale can be considered as one-dimensional.At = 0, the concentrations of two solutes A and B are represented by Heaviside-step functions, forming a sharp interface at = 0: where H ( ) is the Heaviside step function.Additionally, AB ( , 0) = 0.The solution of the ADE (1) in this case is Then, the ideal mixing term is It can be shown (see derivation and generalized expression in Appendix D.2.1) that the departure of the actual mixing from (29), assuming a dual-rate (fast/slow) local mixing parametrization as in (12), is given by: AB ( ) where * is a dimensionless time and is the Dawson integral (D.15).
Figure 3 shows the evolution in time of the mixing metrics.Higher values of the slow-mixing fraction in the dual-rate model accentuate the departure of the actual mixing state (gray solid lines) from the ideal well-mixed case (dotted line).The relative difference between these two quantities, quantified by AB = 2 * AB (orange solid line), is highest at the beginning, and decays at late times ( 1) as the inverse of time.Hence at late times the actual mixing state AB converges towards the ideal well-mixed case c AB and scales with 1/2 .Looking closely at the log-log scale plot (Figure 3(b)), and taking the modeled mixing state AB as a proxy for the amount of reaction, we see that it does reproduce trends observed in mixing-limited systems such as simple Poiseuille flows [e.g., Perez et al., 2019, Figure 7].

Pulse injection
Now let us consider the same simple uniform flow in an infinite-length medium, but with a different initial condition.In this case, there is only one solute A, of which a mass (per cross- section unit area) o is injected over a small region around the origin with a Gaussian distribution characterized by a length o : Here we study the mixing state of A, AA ( ).Note that it has the opposite intuitive meaning than the AB ( ) analyzed in §3.2.1:A more advanced mixing process will be characterized by lower values of AA ( ), and viceversa.Once again, we assume that the initial condition for the fluctuations is AA ( , 0) = 0.
The resulting time-dependent mean-concentration profile is also a Gaussian: with o 2 o /2 .Then, the ideal mixing term is Once again, we characterize the relative deviation from the well-mixed behavior, which has the following closed form (see derivation and generalized expression in Appendix D.2.2) for the dual-rate (fast/slow) local mixing parametrization (12): AA ( ) Figure 4 depicts the evolution of the mixing metrics in dimensionless time, for a small value of initial pulse size, with o = (15 ) −1 .Similar to the case in §3.2.1, higher values of result in reduced mixing, as exhibited here by higher values of the actual mixing state (gray solid lines) compared to the ideal mixing state (dotted line).The ratio between these two quantities, AA = 2 * AA (orange solid line), is zero at = 0, since there is no incomplete mixing, and it starts to grow as the spreading process generates local concentration fluctuations.This increasing trend peaks at ≈ 0.785 (for the specific value of o = 1/15).After that, fluctuation destruction dominates and AA decreases as −1 for 1.At these long times, the actual mixing AA tends to approach the ∝ −1/2 trend of ideal mixing c AA .These features agree with semi-analytical [Bolster et al., 2011] and numerical [de Dreuzy et al., 2012] calculations of the mixing state evolution in fully-resolved porous media flows for a pulse injection of solute.

Reproducing results of reactive transport benchmark experiments
In this Section we use a numerical implementation of the proposed Lagrangian model to reproduce results from two well known laboratory experiments, conducted by Gramling et al. [2002], and Raje and Kapoor [2000].Hereafter, for conciseness, we refer to these two column experiments as G02 and RK00, respectively.

G02: Instantaneous equilibrium reaction
In the G02 experiments, performed in a column with a saturated granular material, a solution of EDTA 4− , initially occupying all the pore space with concentration o , was displaced longitudinally by an invading solution of CuSO 4 with the same molar concentration o .As these two solutes moved through the porous medium, the combination of hydrodynamic dispersion and molecular diffusion allowed them to mix and react forming CuEDTA 2− , among other reaction products.Hereafter, for simplicity, we refer to the three cited compounds as A, B, and C, respectively.The reaction can be expressed as with equilibrium equation, and a reaction rate that can be assumed instantaneous given the time scales of the experiment.Because equilibrium constant eq is very small (practically zero), A and B will always be instantaneously consumed when in contact, until one of them is exhausted locally (i.e., they cannot coexist).Hence, if we define the following conservative components, then the reaction product concentration will be given by The fully-resolved (pore-scale) transport of A and B follows the conservative form of the advectiondiffusion equation, with operator L defined as in (1), v being the heterogeneous velocity field within the saturated pore geometry, and being the molecular diffusion coefficient.The analogous of (41) applies to B ; however, in this particular case, because of the initial condition, A (x, ) + B (x, ) = o , hence we have that B (x, ) = o − A (x, ) and the transport is fully described by just one of the two equilibrium components.
However, in practice, a simple and complete solution is rarely obtainable, because of ( ) the lack of detailed information on the pore geometry and ( ) high computational demands, which is why this problem requires an upscaled approach.As an approximation, we ignore the boundary effect at the inlet, i.e., we assume an infinite medium.Then, the upscaled one-dimensional description of the transport of A and B , under the assumption of Fickian hydrodynamic dispersion, is identical to (28), where the constants and are the cross-section averaged vertical velocity and the upscaled longitudinal hydrodynamic dispersion coefficient, respectively.These constants were quantified by the authors of the cited experiment as = 1.21 × 10 −2 cm/s and = 1.75 × 10 −3 cm 2 /s.In the classical well-mixed upscaled ADRE approach, in which coarse-scale concentrations govern the chemical reactions, the combination of ( 42) with (40) would lead to the following equation for the concentration of C: However, the experimental observations of Gramling et al. [2002] do not agree with (43).Instead, the latter tends to overestimate the amount of reaction product generation, because of the incorrectness of the underlying assumption of full local mixing.
The proposed Lagrangian model is implemented as follows.Particles carry local concentrations of just one of the two conservative components, A, ( ), because the other is defined by B, ( ) = 1− A, ( ).Equal volumes (weights) are assigned to = 10 6 particles, which are initially distributed in space uniformly over an interval [− /2, /2], with = 15 cm, and A, (0) = H (− (0)). (44) As detailed in §2.1, transport and mixing of A, ( ) are decoupled and reproduced by equations ( 2) and ( 9), respectively.In the latter, the local averaging operator is implemented through binning (see Appendix B), with a bin size /300.We use a simple dual-rate mixing model like ( 12), parameterized by a slow mixing fraction and a slow mixing rate coefficient .Note that one does not need to explicitly simulate the evolution of the local concentration fraction corresponding to the fast mixing zone, A, ,1 ( ), which is always in equilibrium with the average at ( ).The coarse-scale reaction product concentration is given by the combination of ( 40) and ( 18), The code implementing this is written in Matlab (version 2016b) and the simulation with one million (10 6 ) particles runs in less than 5 minutes on a conventional laptop computer (Intel ® Core ™ i7-6700HQ, 2.60GHz).Equivalent, noisier results are generated in 2 seconds using ten thousand (10 4 ) particles.
An alternative approach to implement the proposed model is also tested, which does not require to explicitly simulate the Lagrangian particles.In this specific case, we have an analytical solution for A ( , ) = 1 − B ( , ), given by (42), as well as a semi-analytical solution for AB ( , ) = − AA ( , ) = − BB ( , ), given by (D.13).These quantities are, in fact, entries of the mean and the covariance matrix of a bivariate distribution (i.e., probability density function) of local concentrations of components A and B at any (coarse-scale) position and time, F ( A , B , , ).By assuming that F is multiGaussian, it is then fully defined by its mean and covariance matrix, and the local average (45) becomes where we use equation ( 42) for A and the numerical time-integration of (D.13) for AA = − AB .
Note that the multiGaussianity assumption may introduce inaccuracies, including the fact that a portion of F may fall outside the physically meaningful interval [0, o ].We refer to this solution strategy as the two-moment approach.9) and ( 12), with = 0.5 and = 10 −3 s −1 .The curves labeled as Lagrangian correspond to small-volume average performed on particles in a Lagrangian simulation, whereas the curves labeled as Two-moment are drawn using equations ( 46), ( 42), and (D.13).The shaded colored regions cover the range of product concentration results in the Lagrangian implementation for a ±15% variation of the and parameter values.The curves labeled as Well-mixed correspond to equation ( 43).
The dispersion coefficient is set to = 1.3 × 10 −3 cm 2 /s, somewhat lower than the value of = 1.75 × 10 −3 cm 2 /s estimated by Gramling et al. [2002] from the results of non-reactive experiments.This slightly lower than expected longitudinal dispersion for the reactive runs of G02 is common to several prior modeling interpretations [Rubio et al., 2008;Sanchez-Vila et al., 2010;Chiogna and Bellin, 2013].It may be related to several factors, including: ( ) the possible non-Fickianity of the curves from which the dispersivity value is adjusted [Gramling et al., 2002, Figure 4] combined with the intrinsic differences between the conservative and the reactive transport experiments, ( ) the fact that this value is calibrated for the reaction product only, or ( ) the unknown confidence intervals of the dispersivity value, presumably broad considering the reported variability between runs [Gramling et al., 2002, Table 1].
The results of the Lagrangian approach display close agreement with the experimental observations, as shown in Figure 5(a), for manually-adjusted values = 0.5 and = 10 −3 s −1 .A possible interpretation of the former is that in the pore-scale flow and transport conditions of the experiment, early-stage fast mixing controlled about half of the mixing process, whereas diffusive mixing across stable fluctuation structures was responsible for the other half.This is rationalized in Appendix A by idealizing local concentration fluctuations as periodic waves.
Assuming that the slow mixing is essentially two-dimensional (dominated by transverse diffusion between concentration filaments), and approximating both reactants' bulk diffusion coefficients as the value for C reported by the authors of the experiment, = 7.02 × 10 −7 cm/s, then according to (7), where = 0.13 cm is the mean grain size of the granular medium.That is, ℓ is probably close to the typical size of a pore, considering the reported porosity of 0.36.This suggests that the slow mixing process captured by the model corresponds indeed to the diffusive relaxation of pore-scale concentration fluctuations.The inferred value of = 0.5 shows that a single-rate local mixing model (i.e, = 1) would not be able to reproduce the experimental results.Neither would the high value of = 1 − / ≈ 1 that would render our model's local covariance behavior equivalent to Kapoor and Gelhar [1994] (see discussion below equation ( 23)).This is consistent with previous studies on mixing in porous media, both at pore and Darcy scales [e.g., Kapoor and Kitanidis, 1998;de Anna et al., 2014b;Le Borgne et al., 2013, 2015], which show that, after first encounter between two solutions with different composition, the mixing rate coefficient is higher at the beginning and decreases with time.In other words, the dynamics of mixing are subjected to aging, a feature which is effectively reproduced in our model by a parallel multi-rate process (9), without introducing any timedependent parameters.Although the dual-rate simplification appears to capture the general behavior for this case-study, allowing us to reproduce the experimental results, more complicated forms may be needed depending on the characteristics of the flow field, especially for highly heterogeneous porous media.This should be the subject of future research, as new experimental and/or numerical datasets in such settings become available.Looking closely at Figure 5, the main discrepancy between the data and the well-mixed solution ( 43) is the notable decrease in the peak of reaction product concentrations at the coarse-scale mixing interface.In the Lagrangian simulation, this reduction is caused by the anti-correlated fluctuations of A, and B, on particles with respect to the local average, which is also reflected by the negative values of AB , depicted in Figure 5(b).As expected, the spatiotemporal description of the covariance in the two-moment and in the Lagrangian approach are identical, which ratifies the validity of the expressions given in §3.1.However, the reaction product concentration prediction is slightly different, because of the multiGaussian approximation used in the two-moment approach.

RK00: Bilinear kinetic reaction
The laboratory experiments of RK00 consist of a similar setup to the one introduced in 4.1.In this case, the column is cylindrical, filled with spherical glass beads all of the same radius, and the reaction (aniline and NQS forming NQAB, among other reaction products) does not attain thermodynamic equilibrium instantaneously.Denoting once again the reactants as A and B, and the product as C, the reaction rate is bilinear with respect to the reactants' concentrations, with stoichiometric coefficients A = B = −1, AB = 1 in equation ( 15).The experimental observations are reaction product concentrations measured periodically at the outlet, which is located a distance = 18 cm from the inlet.
The authors conducted two experimental runs with different mean velocity and initial/invading reactant concentration o : Run 1, with = 0.096 cm/s and o = 0.5 mM; and Run 2, with = 0.07 cm/s and o = 0.25 mM.An additional difference between the two runs of RK00 is the switching of the roles of the A and B solutions (resident/invading), but this does not have any effect in the model if transport and mixing parameters are assumed to be species-independent.Longitudinal spreading is described by = , with dispersivity = 0.33 cm fitted from nonreactive transport experiments.
The evolution of C predicted by the model at = 18 cm is in reasonably good agreement with the experimental observations, as displayed in Figure 6(a), the main discrepancy being a delayed breakthrough of reaction product with respect to the experiments.Said discrepancy is observed even in the limiting well-mixed case, and therefore it may be due to non-Fickian spreading at early times.The errors in predicting the peak concentrations fall within a region of ±15% change in the mixing parameter values.Overall, we deem the model results as satisfactory considering that the parameter values have not been readjusted.The fact that G02 and RK00 exhibit consistent mixing behaviors when interpreted by the proposed model is promising in terms of future predictive use of the model.C at the outlet ( = 18 cm), to the experimental observations of RK00 in "Run 1" and "Run 2", and (b) corresponding point-covariance of A and B , AB .The local mixing model is given by ( 9) and ( 12), with = 0.5 and = 2.78 × 10 −3 s −1 .The curves labeled as Lagrangian correspond to small-volume average performed on particles in a Lagrangian simulation, whereas the curves labeled as Two-moment depict the results of a simultaneous first-order explicit finite difference solution of differential equations (D.26), (D.27), (D.28), (D.31), (D.32) and (D.33).The shaded colored regions cover the range of product concentration results in the Lagrangian implementation for a ±15% variation of the and parameter values.The curves labeled as Well-mixed correspond to the case = 0, solved either by the proposed particle method or by finite difference.
Like in the previous case, we compare the default Lagrangian model to a Eulerian (two-moment) implementation based on propagating both the mean and the covariance of the concentrations.The equation for the latter, however, is not (23), since in this case one should consider the effect of the kinetic reaction on the covariance.Approximate differential equations are derived in Appendix D.3 (by assuming that the local concentrations follow a multiGaussian distribution), and solved numerically with a first-order explicit finite difference scheme.Once again, the two-moment solution is different than its Lagrangian counterpart because of the multiGaussianity assumption, which exaggerates covariance removal due to chemical reaction (Figure 6(b)).This in turn results in a larger amount of reaction product generation in comparison to the Lagrangian implementation (Figure 6(a)).Nevertheless, both results share most key features and the use of this two-moment implementation approach may be advantageous in some situations, the main drawback being that the derived expressions (D.26), (D.27), (D.28), (D.31), (D.32) and (D.33) are only valid for the specific simple case of a bilinear kinetic reaction described by (48).

Summary and conclusions
We have proposed a Lagrangian mathematical model to represent the transport and mixing of solutes in a dual-scale (coarse/local) framework.Local concentrations carried by individual particles evolve by relaxation towards the coarse-scale concentration values that they perceive along their random path (described by ( 2)).This relaxation or mixing process is characterized by (9) as a parallel multi-rate interaction by exchange with the mean (MRIEM).We derived the differential equation describing the corresponding evolution of the (Eulerian) concentration point-covariance (22), and found solutions corresponding to the mixing state evolution for two simple generic cases.Finally, the proposed model (in its dual-rate form) was successfully implemented to reproduce reaction product concentration data from two well-known laboratory experiments that display incomplete mixing effects.Below, we enumerate additional findings and conclusions: 1.The partial differential equation describing the behavior of the local concentration covariance becomes nearly equivalent to the concentration variance conservation equation suggested by Kapoor and Gelhar [1994], given a dual-rate (fast/slow) parametrization of the mixing process.2. The temporal evolution of the mixing state for a pulse injection shows similar trends to those observed in previous studies of mixing in porous media within fully-resolved systems [Bolster et al., 2011;de Dreuzy et al., 2012], suggesting that the model may be able to accurately upscale local mixing limitations.Both for a pulse and for a continuous injection, the ratio between the mixing state components corresponding to the fluctuating and the averaged concentration terms decays at late times as the inverse of time.3. The analyzed experimental results would not be explicable, from our model's perspective, through a single-rate local mixing process.This agrees with previous knowledge on the complexity of mixing dynamics in porous media [Oates, 2007;Le Borgne et al., 2013;de Anna et al., 2014b], which establish the need to somehow include a temporal decrease of the mixing rate coefficient from the time of first coarse-scale contact between reactants.4. Along the same vein, the = 0.5 value for the dual-rate model that fits the experimental results does not agree with Kapoor and Gelhar's model, that is, with the restriction 2 = 1 − / ≈ 0.9995, if is assumed to be the molecular diffusion.This discrepancy can be attributed to the deformation-dominated stage where mixing occurs at a much faster rate than the stationary , as noted by Kapoor and Kitanidis [1998].5. Considering the characteristic local mixing distance ℓ = 4 / , the fitted value of = 10 −3 s −1 for G02 yields ℓ ≈ 0.05 cm, which is close to the typical pore size.This consistency of scale suggests that the model is properly capturing the physics underlying the mixing process.
Several future avenues for research arise due to this work, including ( ) identifying methods to readily estimate parameter values in ( 8) and ( 9) to upscale mixing within different local-scale heterogeneity patterns of velocity and dispersion, thus making the model scalable and translatable to the diverse range of hydrogeologic settings out there ( ) exploring the use of motion equations other than (2) [e.g.Berkowitz et al., 2006;Le Borgne et al., 2008], which are well-known to better describe coarse-scale transport in complex heterogeneous settings ( ) extending the model to incorporate heterogeneous reactions, and ( ) using the model to study the effects of local concentration fluctuations on realistic complex geochemical reaction systems.

A: Local-scale mixing event model
The model proposed in this work is based on fluid particles that carry local concentrations which may fluctuate with respect to the coarse-scale averaged concentration.These fluctuations are promoted by random walks of the particles, representing hydrodynamic dispersion.In that sense, each change in coarse-scale concentration experienced by an individual particle represents the start of a new mixing event, which consists in the relaxation of the marginal fluctuation through local mixing.The rate at which this relaxation occurs is well known to be subjected to aging, stabilizing after some initially faster, deformation-enhanced mixing [e.g., Kapoor and Kitanidis, 1998].Below we describe a simple mathematical model that would emulate such a marginal mixing event.
Spreading, simulated by the random motion of a particle, promotes marginal local departures from the local mean concentration.Before any mixing occurs, newly generated local fluctuations are perfectly sharp, hence they could be modeled as a square wave in the local space x = [ 1 , . . ., ] T with dimension .During the initial stage in which these fluctuation structures are being generated through fluid deformation, diffusion is enhanced by stretching [Le Borgne et al., 2013, 2015], and therefore the coefficient may appear larger than the actual diffusion.Since the rate of relative interface elongation decreases approximately as the inverse of time after the start of the mixing

B.1.1 Binning
A straightforward approach to compute the averaged concentrations is to discretize the spatial domain into a set of bins.The kernel (x − X ) in ( 3) and ( 4) is then replaced by an indicator function (x, X ) that has a value of 1 when x and X belong to the same bin, and a value of 0 otherwise.Compared to other smoothing techniques, this approach has very low computational demands.Another advantage is that it does not present any mass conservation issues.This is because the sum of all differences with respect to the mean within each individual bin is zero by definition, and therefore so is the sum of all exchanges with the mean given by ( 9).The main potential disadvantage of binning is that, compared to other smoothing techniques, it tends to require higher particle numbers (here, a finer Lagrangian discretization of the fluid mass) in order to converge to a smooth solution [Fernàndez-Garcia and Sanchez-Vila, 2011].This approach is used in the Lagrangian implementation described in §4.

B.1.2 Kernel smoothing
An alternative to binning is to use kernel smoothing on particles, that is, a radially symmetric kernel (x − X ; ℎ) with non-zero smoothing bandwidth ℎ in (3) and (4) to compensate for the finite particle density.This interpolation method is commonly used in smoothed particle hydrodynamics [Monaghan, 2005].This approach may offer better convergence rates with particle number than binning.On the other hand, in this case, exact mass conservation for (9) is not guaranteed by default, and a correction strategy is required.One mass-conserving approach can be obtained by considering symmetric pair-wise particle interactions.Let us first define a smooth interpolator to approximate (3) that is pair-wise symmetric Here, ˜ (x, X ) is some average of (x) and (X ).We redefine (9) by inserting A, , inside the local average operator The right-hand side of (B.2) clearly shows that a symmetric and therefore consistent mass exchange between each pair of particles , , is imposed [Herrera et al., 2009;Sole-Mari et al., 2019b].This expression may also be rewritten as B.3) which elucidates that mass conservation can be achieved in the kernel-based MRIEM through multiplication of A, , by a correcting factor 1 (X ), which converges to 1 as → ∞ and ℎ → 0. However, this approach does come at higher computational cost than binning.

B.2 Numerical dispersion and relation with other formulations
Smoothing tends to artificially spread out particle masses, which may generate some numerical dispersion in the Lagrangian numerical implementation.This can be straightforwardly quantified by comparing the right-hand side of (B.2) to the SPH formulation for simulating dispersion given a multiGaussian [Sole-Mari et al., 2019b, eq. 8].If we denote as SPH the dispersion coefficient in the cited SPH formulation, one can make both expressions equivalent by setting Since, ideally, the interaction given by (B.2) should not produce any dispersion (as explained in §2.1), the identity (B.4) quantifies the numerical dispersion involved in the numerical simulation of a MRIEM mixing with rate using a Gaussian smoothing kernel with bandwidth ℎ for computing the averaged concentrations.The same quadratic scaling should be expected for the bin size when binning is the chosen smoothing approach.Fortunately, as explained in §4.1, the high or virtually instantaneous fraction of mixing rates in the MRIEM formulation does not need to be explicitly simulated, and therefore only the small values of may produce numerical dispersion, which can be controlled by choosing a small-enough smoothing bandwidth ℎ, or by slightly reducing the value of used in the particle motion such that the added total dispersion has the correct value.The latter strategy, in fact, is tightly related to previous works [Benson and Bolster, 2016;Herrera et al., 2017;Sole-Mari et al., 2019b;Engdahl et al., 2019;Benson et al., 2019Benson et al., , 2020] ] in which the total dispersion results from the sum of ( ) random walks and ( ) some form of exchange between particles.In particular, if we look at the approach suggested by Benson and Bolster [2016, eq. 7], with mass transfers based on probabilities of collision between particles, we see that it is equivalent to the kernel form of the MRIEM as given by (B.2) for the specific case of a single mixing zone = 1 with mixing rate coefficient 1 = 1/d , and a multiGaussian with bandwidth where MT is the simulated dispersion.That is, the mass transfer approach by Benson and Bolster [2016] is equivalent to a Gaussian kernel-based single-rate IEM with instantaneous full mixing and "numerical" dispersion (insert (B.5) in (B.4)) SPH = MT .

C: Relation between local averages of Eulerian and Lagrangian derivatives
In this Appendix we provide the derivation for expression (10) given in §2.1.We start from the definition of local average (3).Note that, for incompressible flow, the fluid density (X ) is proportional to the porosity -amount of fluid per unit volume of medium.Hence, maintaining the assumption of constant porosity, the particle estimate of fluid density given by (4) must converge to a constant value as → ∞, being the number of particles.The time-derivative of expression (3) is then where the time-derivative of the Dirac delta function has been written as the variation Δ (x − X ) over Δ , where Δ → 0. A second-order Taylor expansion over the corresponding small particle displacement ΔX writes: where and are, respectively, the gradient and the Hessian matrix of .Then, using expression (C.2), knowing that the weighted summation is equivalent (in the limit → ∞) to an integral over the particle space, and given the evaluation properties of the Dirac delta distributional derivatives, we may rewrite the summation in (C.1) as where we applied well-known identities associated to (2) [Risken, 1989;Salamon et al., 2006], ΔX = vΔ and ΔX T ΔX = 2DΔ (for Δ → 0), and assumed that v and D are spatially constant.Finally, introducing the result from (C.3) in (C.1), we obtain the expression given by ( 10 The first-order Taylor expansion of Δ A, is We have obtained the partial differential equation describing the spatio-temporal evolution of the " " entry of the local concentration covariance of A and B in the absence of reactions.

D.2 Analytical expressions for the mixing state
Here, we derive the closed-form mixing state temporal evolution curves depicted in §3.2.For the derivations, we focus on the dual-rate case (12), which is more illustrative, and the results are then generalized to arbitrary MRIEM parametrizations.The differential equation ( 23) is linked with the solution of (1) through the source term (x, ) 2 2 ∇ T A D∇ B . (D.9) The other terms are exponential decay, advection, and dispersion.Therefore, if (x, ) is known, AB (x, 0) = 0, and the domain is unbounded, the solution for the point-covariance evolution can be obtained through the space-time convolution of the Greens function with the source term, i.where the operator | | applied to a tensor is its determinant.

D.2.1 Continuous injection
In the case introduced in §3.2.1, the source term of the covariance (eq.(D.9)) is:

D.2.2 Pulse injection
In the case introduced in §3.2.2, the source term of the variance (eq.(D.9) with B = A) is:

Figure 1 .
Figure 1.Local concentration variability within a true physical system and its conceptual representation in the proposed Lagrangian model, in which the coarse-scale concentrations are defined on the Eulerian space whereas the local concentrations are defined on the Lagrangian particles.Particles are represented by dark dots, and the colored circles around them show the corresponding local concentrations.

Figure 2 .
Figure 2. Spreading and mixing within a true physical system and their conceptual representation in the proposed Lagrangian model, in which the coarse-scale transport is simulated by the random motion of particles (2), and the local concentration is updated through multi-rate interaction by exchange with the mean (9).Particles are represented by dark dots, and the colored circles around them show the corresponding local concentrations or fluctuations.

Figure 3 .
Figure 3. Temporal evolution of various mixing metrics for the dual-rate model (12), according to the solution of equations (1) and (23), for two solutes A and B initially forming a sharp interface perpendicular to the flow direction, shown in (a) linear and (b) logarithmic scale.The normalizing factor ℓ 2 / is the typical distance traveled by the solute via dispersion within one characteristic mixing time.

Figure 4 .
Figure 4. Temporal evolution of various mixing metrics for the dual-rate model (12), according to the solution of equations (1) and (23), for one solute A initially placed as a Gaussian pulse of longitudinal standard deviation 0.26ℓ, shown in (a) linear and (b) logarithmic scale.The normalizing factor ℓ 2 / is the typical distance traveled by the solute via dispersion within one characteristic mixing time.
the derivative of the Dawson function (D.22).

Figure 5 .
Figure 5. (a) Comparison of various models' predictions of reaction product coarse-scale concentration, C , to the experimental observations of G02, at four different times, and (b) corresponding point-covariance of A and B , AB .The local mixing model is given by (9) and (12), with = 0.5 and = 10 −3 s −1 .The curves

Figure 6 .
Figure 6.(a) Comparison of various models' predictions of reaction product coarse-scale concentration,

Figure A. 1 .
Figure A.1.Illustrative depiction of the idealized local mixing event (A.1) for = 2 at different stages, with transitions controlled by each of the dual-rate mixing parameters and .Transition from (a) to (b) is nonlinear and fast, whereas transition from (b) to (c) is first-order and described by (A.6).Isolines split the images into eleven intervals of identical span.
, the resulting time-integral in (D.13) does not have an exact analytical solution.Nevertheless, a (pseudo-)closed form does exist for its integral in space (i.e., the local mixing term): 14) and (D.15), we see that there is an early-time regime where concentration covariance generation (promoted by spreading) dominates and Σ AB ∝ 1/2 , followed by a late-time regime where concentration covariance destruction (promoted by mixing) dominates and Σ AB ∝ −1/2 .Σ AB is negative meaning that the local covariance reduces the contact between A and B with respect to the ideal value c AB .The maximum negative magnitude of Σ AB ( ) is achieved for = 0.854 −1 .Dividing (D.14) by (29), we finally obtain the temporal description of the relative integrated departure from the well-mixed behavior in the continuous injection case, the beginning of D.2, the summation in (21) allows us to generalize the solution for any choice of mixing parameters as a summation of elementary building blocks: + o ) 3 e − ( − ) 2 2 ( + o ) , (D.19)and expression (D.10) gives the evolution of the variance:AA ( + o ) + ( − ) 2 ( + o ) (2 − + o ) 2.5 ( + o ) continuous injection case (D.13), we could not find a closed-form solution to the timeintegral in (D.20).But again, its spatial integral (the local mixing term) can be expressed in terms of the Dawson function (or more precisely, its derivative): 21) by (34), we obtain the temporal description of the relative integrated departure from the well-mixed behavior in the pulse injection case, we may define an elementary building block * AA to generalize (D.23) into more complex mixing parametrizations than the dual-rate form: