On the validity of effective formulations for transport through heterogeneous porous media

. Geological heterogeneity enhances spreading of solutes and causes transport to be anomalous (i.e., non-Fickian), with much less mixing than suggested by dispersion. This implies that modeling transport requires adopt-ing either stochastic approaches that model heterogeneity ex-plicitly or effective transport formulations that acknowledge the effects of heterogeneity. A number of such formulations have been developed and tested as upscaled representations of enhanced spreading. However, their ability to represent mixing has not been formally tested, which is required for proper reproduction of chemical reactions and which moti-vates our work. We propose that, for an effective transport formulation to be considered a valid representation of transport through heterogeneous porous media (HPM), it should honor mean advection, mixing and spreading. It should also be ﬂexible enough to be applicable to real problems. We test the capacity of the multi-rate mass transfer (MRMT) model to reproduce mixing observed in HPM, as represented by the classical multi-Gaussian log-permeability ﬁeld correlation heterogeneity


Introduction
Transport is anomalous in heterogeneous porous media.Anomalous transport observations include tailing in concentration breakthrough curves and plumes, or the strong increase in the rate of spreading of plumes.Several frameworks have been developed to generalize the advection-dispersion equation (ADE) and overcome its limitations (Frippiat and Holeyman, 2008).All these alternative frameworks share the goal to model complex permeability, velocity and concentration patterns in unified parsimonious effective equations.The limited number of parameters makes them efficient for the limited quantity of data usually available.In fact, they can be parameterized from breakthrough curves.They comply with the broad residence time distributions and non-local transport processes observed in reality (Gjetvaj et al., 2015;Le Borgne and Gouze, 2008;Willmann et al., 2008).They represent the consequences of complex concentration patterns, of simultaneous concentration trapping and fast progress on residence times while averaging out all the fine concentration structures in the upscaling process.These anomalous transport frameworks have proven to be highly effective for residence times, transport time distribution and effective spreading, both phenomenologically and practically (Berkowitz et al., 2006;Neuman and Tartakovsky, 2009).However, their ability to reproduce mixing, which is required for properly reproducing chemical reactions, has not been tested.

Published by Copernicus Publications on behalf of the European Geosciences Union.
We argue that an effective transport formulation should honor not only the mean advection, and spreading observed in heterogeneous porous media (HPM), but also the evolution of mixing.This should not be understood as limiting anomalous transport frameworks but as extending them to handle broader ranges of physical and chemical processes, and at further promoting the approach of effective equations that upscale out the fine-scale structures to retain only their main consequences in terms of transport, reactivity and reactive transport couplings.
Here, we investigate the relevance of multi-rate mass transfer (MRMT) framework to model not only spreading but also mixing.MRMT is taken as a typical anomalous transport framework.Its advantage lies in providing local concentrations, which can be straightforwardly used to evaluate concentration variance, mixing and mixing-induced reactivity (Babey et al., 2014;Carrera et al., 1998;de Dreuzy et al., 2013;Haggerty and Gorelick, 1995), as well as the apparent reduction in the rate of kinetic reactions (Dentz et al., 2011).The question is whether its validity as a representation of transport through HPM can be extended to reproduce the effects of the evolution of mixing rates resulting from the stretching and folding associated with complex velocity structures (de Anna et al., 2014b;Jimenez-Martinez et al., 2015;Le Borgne et al., 2015).
Some assessment of MRMT to model reactivity in HPM has been made in former works (Willmann et al., 2010).Equivalent reactivity has been evaluated at some welldefined travel distances on MRMT calibrated on residence time distributions.Here, we follow a different approach by analyzing the temporal development of spreading and mixing.We extend the integrated assessment of mixing-induced reactivity at given travel distances to its temporal development.
Our contribution concerns the comparison of different models much more than the HPM and MRMT models themselves.For the sake of completeness, we recall model equations and simulation methods in Sect. 2 (models and methods) and measures of spreading and mixing in Sect.3. We use these measures to propose the conditions that should be met by effective (upscaled) transport formulations to be considered valid representations of transport through HPM (Sect.4).We then test whether MRMT formulations meet the proposed conditions (Sect.5).While this last section depends on the specific choice of the MRMT framework as an equivalent transport model, the comparison methodology is independent of it and can be used to assess transport equations respecting both spreading and mixing.

Model and methods
We present the MRMT and HPM models sequentially.As they are both well known, we present only the main equations and highlight the critical assumptions of importance in this study.

Multi-rate mass transfer model (MRMT)
MRMT models express anomalous transport by the interaction between transport in a mobile zone and a series of immobile zones (Carrera et al., 1998;Haggerty and Gorelick, 1995).Transport in the mobile zone is advective and dispersive with a mean solute velocity v (water flux q divided by mobile porosity, ϕ) and a dispersion coefficient d.Each immobile zone i is parameterized by a characteristic rate α i (inverse of a characteristic exchange time) and an immobile porosity ϕ i .The concentrations c and c i (i = 1. ..N ) in the mobile and immobile zones, respectively, are determined by the following set of equations: The ratio of immobile to mobile water volumes is rated by the total capacity ratio β = ϕ i ϕ .The term capacity derives from the fact that MRMT formulations were originally devised to represent trapping by sorption in hard-to-reach sorption sites, which were characterized by capacity (including both dissolved and sorbed solute mass) (see, e.g., Haggerty and Gorelick, 1995).We use here an equivalent MRMT formulation for non-sorbing solutes, so as to facilitate comparison with HPM.Initial and boundary conditions will be described later for both MRMT and HPM models.MRMT models differ by the distributions of characteristic rates α i and immobile porosities ϕ i .Among the available models (Cvetkovic, 2012;Haggerty et al., 2000), we choose a uniform distribution for characteristic times (1/α i ) bounded by the two extreme rates α 1 = 1/t 1 and α N = 1/t N (t 1 < t N ) and a power-law distribution for ϕ i : ( The power-law distribution is consistent with observed breakthrough curves in HPM, which often display long tails that appear linear in log(c) versus log(t) (Gouze et al., 2008;Haggerty et al., 2004;Li et al., 2011;Silva et al., 2009;Willmann et al., 2008).This tailing is well modeled by a power law, such that the breakthrough concentration c evolves as c ∼ t −m .Haggerty et al. (2000) showed that the slope m relates to the exponent of the power-law distribution of the MRMT rates (Eq.3).m is generally found to be in the interval [1.5, 2.5] but little is known about its relationship to the geological heterogeneity.Willmann et al. (2008) found some correlation between the degree of connectivity and the slope.The more connected the field, the smaller the slope.
In this context, fracture-matrix exchanges in fractured media represent the lowest bound (m = 1.5), which is controlled by diffusion into immobile regions (Haggerty and Gorelick, 1995).On the contrary, a slope m of 2.5 may represent a heterogeneous but poorly connected hydraulic conductivity field, where late time arrival is controlled by slow advection.We simulate MRMT models with a standard time-and space-adaptative method that preserves mass (de Dreuzy et al., 2013) and always complies with the CFL conditions (Daus et al., 1985).The advective and the diffusive processes in the mobile zone, as well as the exchange with the immobile zones, are treated with a sequential non-iterative coupling method.These methods lead to efficient simulations of large spatial domains and extended times with initial refined resolutions.We have successfully compared them with a more classical fixed-time Galerkin finite element method, integrated with the fourth order Runge-Kutta method (ode45 function of Matlab) and found relative differences less than 10 −3 %.Simulations have been performed over the time required for transport to reach its asymptotic regime.

Heterogeneous porous media (HPM)
For reference purposes, we restrict the analysis to heterogeneity of hydraulic conductivity (K) as represented by the classical 2-D Gaussian correlated multi-Gaussian logK fields.These are characterized by their isotropic correlation function: with r the distance, λ the correlation length which is used to scale distances, and σ 2 Y the variance of the logarithm of Y = logK.We use simulation results performed in previous studies (de Dreuzy et al., 2012) obtained on 2-D domains of sizes L L and L T in parallel and orthogonal directions, respectively, to the mean flux.L L is large enough to avoid any finite-size effects (from 10 2 to 10 3 correlation lengths λ).Boundary conditions for flow and transport are periodic in the transverse direction to minimize boundary effects.L T is of the order of 100 times λ to ensure initially ergodic transport conditions.Under such uniform extended injection con-ditions, transport in HPM can be considered ergodic and can be fundamentally compared with a 1-D MRMT model.The immobile zones of MRMT can be viewed as representing the low velocity zones of HPM, so that the mobile zone may represent the high-velocity channels.
Flow is solved with a finite volume scheme with permeameter-like boundary conditions under a unit head gradient.Transport is simulated using the ADE, with heterogeneous advection and homogeneous diffusion.Therefore, it is characterized by the Peclet number Pe, equal to the mean velocity multiplied by the correlation length, divided by the diffusion coefficient.Transport is simulated with a random walk Lagrangian method.Numerical methods are exhaustively described in several previous papers (Beaudoin et al., 2006(Beaudoin et al., , 2007(Beaudoin et al., , 2011)).

Injection and boundary conditions
The same type of injection and boundary conditions are used for both models.Flow has a major flow direction imposed in HPM by a head gradient in the longitudinal direction and periodic boundary direction in the transverse direction.For transport, reflecting and absorbing boundary conditions are used upstream and downstream, respectively (Beaudoin and de Dreuzy, 2013).Injection is performed downstream to the inlet boundary to minimize boundary effects.
Extended injection conditions are used for the HPM and MRMT models.Concentrations are homogeneous orthogonally to the main flow direction within a square wave of longitudinal and transverse widths L 0 and T 0 , respectively.In the HPM case, concentration is a sole function of the coordinate x L along the flow direction: with c 0 given by: ϕ T is the total porosity.To ensure that the same mass m 0 is injected in the HPM and MRMT cases, we adapt the initial state of the MRMT model to: Spreading becomes independent of the injection length when the longitudinal plume size becomes significantly larger than L 0 .Mixing depends on the injection conditions more critically than spreading, as the initial concentration value depends on the injection width L 0 (Eq.6).
Concentration fields normalized by their maximal value c(x,t) / max(c(x,t)) and their related Gaussian profile concentrations c D (x,t) / max(c(x,t)) in the bar over them at the four evolving times indicated in Fig. 2. In this case, the time at which the non-dispersive mixing reaches its maximum t γ max is of the same order as the advection time.
3 Measures of spreading and mixing

Spreading
For an extended plume, spreading is generally measured by the square root of the second centered moment of the spatial distribution of concentration σ L : where m (k) L (t) is the kth order moment of the concentration distribution with x L the coordinate of x in the direction parallel to the main flow direction (longitudinal direction) and the flow domain.With this definition, σ L can be viewed as the longitudinal extent of the plume (i.e., how far it spreads).Dispersion is the rate of spreading (i.e., time derivative of σ 2 L ), usually characterized by the longitudinal dispersivity α L : where v is the plume velocity equal to the time derivative of the mean position plume m (1) L (t). α L increases until it converges to an asymptotic value α LA , thus defining in turn the asymptotic regime (Dagan, 1990;Gelhar, 1993).
In MRMT, spreading comes from the exchanges to the mobile zone.That is, spreading results from trapping.Solutes are slowed down and dispersed by the exchanges with the immobile zones.The resulting dispersivity is a monotonously increasing function of the residence times in immobile zones (both their mean τ MRMT and range (t N − t 1 )).The dispersivity induced by the dispersive and diffusive processes in the mobile zone is comparatively negligible and could be disregarded.
In HPM, spreading comes both from diffusive exchanges with low velocity zones and from spatial fluctuations of the velocity field (de Dreuzy et al., 2007;Salandin and Fiorotto, 1998).The asymptotic dispersivity increases both with the correlation length λ and with the logK variance σ 2 Y : where g is either a linear function for small values of σ 2 Y (σ 2 Y < 1) or a quadratic function at larger values (de Dreuzy et al., 2007).
h σ 2 Y , Pe is a correction factor accounting for diffusion (Beaudoin et al., 2010).Local diffusion reduces the effective dispersivity in the high-heterogeneity cases by releasing solutes from the low velocity zone and truncating the trapping times induced by slow advection.
Any concentration plume can be approximated by a Gaussian concentration profile c D (x,t), defined by the two first moments, m (1) L (t) as the mean and σ 2 L (t) as variance.It is the smoothest equivalent profile.Both MRMT and HPM converge asymptotically to this profile.However, it is far away from the full concentration profile c(x,t) at any time, as shown by the comparison of Fig. 1.At early times (left snapshots in Fig. 1), the concentration profile remains heterogeneous especially in the transverse direction with both higher and lower concentrations.Around the advection time, defined as the correlation length λ divided by the plume velocity v, the deviation reaches its maximum.At this point, the Gaussian concentration profile has become much more diluted than the real concentration field (second from the left snapshot of Fig. 1).Concentration inhomogeneities decrease very slowly and remain over very long periods of time even though the range of concentration values decreases (two right-most snapshots of Fig. 1) (de Anna et al., 2014a;Jimenez-Martinez et al., 2015;Le Borgne et al., 2011).
In summary, in HPM, dispersivity comes primarily form the velocity structure, which drives the generation of gradients in concentration, and thus, mixing.Instead, in MRMT, effective dispersivity is controlled by mobile-immobile exchanges and delays the actual mixing between the immobile and mobile solute concentrations.

Mixing
The Gaussian profile only gives a crude approximation of the concentration field with a strong deviation on the distribution of concentration values, especially at early times when diffusion has not homogenized the concentration field in the transverse direction (Fig. 1).Actual concentrations remain much higher and closer to the initial concentration value than in the Gaussian profile prediction.That is, the initial concentrations are much less diluted (i.e., mixed) than in the maximum entropy Gaussian distribution.The Gaussian profile c D (x,t) thus sets a lower bound to the effective concentration variability.Therefore, it is most natural to compare the actual distribution of concentration values to that of the Gaussian profile in order to describe the mixing state.Notice that, contrary to spreading, we are not concerned here with the spatial distribution, but only with the values of concentration and their time evolution, which are most simply characterized by the second moment.We quantified the deviation from the Gaussian mixing regime as the ratio of the actual concentration second moment M(t) to the second moment M D (t) of the Gaussian profile concentration c D (x,t) minus 1 (de Dreuzy et al., 2012): with and the second moment of the reference Gaussian concentration: M D is directly the square of the injected mass m 2 0 , divided by an effective area occupied by the plume 2 √ π T 0 σ L .As M(t) is always larger than M D (t), γ is always positive.γ is initially and asymptotically very close to 0. It is however significantly positive while the concentration distribution is far from the Gaussian profile.M(t), introduced here as a measure of global concentration variability, is widely used because its derivative is giving the dissipation rate and determines the physical constrains of chemical reactivity (de Simoni et al., 2005;Le Borgne et al., 2010).The dissipation is also closely related to the dilution index, which is another measure of mixing (Kitanidis, 1994;Rolle et al., 2009).It should be finally noticed that γ and M D fully characterize the mixing state given by M: In HPM models, resistance to dispersive mixing, as we can also call γ , is enhanced by heterogeneity and reduced by large diffusion rates (smaller Peclet number) (de Dreuzy et al., 2012).γ sharply increases at initial times to a maximum value γ max , at a time t γ max close to the advection time, and slowly decreases back to 0 (Fig. 2).The time range, over which γ is significantly non-zero, can be characterized by r tγ , which is the ratio of the upper and lower times at which γ is equal to a quarter of its maximal value γ max .While the amplitude of γ depends on the variability of the velocities and on the rate of advection to diffusion, the shape of the function γ remains unchanged by the K field heterogeneity (σ 2 Y ), the ratio of advection to diffusion (Pe) and the width of the initial conditions ( L 0 ).The time range r tγ , over which γ is non-negligible, also remains constant (Fig. 2).Therefore, t γ max can be used for scaling time, so that γ can be written as: where f is the characteristic scaling function (Fig. 2, insert).A similar constant shape behavior has been noted for viscous fingering in heterogeneous velocity fields (Jha et al., 2011a, b).

Conditions for effective formalisms of transport through HPM
We propose four conditions for any effective transport formulation to be considered as a valid representation of transport through heterogeneous media.In essence, an effective transport equation should yield the same mean advection, spreading and mixing as the HPM and be sufficiently flexible to represent real problems.Evaluation of these conditions can be done as follows: 1. Mean advection simply requires mean water velocity (i.e., mean plume velocity for non-reactive solutes) to equal v = q/ϕ T .This condition can be met by all published upscaled transport equations, by imposing some simple constraints on their parameters.In MRMT, it is sufficient to impose ϕ T = ϕ + ϕ i = ϕ (1 + β).
2. Spreading is characterized by dispersivity, which measures the rate of growth of plume size (Eq.10).In cases   12) in HPM for evolving logK variances, σ 2 Y , under a small-width injection window ( L 0 / α LA = 0.075), flux-weighted injection conditions and Pe = 100 (adapted from de Dreuzy et al., 2012).The similarity of function shapes is highlighted in the insert by the scaling function f of Eq. ( 16) where the thick black line is the average of the displayed functions.Note that the time of maximum deviation, t γ max , is hardly affected by σ 2 Y and falls around the characteristic advection time λ/v.The four dashed lines indicate the times displayed in Fig. 1.
where asymptotic dispersion is reached, this condition implies that dispersivity of the effective equation should tend to the asymptotic dispersivity of the HPM.Otherwise, dispersion (or directly, spread, as measured by σ L ) can be compared to a spatial scale comparable to the problem dimension (e.g., size of the aquifer or distance covered by the plume).
In addition, the time required to reach the above dispersion value should also be honored by the effective formulation to ensure that the rate of growth of the plume is reproduced.In our case, where asymptotic dispersion is reached, we propose to define this criterion in terms of r α , the mean distance covered by the plume at the time t α LA /2 where dispersivity reaches half of its asymptotic value normalized by the asymptotic dispersivity α LA : where t α LA /2 is implicitly defined by r α can also be interpreted as the ratio of advective and dispersive scales like in the definition of the Peclet number.
3. Mixing is required for properly reproducing fast reactions (slow reactions should be properly reproduced if the resident time distribution is honored, which is assured if mean advection and dispersion are reproduced).
As discussed above, mixing is essentially dispersive and well characterized by M D (Eq.14) for late times.Therefore, assuming dispersion to be well reproduced, an effective transport formulation only needs to reproduce the deviation from dispersive mixing, characterized by γ (Eq.12).In the first stage, the comparison can be restricted to the amplitude of the deviation γ max and the time range over which it extends r tγ .In a more advanced stage, the characteristic shape of the γ function, f , can be used for comparison.
To compare the timings of spreading and mixing, we define the additional criterion r MT as the ratio of the characteristic spreading time t α LA /2 to the characteristic mixing time t γ max r MT of the development of the resistance to mixing and rates the lag between the timing of mixing and spreading.
4. Most of the work on effective transport is of a theoretical nature, but the ultimate goal should be applica-tion to real problems.This implies that a valid transport formulation should be able to accommodate different types of boundary conditions and flow regimes (i.e., transient flow) and dimensions.Most importantly, it should accommodate characterization.Dispersion usually includes the effects of heterogeneity and uncertainty.Whereas the latter is reduced by aquifer characterization, the former is not.Specifically, hydrologists use geology, hydraulics, geophysics, hydrochemistry and isotopes to figure out, among other things, the patterns of spatial variability of hydraulic conductivity.
The resulting models display variability not only in the mean logK but also on their correlation distance and variance.An effective transport formulation should be able to honor this variability.

Results and discussion
We consider it well established that MRMT, and other nonlocal in-time formulations, can reproduce mean advection and spreading, as discussed in the introduction.Mean advection in the MRMT approach is equivalent to that of the HPM provided that flux and total porosity are equivalent.
Additionally, the distribution of residence times in immobile zones can be adapted so that the asymptotic dispersivity of the MRMT model be equal to that of the HPM model in Eq. ( 11).It is always possible, as dispersivity is an increasing function of the residence times.This imposes a condition on the temporal range of t 1 , t N or equivalently on their mean residence time τ MRMT .As trapping in the immobile zones is the main dispersive mechanism, the mean residence time is logically adapted to calibrate the asymptotic dispersivity.With the total flow imposed to be set by the HPM, the characteristic spatial scale is the typical plume position at τ MRMT .As the characteristic spatial and temporal scales are interrelated to ensure consistent asymptotic behaviors, comparison of results can be performed on dimensionless terms and should ensure consistent preasymptotic regimes.In fact, MRMT models are calibrated on tracer tests and breakthrough information, but this does not ensure a good reproduction of mixing (Luo and Cirpka, 2011).Therefore, we restrict our comparison to mixing criteria and sensitivity to initial conditions.

Comparison of mixing in HPM and MRMT
In HPM, the temporal extension of the deviation from the dispersive mixing regime r tγ does not depend significantly on the permeability heterogeneity, as also expressed by the constancy of the shape of γ (Fig. 2).We thus compare the shape of γ obtained for the HPM with σ 2 Y = 9 (f function of Eq. 16) to shapes of γ obtained for various MRMT models obtained under consistent injection conditions (Eqs.6 and 7).For MRMT, extreme values have been investigated to get the possible range of behaviors.For slopes m, we adopted the range observed in nature, as discussed in Sect.2.1, with m varying between m = 1.5 (typical fracture-matrix case) and m = 2.5.The single rate mass transfer (SRMT) is also shown for comparison.The porosity ratio β does not have an upper bound.In fact, ideally, the mobile porosity could be 0. We adopted β = 150 as a large upper value.Larger upper values would not affect results and might cause numerical difficulties.The same can be said for t N /t 1 , for which we took t N /t 1 = 10 3 as the upper bound.The analysis presented hereafter has been made for different combinations of parameters within these bounds.As all models lead to consistent conclusions, we only present the most characteristic results.
All MRMT models capture the sharp rise of γ at times smaller than t γ max (Fig. 3).The sharp rise comes from a strong initial divergence from the equivalent Gaussian concentration profile.Initial behavior is dominated by the contrast of the quickly progressing concentrations in the mobile zones and trailing concentrations in the low-flow or immobile zones.On one hand, dispersion induces a sharp decrease of M D , which is inversely proportional to σ L (Eq. 14).On the other hand, trapping maintains high concentrations in the immobile zone and high values for the second moment of the concentration distribution M. Divergence of M from M D increases until it reaches its maximum at t γ max .
At larger times, progressive release of solute mass from the immobile zone and equilibration with the concentration values in the mobile zone allow f to decrease.The insert of Fig. 3 highlights, in a lin-lin graph, the differences of the decrease stage.The MRMT model that best matches the scaling function f is obtained for m = 2.25.For power-law slopes m larger than 2, the decrease of f is qualitatively similar to that of the single-rate mass transfer model.For m values closer to 2, the deviation is more sustained but displays the same decreasing trend.The behavior changes significantly when m becomes smaller than 2, with a very slow decrease of f coming directly from the effect of the slow drainage of the immobile zones with the smallest rates (long exchange times).At least in MRMT models with m lower than 2, trapping displays a much longer memory effect in MRMT than in HPM.
While they are similar to HPM for the extension of the non-dispersive mixing regime, MRMT models with slopes m larger than 2 converge to less anomalous single rate mass transfer.In terms of mixing, this translates in small amplitudes in γ (small γ max values).For m = 2.25 and L 0 / α LA = 0.075, we have computed γ max for a large variety of t N /t 1 and β values, the only two remaining parameters of the MRMT model.γ max first increases with t N /t 1 and β and quickly saturates to a maximal value of 2.85 (Table 1).γ max varies between 0.57 and 2.85 by a maximum factor of 5.In the HPM case, however, γ max is always larger than 3, reaches values of 15 for a Peclet number Pe of 100, scales like the square root of Pe and is thus not limited in amplitude.MRMT models cannot match both the amplitude and the timing of γ .For MRMT models with m-slopes larger than 2, mixing is far more dispersive in MRMT than in HPM (smaller deviation values γ in the MRMT models).MRMT models with m-slopes smaller than 2 induce larger but much too sustainable deviations.As a result, MRMT models have a stronger memory of trapping or display less non-dispersive mixing, without excluding to display both differences simultaneously.The difficulties of MRMT models to capture mixing might be linked to the existence of the structure of concentrations in lamellas where stretching and folding extends the concentration forward and enhances the eventual mixing by diffusion (Fig. 1) (de Anna et al., 2014b;Le Borgne et al., 2015).
Concerning the non-dispersive mixing shapes of the scaled function f , MRMT models display a much broader range of shapes than HPM.The invariant shape property of HPM is not recovered in MRMT.On the contrary, MRMT shapes depend strongly on the distribution of transfer rates.This is an advantage to match a wider range of cases issuing possibly different γ functions.However, it is a drawback to fit just one case, as it restricts the MRMT models that can match HPM simulations with broad ranges of Pe and σ 2 Y values.

Influence of initial injection size
To qualify the memory effect in MRMT and HPM, we analyze their sensitivity to the initial injection width L 0 .Spreading, as defined by the characteristic longitudinal plume extension (Eqs.8, 9 and 10), does not depend on the initial concentration c 0 (Eq.6).σ L is initially of the order of L 0 and quickly becomes larger.Spreading quickly loses any memory of the initial conditions.We note that this is the case because sampling effects do not intervene as the transverse injection scale is assumed large enough to ensure ergodic sampling by itself.Like spreading, the second moment of the reference Gaussian concentration M D (t) does not depend on the initial concentration but only on the in- jected mass divided by the characteristic area occupied by the plume ( T 0 σ L ) (Eq. 14).
The concentration second moment M(t) (Eq.13), however, depends critically on the injection width through the relation between injected mass and concentration (Eq.6).At initial times, the concentration second moment is proportional to the injected concentration value.At later times, the concentration second moment has lost the memory of the initial concentration and is only the function of the injected mass m 0 .As a result, we expect that the concentration sec-ond moment and the deviation towards the dispersive mixing regime γ depends on the injection conditions, here represented by the injection width L 0 .On the basis of numerical simulations, we compare the evolution of γ with L 0 for HPM and MRMT models.
Results of the γ function for both MRMT and HPM models are displayed in Fig. 4a and b, respectively, for different injection sizes.We have performed simulations for comparable ranges of L 0 / α LA values (0.005-0.1 for HPM and 0.02-0.4for MRMT).We have checked numerically that the results for the two specific MRMT and HPM cases display generic tendencies.In both models, injection width has a critical influence on γ (Fig. 4).Smaller injection windows let the initial concentration increase and enhance the deviation towards the dispersive mixing regime.We use the maximum deviation γ max to characterize the overall influence of L 0 .In the MRMT and HPM models, maximum deviations γ max have different scaling (Fig. 4a and b, inserts): For MRMT, γ max evolves like the initial concentration level (Eqs.6 and 7).For HPM, γ max has the same scaling for the initial conditions γ max ∼ 1 √ L 0 and for the diffusion coefficient as γ max ∼ √ Pe (de Dreuzy et al., 2012).Doubling the injection width has a comparable effect to doubling the diffusion coefficient.Initial dilution over L 0 and dilution induced by diffusive-dispersive processes reduce the overall deviation to the dispersive mixing regime γ in the same proportion.The reduction of concentration in HPM comes from the diffusive-dispersive processes while, in MRMT, it comes from the progressive release of solutes with high concentrations close to c 0 trapped in the immobile zone.Due to their differing signatures, both processes cannot be compared and the dispersive-diffusive processes of HPM cannot be modeled as trapping-release mechanisms.

Conclusion
We propose conditions to test anomalous transport frameworks not only on spreading but also on mixing.We define a minimum set of six essential constraints that should be respected in order to retain the main transport, reactivity and reactive transport couplings.These constraints involve the conservation of (1) the mean advection, (2) dispersivity amplitude and (3) timing generally imposed.Beyond these flow and spreading metrics, (4) amplitude and (5) timing of the deviation towards the dispersive mixing regime should be respected.The last condition concerns (6) the respective timings for mixing and spreading.Under ergodic injection conditions, spreading is characterized by the standard dispersivity describing the evolution of the plume size along the main flow direction.Mixing is characterized by the deviation from the dispersive mixing regime γ , defined as the second moment of the concentration distribution of a conservative tracer, divided by the second moment of a Gaussian concentration pattern with the same spread minus 1 (Eq.12).Zero initially and asymptotically, γ traduces the macroscopic effect on mixing of the concentration structures within the solute plume.
We use these criteria to evaluate MRMT models by comparison to advective-diffusive transport simulations through HPM, represented by the classical isotropic 2-D Gaussian correlated multi-Gaussian log-permeability fields, characterized by variances between 1 and 9.A broad range of MRMT models are considered.We conclude that MRMT models cannot match both the amplitude and the timing of γ .MRMT models can reproduce observed spreading rates and some non-dispersive mixing.However, they tend to induce larger and too sustained deviations from dispersive mixing.As a result, MRMT models display a longer memory but less nondispersive mixing than HPM.We attribute this divergence to the fact that MRMT represents non-dispersive mixing through trapping mechanisms, whereas non-dispersive mixing is controlled by stretching and folding in HPM.Divergent sensitivities to initial conditions confirm that dispersivediffusive induced mixing in HPM cannot be modeled by mobile-immobile models.
Our study does not preclude, however, the existence of effective transport equations consistent with spreading and mixing of HPM.Nonetheless, we argue that the proposed criteria and existing results of HPM should be used as guidelines to set up effective transport equations that respect spreading, mixing and eventually reactive transport.

Figure 2 .
Figure 2. Time evolution of the deviation from dispersive mixing γ (t) defined by Eq. (12) in HPM for evolving logK variances, σ 2 Y , under a small-width injection window ( L 0 / α LA = 0.075), flux-weighted injection conditions and Pe = 100 (adapted from deDreuzy et al., 2012).The similarity of function shapes is highlighted in the insert by the scaling function f of Eq. (16) where the thick black line is the average of the displayed functions.Note that the time of maximum deviation, t γ max , is hardly affected by σ 2 Y and falls around the characteristic advection time λ/v.The four dashed lines indicate the times displayed in Fig.1.

Figure 4 .
Figure 4. Dependence of the deviation from dispersive function γ on the injection width L 0 /α LA for (a) MRMT and (b) HPM models.Insert shows the dependence of γ max on L 0 /α LA .

Table 1 .
Comparison of γ shapes (f -scaling functions defined in Eq. 16) for HPM and MRMT simulations with slopes of the power-law distributions of MRMT transfer rates, m, between 1.575 and 2.5.HPM is represented by the broad black curve obtained from the insert of Fig.2.The insert displays the same curves in arithmetic scale.Values of the maximum deviation to the dispersive mixing regime γ max for MRMT models with a power-law exponent of the rate distribution m = 2.25 and L 0 / α LA = 0.075.