Joaquim Soler-Sagarra

Vivien Hakoun
email: vivien.hakoun@gmail.com

Marco Dentz
email: marco.dentz@csic.es

Jesus Carrera
email: jesus.carrera.ramirez@gmail.com

The Multi-Advective Water Mixing Approach for Transport through Heterogeneous Media

Keywords: MAWMA, mixing, heterogeneity

Finding a numerical method to model solute transport in porous media with high heterogeneity is crucial, especially when chemical reactions are involved. The phase space formulation termed the multi-advective water mixing approach (MAWMA) was proposed to address this issue. The water parcel method (WP) may be obtained by discretizing MAWMA in space, time, and velocity. WP needs two transition matrices of velocity to reproduce advection (Markovian in space) and mixing (Markovian in time), separately. The matrices express the transition probability of water instead of individual solute concentration. This entails a change in concept, since the entire transport phenomenon is defined by the water phase. Concentration is reduced to a chemical attribute. The water transition matrix is obtained and is demonstrated to be constant in time. Moreover, the WP method is compared with the classic random walk method (RW) in a high heterogeneous domain. Results show that the WP adequately reproduces advection and dispersion, but overestimates mixing because mixing is a sub-velocity phase process. The WP method must, therefore, be extended to take into account incomplete mixing within velocity classes.

Introduction

The advection-dispersion equation (ADE) is the most widely used formulation to model solute transport through porous media. Yet, it does not adequately characterize transport in heterogeneous media [START_REF] Gjetvaj | Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone[END_REF][START_REF] Le Borgne | Non-Fickian dispersion in porous media: 2. Model validation from measurements at different scales[END_REF][START_REF] Willmann | Transport upscaling in heterogeneous aquifers: What physical parameters control memory functions?[END_REF] where dispersion grows with scale [START_REF] Gelhar | A Review of Field-Scale Physical Solute Transport Processes in Saturated and Unsaturated Porous Media[END_REF][START_REF] Neuman | Universal scaling of hydraulic conductivities and dispersivities in geologic media[END_REF], nonequilibrium occurs [START_REF] Vogel | The dominant role of structure for solute transport in soil: Experimental evidence and modelling of structure and transport in a field experiment[END_REF][START_REF] Alcolea | Regularized pilot points method for reproducing the effect of small scale variability: Application to simulations of contaminant transport[END_REF], or breakthrough curves display tailing [START_REF] Valocchi | Validity of the Local Equilibrium Assumption for Modeling Sorbing Solute Transport Through Homogeneous Soils[END_REF][START_REF] Carrera | An overview of uncertainties in modelling groundwater solute transport[END_REF]. These features are not well represented by the ADE. Therefore, transport through heterogeneous media is called anomalous (i.e., non-Fickian transport). The problem becomes critical when chemical reactions are involved [START_REF] Soler-Sagarra | Simulation of chemical reaction localization using a multi-porosity reactive transport approach[END_REF][START_REF] Battiato | On breakdown of macroscopic models of mixing-controlled heterogeneous reactions in porous media[END_REF][START_REF] Sadhukhan | A simulation study of reactive flow in 2-D involving dissolution and precipitation in sedimentary rocks[END_REF][START_REF] Scheibe | Hybrid multiscale simulation of a mixing-controlled reaction[END_REF][START_REF] Tartakovsky | Effects of incomplete mixing on multicomponent reactive transport[END_REF][START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF]. Anomalous transport can be observed at different scales: from pore [START_REF] Seymour | Anomalous Fluid Transport in Porous Media Induced by Biofilm Growth[END_REF][START_REF] Bijeljic | Signature of Non-Fickian Solute Transport in Complex Heterogeneous Porous Media[END_REF][START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-D porous media[END_REF], column [START_REF] Hatano | Dispersive transport of ions in column experiments: An explanation of long-tailed profiles[END_REF][START_REF] Rinaldo | Solute transport in low-heterogeneity sandboxes: The role of correlation length and permeability variance[END_REF] to field scale [START_REF] Garabedian | Large-Scale natural gradient tracer test in sand and gravel, Cape Cod, Massachusetts 2. Analysis of spatial moments for a nonreactive tracer[END_REF][START_REF] Becker | Tracer transport in fractured crystalline rock: Evidence of nondiffusive breakthrough tailing[END_REF][START_REF] Mckenna | Tracer tests in a fractured dolomite: 3. Double-porosity, multiple-rate mass transfer processes in convergent flow tracer tests[END_REF][START_REF] Kang | Impact of velocity correlation and distribution on transport in fractured media: Field evidence and theoretical model[END_REF][START_REF] Zech | Is unique scaling of aquifer macrodispersivity supported by field data?[END_REF]. An alternative formulation for anomalous transport that takes advection, dispersion, and mixing into consideration is, therefore, warranted [START_REF] De Dreuzy | Time evolution of mixing in heterogeneous porous media[END_REF][START_REF] De Dreuzy | On the validity of effective formulations for transport through heterogeneous porous media[END_REF].

The structure of heterogeneity is usually well known at large (km) scales from geological understanding. Thus, below a certain scale, it is necessary to address the absence of structure definition. At this point, the velocity development shows Markovianity in space rather than time [START_REF] Montroll | Random Walks on Lattices[END_REF][START_REF] Scher | Stochastic Transport in a Disordered Solid. I. Theory[END_REF][START_REF] Metzler | The random walk's guide to anomalous diffusion: A fractional dynamics approach[END_REF][START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Neuman | Perspective on theories of non-Fickian transport in heterogeneous media[END_REF][START_REF] Scher | Anomalous transit-time dispersion in amorphous solids[END_REF][START_REF] Klafter | Derivation of the Continuous-Time Random-Walk Equation[END_REF][START_REF] Le Borgne | Spatial Markov processes for modeling Lagrangian particle dynamics in heterogeneous porous media[END_REF]. A number of alternatives to the ADE have been proposed to address anomalous transport. The most widely extended is the continuous time random walk (CTRW). It consists of random velocity transitions once the solute has travelled a certain space step. However, these transitions are not fully random processes, but correlated ones [START_REF] Le Borgne | Lagrangian Statistical Model for Transport in Highly Heterogeneous Velocity Fields[END_REF]. Thus, a transition matrix M vs is needed [START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-D porous media[END_REF][START_REF] Kang | Spatial Markov Model of Anomalous Transport Through Random Lattice Networks[END_REF][START_REF] Benke | Modeling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF][START_REF] De Anna | Flow Intermittency, Dispersion, and Correlated Continuous Time Random Walks in Porous Media[END_REF][START_REF] Kang | Anomalous transport on regular fracture networks: Impact of conductivity heterogeneity and mixing at fracture intersections[END_REF]. M vs stands for the Energies 2021, 14, 6562 2 of 18 matrix probability to change the velocity state v given a fixed space phase s step. In this context, velocity becomes a variable of concentration such as space and time. The solute dependency of velocity was expressed in a phase space formulation proposed by [START_REF] Kang | Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes[END_REF]. However, this formulation does not take into account mixing.

Mixing is a consequence of diffusion among water bodies at a given time. Therefore, mixing is Markovian in time in contrast to dispersion. This observation suggests that solute transport should be localized not only in space and time, but also in velocity. The localization in time is of paramount importance, since mixing has a direct impact on chemical reactions in fluids [START_REF] Cirpka | Reply to comments on "Two-dimensional concentration distribution for mixing-controlled bioreactive transport in steady state[END_REF][START_REF] Rezaei | Reactive transport modeling of calcite dissolution in the fresh-salt water mixing zone[END_REF][START_REF] De Simoni | A procedure for the solution of multicomponent reactive transport problems[END_REF][START_REF] De Simoni | A mixing ratios-based formulation for multicomponent reactive transport[END_REF][START_REF] Tartakovsky | Mixing-induced precipitation: Experimental study and multiscale numerical analysis[END_REF]. The classic definition of mixing is proportional to the concentration gradient [START_REF] Fick | On liquid diffusion[END_REF]. In fast chemical reactions, the reaction rate may also be calculated with the concentration gradient [START_REF] De Simoni | A procedure for the solution of multicomponent reactive transport problems[END_REF]. However, [START_REF] Einstein | Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen[END_REF] demonstrated that mixing is due to the Brownian movement of single particles. Although this is a downscaling formulation of [START_REF] Fick | On liquid diffusion[END_REF] and leads to the same results, the comparison of these two expressions is inconsistent [START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF]. This inconsistency is evidenced by the fact that domains with constant concentration (where the concentration gradients are equal to zero) are defined as immobile zones by [START_REF] Fick | On liquid diffusion[END_REF], even though they are full of diffusing solute particles. As an alternative, the water mixing approach (WMA) has recently been proposed [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF][START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF]. WMA defines the mixing flux as an exchange of waters that carries the mean concentrations of the cells. Thus, mixing is dependent on the concentrations (i.e., number of particles that cross the cell interface) instead of on their gradients.

Dispersion and mixing are different processes, but they are closely linked [START_REF] Kitanidis | Prediction by the method of moments of transport in a heterogeneous formation[END_REF][START_REF] Kitanidis | The concept of the Dilution Index[END_REF]. Spreading is essentially driven by advection variability and tends to enhance the concentration contrast, which in turn enhances mixing [START_REF] Chiogna | Transverse mixing of conservative and reactive tracers in porous media: Quantification through the concepts of flux-related and critical dilution indices[END_REF][START_REF] Le Borgne | Non-Fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media[END_REF][START_REF] Rolle | Enhancement of dilution and transverse reactive mixing in porous media: Experiments and model-based interpretation[END_REF][START_REF] Tartakovsky | Stochastic Langevin Model for Flow and Transport in Porous Media[END_REF]. The link is evidenced by the stretching and folding processes [START_REF] De Anna | Mixing and Reaction Kinetics in Porous Media: An Experimental Pore Scale Quantification[END_REF][START_REF] Jiménez-Martínez | Pore-scale mechanisms for the enhancement of mixing in unsaturated porous media and implications for chemical reactions[END_REF][START_REF] Le Borgne | The lamellar description of mixing in porous media[END_REF]. This link leads to a non-Fickian mixing at earlier times over a considerable period [START_REF] Le Borgne | Non-Fickian dispersion in porous media: 2. Model validation from measurements at different scales[END_REF][START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Neuman | Perspective on theories of non-Fickian transport in heterogeneous media[END_REF][START_REF] Le Borgne | Lagrangian Statistical Model for Transport in Highly Heterogeneous Velocity Fields[END_REF][START_REF] Zhang | Time and space nonlocalities underlying fractional-derivative models: Distinction and literature review of field applications[END_REF]. Fickian mixing at later times is attributed to the spreading rate [START_REF] Le Borgne | Non-Fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media[END_REF].

Several formulations have been put forward to overcome the problems of the ADE [START_REF] Frippiat | A comparative review of upscaling methods for solute transport in heterogeneous porous media[END_REF]. At the continuum scale, alternative methods include CTRW [START_REF] Bijeljic | Signature of Non-Fickian Solute Transport in Complex Heterogeneous Porous Media[END_REF][START_REF] Le Borgne | Spatial Markov processes for modeling Lagrangian particle dynamics in heterogeneous porous media[END_REF][START_REF] Le Borgne | Lagrangian Statistical Model for Transport in Highly Heterogeneous Velocity Fields[END_REF][START_REF] Kang | Spatial Markov Model of Anomalous Transport Through Random Lattice Networks[END_REF][START_REF] Berkowitz | Anomalous Transport in Random Fracture Networks[END_REF][START_REF] Dentz | Time behavior of solute transport in heterogeneous media: Transition from anomalous to normal transport[END_REF][START_REF] Geiger | Upscaling solute transport in naturally fractured porous media with the continuous time random walk method[END_REF][START_REF] Selker | Capabilities and limitations of tracing spatial temperature patterns by fiber-optic distributed temperature sensing[END_REF][START_REF] Dentz | Continuous time random walks for non-local radial solute transport[END_REF][START_REF] Edery | Origins of anomalous transport in heterogeneous media: Structural and dynamic controls[END_REF][START_REF] Aquino | Chemical Continuous Time Random Walks[END_REF], fractional derivatives in ADE [START_REF] Cushman | Fractional advection-dispersion equation: A classical mass balance with convolution-Fickian Flux[END_REF][START_REF] Benson | Application of a fractional advection-dispersion equation[END_REF], stochastic convective stream tube [START_REF] Becker | Interpreting tracer breakthrough tailing from different forced-gradient tracer experiment configurations in fractured bedrock[END_REF], or multi-rate mass transfer [START_REF] Soler-Sagarra | Simulation of chemical reaction localization using a multi-porosity reactive transport approach[END_REF][START_REF] Babey | Multi-Rate Mass Transfer (MRMT) models for general diffusive porosity structures[END_REF][START_REF] Carrera | On matrix diffusion: Formulations, solution methods and qualitative effects[END_REF][START_REF] De Dreuzy | Influence of porosity structures on mixing-induced reactivity at chemical equilibrium in mobile/immobile Multi-Rate Mass Transfer (MRMT) and Multiple INteracting Continua (MINC) models[END_REF][START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF][START_REF] Fernandez-Garcia | Mathematical equivalence between time-dependent single-rate and multiratemass transfer models[END_REF]. Given that chemical reaction occurs at pore scale [START_REF] Steefel | Reactive transport modeling: An essential tool and a new research approach for the Earth sciences[END_REF], some pore-scale methods such as the lattice Boltzmann equation [START_REF] Benzi | The lattice Boltzmann equation: Theory and applications[END_REF][START_REF] Chen | Lattice Boltzmann Method for fluids flows[END_REF][START_REF] Kang | Lattice Boltzmann pore-scale model for multicomponent reactive transport in porous media[END_REF][START_REF] Acharya | Pore-scale simulation of dispersion and reaction along a transverse mixing zone in two-dimensional porous media[END_REF][START_REF] Willingham | Evaluation of the Effects of Porous Media Structure on Mixing-Controlled Reactions Using Pore-Scale Modeling and Micromodel Experiments[END_REF], smoothed particle hydrodynamics [START_REF] Tartakovsky | Effects of incomplete mixing on multicomponent reactive transport[END_REF][START_REF] Tartakovsky | Simulations of reactive transport and precipitation with smoothed particle hydrodynamics[END_REF][START_REF] Tartakovsky | Smoothed particle hydrodynamics and its applications for multiphase flow and reactive transport in porous media[END_REF][START_REF] Liu | Smoothed Particle Hydrodynamics (SPH): An Overview and Recent Developments[END_REF], pore network models [START_REF] Meile | Scale dependence of reaction rates in porous media[END_REF][START_REF] Li | Upscaling geochemical reaction rates using pore-scale network modeling[END_REF][START_REF] Blunt | Flow in porous media-Pore-network models and multiphase flow[END_REF][START_REF] Blunt | Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow[END_REF][START_REF] Raoof | Upscaling transport of adsorbing solutes in porous media: Pore-network modeling[END_REF][START_REF] Raoof | A new formulation for pore-network modeling of two-phase flow[END_REF][START_REF] Varloteaux | Modélisation Multi-Echelles des Mécanismes de Transport Réactif[END_REF], or density kernel methods [START_REF] Sole-Mari | Particle density estimation with grid-projected and boundarycorrected adaptive kernels[END_REF][START_REF] Schmidt | A Kernel-based Lagrangian method for imperfectly-mixed chemical reactions[END_REF] have been studied. Hybrid continuum-pore-scale methods have also been proposed [START_REF] Battiato | Hybrid models of reactive transport in porous and fractured media[END_REF][START_REF] Tartakovsky | Hybrid simulations of reaction-diffusion systems in porous media[END_REF][START_REF] Van Leemput | Accuracy of Hybrid Lattice Boltzmann/Finite Difference Schemes for Reaction-Diffusion Systems[END_REF]. All these approaches partially meet the requirements listed by [START_REF] De Dreuzy | On the validity of effective formulations for transport through heterogeneous porous media[END_REF].

Recently, the multi-advective water mixing approach (MAWMA) was advanced by [START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF][START_REF] Soler-Sagarra | Modeling Mixing in Stratified Heterogeneous Media: The Role of Water Velocity Discretization in Phase Space Formulation[END_REF] to simulate anomalous transport. MAWMA can be taken as a WMA extension. The MAWMA is a phase space formulation (detailed in Section 2) where velocity is a new dimension of the state variable. The anomalous nature of transport is addressed by acknowledging that concentrations will be variable, within a representative elementary volume, and that this variability correlates with velocity (e.g., concentration in an advancing front will tend to be larger in high velocity portions of the medium than in nearly stagnant portions). The water parcel method (WP) can be derived by discretizing space, time, and velocity (Section 3.1). The question is how to exchange solute mass among water parcels. Exchange will depend on the properties of the medium (notably hydraulic conductivity and porosity) and, especially, on heterogeneity. In this paper, we test the capacity of the WP model and the proposed mass exchange methods to reproduce transport through heterogeneous porous media.

Governing Equations

Soler-Sagarra et al. [START_REF] Soler-Sagarra | Modeling Mixing in Stratified Heterogeneous Media: The Role of Water Velocity Discretization in Phase Space Formulation[END_REF] proposed a phase space formulation (MAWMA) to describe solute transport through porous media while meeting the requirements of [START_REF] De Dreuzy | On the validity of effective formulations for transport through heterogeneous porous media[END_REF]. As a phase space formulation, MAWMA considers that concentration depends not only on space and time, but also on velocity, i.e., c = c(x,v,t) [M/L 3 ]. Using velocity as a new dimension facilitates describing dispersion, because spreading results naturally from Energies 2021, 14, 6562 3 of 18 velocity variability, so that no explicit accounting is needed for dispersion. Advection and mixing processes are defined by fluxes f in both the space, s, and velocity, v, domains. That is,

φ ∂c(x, v, t) ∂t = f adv,s + f di f f ,s + f adv,v + f di f f ,v + r (1)
where φ [L 3 /L 3 ] is porosity, t [T] is time, and r [M/L 3 /T] is a sink/source term, possibly reflecting chemical reactions. The first term on the right hand side (rhs) represents advection within a velocity class, traditionally expressed in terms of Darcy flux, which we prefer to write here as a function of velocity as q = φv [L 3 /L 2 /T] (blue arrow in Figure 1) as follows

f adv,s = -φv • ∇c (2)
reflecting chemical reactions. The first term on the right hand side (rhs) represents adv tion within a velocity class, traditionally expressed in terms of Darcy flux, which we pr to write here as a function of velocity as q = 𝜙v [L 3 /L 2 /T] (blue arrow in Figure 1) as follo 𝑓 𝑎𝑑𝑣,𝑠 = -𝜙𝑣 ⋅ ∇𝑐

The second term of the rhs represents diffusion within a velocity class, classic defined by Fick's law (orange arrow in Figure 1). However, the water mixing appro (WMA) was proposed as an alternative by [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF][START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF])

𝑓 𝑑𝑖𝑓𝑓,𝑠 = -𝛻 ⋅ (𝑞 𝐷 𝑐 ̅̅̅̅̅)
where the term 𝑞 𝐷 𝑐 ̅̅̅̅̅ is used to express solute exchanges associated to water mass changes, 𝑞 𝐷 = 𝜙𝐷 𝑤 /𝐿 𝐷 [L 3 /L 2 /T] is the water diffusion flux exchanged, 𝐷 𝑤 [L 2 /T] is the ter molecular diffusion coefficient, and 𝐿 𝐷 [L] is a characteristic diffusion scale. Equa (3) expresses diffusion as the exchange of water depending on the concentration inst of on its gradient. The third term on the rhs of Equation (1) represents changes in c(x,v,t) due to chan in velocity (green arrow in Figure 1). [START_REF] Kang | Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes[END_REF]) expressed it in terms of solute mass probab p. We express it in terms of concentration considering that 𝑝 = 𝜙𝑐 𝑀 ⁄ (𝑀 being the t solute mass), which yields The second term of the rhs represents diffusion within a velocity class, classically defined by Fick's law (orange arrow in Figure 1). However, the water mixing approach (WMA) was proposed as an alternative by [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF][START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF])

f di f f ,s = -∇ • (q D c) (3)
where the term q D c is used to express solute exchanges associated to water mass exchanges, q D = φD w /L D [L 3 /L 2 /T] is the water diffusion flux exchanged, D w [L 2 /T] is the water molecular diffusion coefficient, and L D [L] is a characteristic diffusion scale. Equation (3) expresses diffusion as the exchange of water depending on the concentration instead of on its gradient. The third term on the rhs of Equation (1) represents changes in c(x,v,t) due to changes in velocity (green arrow in Figure 1). [START_REF] Kang | Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes[END_REF]) expressed it in terms of solute mass probability p. We express it in terms of concentration considering that p = φc/M (M being the total solute mass), which yields

f adv,v = - v L A φc + v g vs v v v L A φc dv + r (4)
where L A [L] is the characteristic advection length, g vs [TL -1 ] is probability density of a velocity transition after covering a step L A in space, and c = c(x,v ,t). Note that the first term on the rhs of Equation ( 4) refers to transitions to velocity v. It does not involve any velocity integration because v g vs (v |v)dv = 1.

Finally, [START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF][START_REF] Soler-Sagarra | Modeling Mixing in Stratified Heterogeneous Media: The Role of Water Velocity Discretization in Phase Space Formulation[END_REF] proposed expressing diffusive transitions in velocity (purple arrow in Figure 1) as

f di f f ,v = v φ f vt v v c dv - v φ f vt v v cdv (5)
where f vt [L -1 ] is the probability density, per unit time, of diffusive transitions between velocity states v. The expression φ f vt has units of water flux. Note that, as in WMA, all fluxes are expressed in terms of water instead of solute concentrations, which become a mere attribute of water. This is why we termed the formulation of Equation ( 4) multi-advective water mixing approach (MAWMA).

Equation ( 4) could also be expressed as a Lagrangian formulation. This requires revising the definition of a material derivative. D(•)/dt reflects all changes in a flowing element of water. Therefore, it is expressed as the partial derivative minus the changes in c caused by advection. Since we are defining f adv,v to represent advective velocity transitions, we can write the material derivative as

φ dc dt = φ ∂c ∂t + φv • ∇c + v L A φc - v g vs v v v L A φc dv (6)
This definition acknowledges that velocity transitions due to heterogeneity do not cause mixing, which helps us to focus on the impact of mixing, which depends exclusively on diffusive processes:

φ dc dt = -∇ • (q D c) + v φ f vt v v c dv - v φ f vt v v cdv (7)

Solution Method

The equation presented in the previous section can be solved with any numerical methods. Here, we present a modeling option (Section 3.1), termed the water parcel (WP) method, which is an extension of the one proposed by [START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF][START_REF] Soler-Sagarra | Modeling Mixing in Stratified Heterogeneous Media: The Role of Water Velocity Discretization in Phase Space Formulation[END_REF]. We present first this new extension, which require two transition probability matrices. We then describe the computation of these matrices and their properties.

Water Parcel Method

Time, space, and velocity must be discretized to solve Equation [START_REF] Alcolea | Regularized pilot points method for reproducing the effect of small scale variability: Application to simulations of contaminant transport[END_REF]. For simplicity, the discretization procedure is the one used by [START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF][START_REF] Soler-Sagarra | Modeling Mixing in Stratified Heterogeneous Media: The Role of Water Velocity Discretization in Phase Space Formulation[END_REF]. We reduce the spatial dimensions to 1 by integrating the dimensions perpendicular to the mean flow. The velocity field of the entire domain is discretized in classes with the same flux (i.e., equally probable in the Lagrangian velocity distribution). Since fast velocities concentrate most of the flux, they are less probable than slow velocities from an Eulerian point of view (they occupy less volume than slow velocities). Details of this velocity discretization are reported in [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]).

The space and velocity domains are discretized in parcels with the same water volume. Each water parcel is associated to its centroid, which determines the position (in space and velocity) at a given time. The length of a single parcel (i.e., its space extension along the x coordinate) is proportional to its velocity as reported by [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF][START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF]. As a result, slow velocities have more parcels than fast velocities (see Figure 2a), which is consistent with their higher Eulerian probability (i.e., p e in [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]). The width is proportional to its probability (inversely proportional to its velocity). Equation ( 7) is integrated into the continuum space-velocity by using shape functions associated with each parcel. As in the finite volume method, the shape function equals 1 when (x,v) exists in the parcel domain. Otherwise, it equals to 0. Therefore, all attributes of water parcels (e.g., concentration) will be regarded as homogeneous within each parcel.

Energies 2021, 14, x FOR PEER REVIEW 5 of 19 integrated into the continuum space-velocity by using shape functions associated with each parcel. As in the finite volume method, the shape function equals 1 when (x,v) exists in the parcel domain. Otherwise, it equals to 0. Therefore, all attributes of water parcels (e.g., concentration) will be regarded as homogeneous within each parcel. Parcels are injected and advected through the domain similar to solute particles. The injection velocity class is assigned randomly with equal probability for all classes to represent a uniform flux averaged injection.

Advection (Equations ( 2) and ( 4)) is simulated by simply displacing the parcel centroid with its associated velocity until it has covered the distance 𝐿 𝐴 . Then, a random event is performed to assign a new velocity for the next space step 𝐿 𝐴 according to transition probabilities given by the transition matrix 𝑴 𝒂𝒅𝒗 𝒗𝒔 (𝑀 𝑖𝑗 is the probability of jumping from velocity 𝑣 𝑗 to 𝑣 𝑖 ). This transition matrix is similar to the classic solute transition matrix 𝑴 𝒗𝒔 of [START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-D porous media[END_REF][START_REF] Kang | Impact of velocity correlation and distribution on transport in fractured media: Field evidence and theoretical model[END_REF][START_REF] Le Borgne | Spatial Markov processes for modeling Lagrangian particle dynamics in heterogeneous porous media[END_REF][START_REF] Kang | Spatial Markov Model of Anomalous Transport Through Random Lattice Networks[END_REF][START_REF] De Anna | Flow Intermittency, Dispersion, and Correlated Continuous Time Random Walks in Porous Media[END_REF][START_REF] Kang | Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes[END_REF] except that it only accounts for advection transitions. Since the simulation takes place with fixed time steps, each parcel will take a different number of steps to perform the next random event. The advection scheme is plotted in Figure 2b. The discretized form of Equation [START_REF] Alcolea | Regularized pilot points method for reproducing the effect of small scale variability: Application to simulations of contaminant transport[END_REF]

where Vw [L 3 ] is the water volume of each parcel, ∆𝑡 [T] is the time step, Nli is the number of parcels 𝑦 with velocity l spatially connected to i, F is the water volumetric flux diffused, Parcels are injected and advected through the domain similar to solute particles. The injection velocity class is assigned randomly with equal probability for all classes to represent a uniform flux averaged injection.

Advection (Equations ( 2) and ( 4)) is simulated by simply displacing the parcel centroid with its associated velocity until it has covered the distance L A . Then, a random event is performed to assign a new velocity for the next space step L A according to transition probabilities given by the transition matrix M vs adv (M ij is the probability of jumping from velocity v j to v i ). This transition matrix is similar to the classic solute transition matrix M vs of [START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-D porous media[END_REF][START_REF] Kang | Impact of velocity correlation and distribution on transport in fractured media: Field evidence and theoretical model[END_REF][START_REF] Le Borgne | Spatial Markov processes for modeling Lagrangian particle dynamics in heterogeneous porous media[END_REF][START_REF] Kang | Spatial Markov Model of Anomalous Transport Through Random Lattice Networks[END_REF][START_REF] De Anna | Flow Intermittency, Dispersion, and Correlated Continuous Time Random Walks in Porous Media[END_REF][START_REF] Kang | Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes[END_REF] except that it only accounts for advection transitions. Since the simulation takes place with fixed time steps, each parcel will take a different number of steps to perform the next random event. The advection scheme is plotted in Figure 2b.

The discretized form of Equation ( 7) in a single parcel I of velocity class l is given by

V w c k+1 i -c k i ∆t = N li ∑ y =i F iy c y -c i + N v ∑ m =l N mi ∑ j a ij F vt lm c k j - N v ∑ m =l N mi ∑ j a ji F vt ml c k i i ∈ I l (8)
where V w [L 3 ] is the water volume of each parcel, ∆t [T] is the time step, N li is the number of parcels y with velocity l spatially connected to i, F is the water volumetric flux diffused, N v is the number of velocity classes, N mi is the number of parcels of velocity class m connected with parcel i, and a ij is the fraction of the diffusive flux assigned to velocity m, which will be exchanged with the j parcel. Recall that we interpret mixing as a water exchange process derived from water diffusion, which implies symmetry (i.e., F iy = F yi , F vt lm = F vt ml and a ij = a ji ) (the latter requires post processing). We can rewrite Equation ( 8) using the concept of mixing ratio λ = a ∆t F/V w [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF], which leads to

c k+1 i = λ ii c k i + N li ∑ y =i λ iy c k y + N v ∑ m =l N mi ∑ j λ ij c k j ( 9
)
The sum of all λ equals 1, because the coefficients multiplying c k j (∀j = i) on the rhs of Equation ( 8) are always compensated by the same coefficient multiplying -c k i . Therefore,

λ ii = 1 -∑ N li y λ iy -∑ N v m ∑ N mi j λ ij .
A necessary and sufficient condition for stability is that λ ii > 0 ∀i.

Equation ( 9) represents a mixing equation with mixing ratios that are independent of the species, which is consistent with the fact that it has been derived from the mixing of waters. Note that all transport processes described above involve water transfers. [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF][START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF] explain the extension from Equation ( 9) to reactive transport.

Random Walk Method

The WP method described in Section 3.1 requires the velocity distribution, the diffusive transition matrix M vt mix for mixing, and the advective transition matrix M vs adv . Here, we compute these from random walk simulations of flow and transport, which we will also use to test the proposed approach. The model to simulate flow is essentially that of [START_REF] Hakoun | Upscaling and Prediction of Lagrangian Velocity Dynamics in Heterogeneous Porous Media[END_REF]. We summarize it for the sake of completeness. A 2D multi-lognormal random permeability field K(x) is generated with an isotropic Gaussian covariance function

Y (x)Y x = σ 2 Y exp - |x| L A 2 ( 10
)
where

Y (x) = Y(x) -Y(x), Y(x) = ln(K(x)),
v(x) = - K(x)∇h(x) φ ( 11
)
where h is the hydraulic head. Fixed head boundary conditions are imposed to the upstream and downstream boundaries. No-flow conditions are imposed to the top and bottom boundaries. The flow model is solved by using the finite volume method to first obtain heads and, then, using Equation [START_REF] Battiato | On breakdown of macroscopic models of mixing-controlled heterogeneous reactions in porous media[END_REF] for the velocity field v(x).

A Python code is employed to solve particle transport. Particle advection is calculated as in [102]. Diffusion displacement at a given time step is L D ξ, where ξ∼ N(0, 1).

Algebra of Mixing Matrices

The water mixing ratio λ ij defined in Section 3.1 can be understood as the ij-th position of the water transition matrix M vt mix applied to time step. Here, we describe how to compute the mixing matrix M vt mix from an RW simulation above the procedure as general, in the sense that it can be employed for advection transition matrices (applied after space steps) or mixing transition matrices (e.g., [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF].

Transition matrices are obtained directly from their Markovian probability definition (i.e., M ij is the probability of a particle to end in velocity class i, given that it started in class j, which implicitly carries the Markovian statement that the next state solely depends on the current state). Therefore, M ij = N k+1 ij /N k j , where N k+1 ij is the number of particles that ended in velocity class i at time t k+1 after a diffusion step (to avoid advection transitions) having started in class j at time t k and N k j . M vs adv is computed analogously, except that accounting is made not at every time step, but after the particle has covered the characteristic advection scale.

Two issues need to be addressed. First, the above definition refers to probabilities, while we need volumetric water exchanges. Second, Markovianity needs to be tested. It was demonstrated by [START_REF] Le Borgne | Spatial Markov processes for modeling Lagrangian particle dynamics in heterogeneous porous media[END_REF] for advection transitions, and it would be trivially true for mixing transitions in the absence of advection. However, it is not so clear when coupling advection and diffusion, especially when considering that low velocities occupy a much larger volume than high velocities. We will test Markovianity as part of the example in Section 4. However, we need first to clarify the relationship between transition probabilities and mixing matrices.

The relationship between the vector of solute probabilities p (p t according to [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]) and velocity class concentration is expressed as

p = m T -1 Sc ( 12
)
where m T is the total mass, and S is the storage matrix containing the volume of each class. S is not expressed in Eularian processes (time steps) as in Lagrangian processes (space steps). The procedure explained above to obtain the probability transition matrix M p counting particle transitions in RW simulations can be expressed as

p k+1 = M p p k (13)
Combining ( 12) and ( 13)

Sc k+1 = M p Sc k (14)
Then, the concentrations for the next time step are as follows

c k+1 = S -1 M p Sc k (15)
Therefore, the transition matrix for transport simulations M c can be obtained from the RW matrix

M p M c = S -1 M p S (16)
A well-known property of Markov probability transition matrices is that the sum of the columns of M p equals 1 (a particle in any velocity class must end in some class). However, the rows of the M c must add up to 1 to express that concentrations do not change if equal in all velocity classes. In fact, component M Cij can be viewed that as the volume of water received by class i from class j, expressed as a fraction of the volume in i (i.e., a mixing ratio), so that the rows must add up to 1 to also ensure that the class volume does not change. Therefore, the volume of water exchanged is expressed as

M V = M p S (17)
To satisfy mass conservation, water volume exchanged between velocity classes i and j must be equal (i.e., M V must be symmetric). The computation procedure (starting with the probability transition matrix) does not ensure symmetry. In practice, it is nearly so. So, symmetry is imposed by setting

M V = M V t + M V 2 (18)
Energies 2021, 14, 6562

8 of 18
where M V is the volume exchange matrix computed initially from Equation [START_REF] Bijeljic | Signature of Non-Fickian Solute Transport in Complex Heterogeneous Porous Media[END_REF]. Finally, the water transition matrix is expressed as

M w = S -1 M V (19)
While the properties of transition matrices are clear, and the equations in this section facilitate generalizing them for water parcels with varying volumes, which may be of interest for future applications, the question remains about the robustness of the proposed calculation method. Specifically, is it necessary to do a full RW simulation to compute the transition matrix?

Applications

Although WP is a continuum scale method, its solute evolution must reproduce the particle-based behavior [START_REF] Einstein | Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen[END_REF]. This is why it is tested with the classic RW presented in Section 3.2. The model parameters are detailed in Table 1. Three different Peclet number simulations are defined: ∞, 1000, and 50. The Peclet number is defined as follows Initial concentrations are defined in both methods. The WP method employs the initial flux weighted distribution of solute mass (Figure 2a)

Pe = L A v /D w (20)
c i (t = 0) = c re f v i v (21)
where c re f is the initial concentration reference, the angular bracket • denotes the mean injection velocity (mean of the Lagrangian distribution), and v i is the parcel velocity.

In order to simulate a water parcel distribution, each particle of the RW method has an initial time step with a random definition ∆t 0 = ∆t•γ, being γ∼unif (0,1). This definition provides an innovative way to simulate transport, since it differs from the classic Dirac delta. We believe it is a realistic situation, as it reproduces water injection. Other improvements from the classic RW method have been proposed [103]. Here, an initial number of particles N p are placed along the domain width L y at the fixed x 0 coordinate position. The particles have an initial flux-weighted distribution. This means that each cell in x 0 has N c particles, which is a function of the cell velocity v c expressed as follows

N c = N p L y /∆y v c v ( 22
)
where ∆y is the cell width. In order to simulate an injected concentration equal to 1, the mass of a single particle m p is

m p = c re f φ v ∆tL y N p ( 23
)
Energies 2021, 14, 6562

9 of 18
The WP method should reproduce mean advection, dispersion, mixing, and "be flexible enough to be applicable to real problems" [START_REF] De Dreuzy | On the validity of effective formulations for transport through heterogeneous porous media[END_REF]. This latter condition is somewhat subjective and will not be considered here, but we believe that WP may be applied to field cases because (i) it is defined at the continuum scale, so that it can benefit from traditional hydrogeological characterization methods, (ii) it localizes concentration in the (x,v,t) continuum domain, and (iii) it is easily extended to reactive transport [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF][START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF]. Still, a number of developments are needed to address the real cases with a level of maturity comparable to stochastic methods [START_REF] Neuman | Perspective on theories of non-Fickian transport in heterogeneous media[END_REF]104]. Therefore, we restrict ourselves to test advection, dispersion, and mixing on the synthetic case for stationary conditions and mean uniform flow. The spatial distribution of solute concentration for times 5, 20, and 100 s are shown in Figure 3.

gies 2021, 14, x FOR PEER REVIEW 9 of subjective and will not be considered here, but we believe that WP may be applied to fi cases because (i) it is defined at the continuum scale, so that it can benefit from traditio hydrogeological characterization methods, (ii) it localizes concentration in the (x,v,t) co tinuum domain, and (iii) it is easily extended to reactive transport [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF][START_REF] Soler-Sagarra | Mathematical Formulations of Water Mixing for Reactive Transport through Heterogeneous Media[END_REF]. Still, a num of developments are needed to address the real cases with a level of maturity compara to stochastic methods [START_REF] Neuman | Perspective on theories of non-Fickian transport in heterogeneous media[END_REF]104]. Therefore, we restrict ourselves to test advection, disp sion, and mixing on the synthetic case for stationary conditions and mean uniform flo The spatial distribution of solute concentration for times 5, 20, and 100 s are shown Figure 3. The mean advection is characterized by the mean position µ defined as

µ(t) = m (1) x (t) -x 0 ( 24
)
where m

(k)

x is the k-th order moment of the solute distribution in space

m (k) x (t) = Ω cx k dΩ/ Ω cdΩ ( 25
)
where Ω is the flow domain. From the above definition, we can express dispersion by the standard deviation of spatial solute distribution σ 2 x , which is described as

σ 2 x (t) = m (2) x (t) -µ(t) 2 (26)
Global mixing G [105] is defined as

G(t) = Ω c 2 dΩ ( 27
)
Energies 2021, 14, 6562 10 of 18

Note that we can also define global mixing G in terms of the velocity domain such as

G (t) = v c 2 dv (28)
where c = c(v,t) is the mean concentration of an velocity class.

Transition Matrix Validation with Markovian Models

We defined three transition matrices in Section 3.3: M p , M c , and M w . We test here the validity of their computation using a Markov chain model [106]. We first compute the probability transition matrix M p from RW (at σ 2 Y = 1 and Pe = 50) simulations at three different times: t = 1, 5, and 250. The last time corresponds to the characteristic diffusive time (λ 2 /2/D w ), so that we can assume that injected particles have sampled exhaustively the whole velocity space (recall that we are using flux averaged injection, so that the slow velocities volume is less exhaustively sampled than the fast velocities volume).

Matrices M c and M w are calculated as explained in Section 3.3. Equation ( 13) defines the step computation for the matrices M p . The Markovian models that employ M c and M w use concentration c instead of p as the state variable. Concentration c is readily converted to p by using [START_REF] Sadhukhan | A simulation study of reactive flow in 2-D involving dissolution and precipitation in sedimentary rocks[END_REF]. The initial solute probability distribution for any velocity class i is

p 0 i = 1/N v .
The computed evolution of G (Equation ( 28)) in time is shown in Figure 4. The first observation from Figure 4 is that σ decreases in time, which reflects that a uniform flow averaged probability leads to a non-uniform initial concentration. That is, the same mass flux occurs in all velocity classes, but concentration is much longer in the high velocity classes. Mixing causes the concentration to become uniform in all classes, so that G tends to decrease. Second, the M p and M c models evolve identically in time, although they are sensitive to the RW step in which they were computed. Third, we observe that the M w models always give identical results even though they are obtained from different M p matrices. The evolution computed with any of the M w matrices is identical to the one computed with the M p and M c matrices obtained after 250 RW steps. The evaluation method proposed for the M w matrix is indeed robust.

𝐺′(𝑡) = ∫ 𝑐′ 2 𝑑𝑣 𝑣 [START_REF] Montroll | Random Walks on Lattices[END_REF] where c' = c(v,t) is the mean concentration of an velocity class.

Transition Matrix Validation with Markovian Models

We defined three transition matrices in Section 3.3: 𝑴 𝒑 , 𝑴 𝒄 , and 𝑴 𝒘 . We test here the validity of their computation using a Markov chain model [106]. We first compute the probability transition matrix 𝑴 𝒑 from RW (at 𝜎 𝑌 2 = 1 and Pe = 50) simulations at three different times: t = 1, 5, and 250. The last time corresponds to the characteristic diffusive time (𝜆 2 2/𝐷 𝑤 ⁄ ), so that we can assume that injected particles have sampled exhaustively the whole velocity space (recall that we are using flux averaged injection, so that the slow velocities volume is less exhaustively sampled than the fast velocities volume).

Matrices 𝑴 𝒄 and 𝑴 𝒘 are calculated as explained in Section 3.3. Equation ( 13) defines the step computation for the matrices 𝑴 𝒑 . The Markovian models that employ M and Mw use concentration c instead of p as the state variable. Concentration c is readily converted to p by using [START_REF] Sadhukhan | A simulation study of reactive flow in 2-D involving dissolution and precipitation in sedimentary rocks[END_REF]. The initial solute probability distribution for any velocity class i is 𝑝 𝑖 0 = 1/𝑁 𝑣 . The computed evolution of G' (Equation ( 28)) in time is shown in Figure 4. The first observation from Figure 4 is that σ decreases in time, which reflects that a uniform flow averaged probability leads to a non-uniform initial concentration. That is, the same mass flux occurs in all velocity classes, but concentration is much longer in the high velocity classes. Mixing causes the concentration to become uniform in all classes, so that G' tends to decrease. Second, the Mp and Mc models evolve identically in time, although they are sensitive to the RW step in which they were computed. Third, we observe that the 𝑴 𝒘 models always give identical results even though they are obtained from different 𝑴 𝒑 matrices. The evolution computed with any of the 𝑴 𝒘 matrices is identical to the one computed with the Mp and Mc matrices obtained after 250 RW steps. The evaluation method proposed for the 𝑴 𝒘 matrix is indeed robust. The mixing state decreases (higher G' in Equation ( 28)) in the 𝑴 𝒑 models from time t = 1 until reaching t = 5, when the poorest mixing state is attained. A state identical to 𝑴 𝒘 The mixing state decreases (higher G in Equation ( 28)) in the M p models from time t = 1 until reaching t = 5, when the poorest mixing state is attained. A state identical to M w is reached at the characteristic time of diffusion (t = 250), confirming that the water transitions are always constant. This occurs despite the heterogeneity of solute distributions within the velocity class and is of major significance because mixing can be defined in a constant water transition matrix during the entire simulation, which is not the case with the solute matrices.

Comparison between RW and WP in Transport through Heterogeneous Porous Media

Mean position, spreading, and mixing results for MAWMA and RW and σ 2 Y = 1 are shown in Figures 5 and6. A good fit of mean position µ (Equation ( 24)) can be observed for all cases in Figure 5a. Regarding spreading, the evolution of σ 2

x (Equation ( 26)), using the RW, is consistent with those of [107,108] also showed that early time spreading of a Dirac Delta is controlled by velocity heterogeneities, thus becoming ballistic, so that the spatial variance is proportional to t 2 . This explains differences in the WP and RW results in Figure 5b. The WP grid is too coarse to reproduce early time advection of a Dirac pulse. This is not of great concern for practical applications, because this pulse initial condition is unusual in practice. The RW and WP results converge at later times, during the intermediate regime observed by [108] prior to the Fickian regime.

Energies 2021, 14, x FOR PEER REVIEW 11 of 19 is reached at the characteristic time of diffusion (t = 250), confirming that the water transitions are always constant. This occurs despite the heterogeneity of solute distributions within the velocity class and is of major significance because mixing can be defined in a constant water transition matrix during the entire simulation, which is not the case with the solute matrices.

Comparison between RW and WP in Transport through Heterogeneous Porous Media

Mean position, spreading, and mixing results for MAWMA and RW and 𝜎 𝑌 2 = 1 are shown in Figure 5 and Figure 6. A good fit of mean position 𝜇 (Equation ( 24)) can be observed for all cases in Figure 5a. Regarding spreading, the evolution of 𝜎 𝑥 2 (Equation ( 26)), using the RW, is consistent with those of [107,108].also showed that early time spreading of a Dirac Delta is controlled by velocity heterogeneities, thus becoming ballistic, so that the spatial variance is proportional to t 2 . This explains differences in the 𝑊𝑃 and 𝑅𝑊 results in Figure 5b. The 𝑊𝑃 grid is too coarse to reproduce early time advection of a Dirac pulse. This is not of great concern for practical applications, because this pulse initial condition is unusual in practice. The 𝑅𝑊 and 𝑊𝑃 results converge at later times, during the intermediate regime observed by [108] prior to the Fickian regime. We compare the results of global mix G plotted in Figure 6. We first evaluate t 𝑃𝑒 = ∞ case. The mismatch observed is due to the mesh evaluation. WP displays the c rect constant G because the entire solute remains in the initial parcel. In other words, mixing transition of solute occurs between parcels. By contrast, RW uses a structur mesh that is fixed for the evaluation of concentration. The number of concentrated e ments increases with time due to stretching [109], which implies a reduction in the co puted concentration; this is a common problem when comparing Eulerian (RW) and L grangian (WP) evaluations of concentration. For a more accurate comparison, RW shou therefore, be performed with a Lagrangian mesh (such as the one proposed by [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF]).

In the other Pe cases, RW shows a monotonic decrease in G. However, WP under timates mixing at early times. This mismatch is attributed to the mesh distinction giv that similar discrepancies are again observed. In contrast, WP overestimates mixing w respect to RW at late times. This is consistent with the WP standard deviation behav observed previously. We suspect this overestimation is related to incomplete mixi [110,111]. In other words, concentration variability remains at scales lower than veloc discretization, whereas WP changes the entire parcel concentration in the mixing proce (i.e., homogenization), which confirms the mixing results observed in the above sectio We compare the results of global mix G plotted in Figure 6. We first evaluate the Pe = ∞ case. The mismatch observed is due to the mesh evaluation. WP displays the correct constant G because the entire solute remains in the initial parcel. In other words, no mixing transition of solute occurs between parcels. By contrast, RW uses a structured mesh that is fixed for the evaluation of concentration. The number of concentrated elements increases with time due to stretching [109], which implies a reduction in the computed concentration; this is a common problem when comparing Eulerian (RW) and Lagrangian (WP) evaluations of concentration. For a more accurate comparison, RW should, therefore, be performed with a Lagrangian mesh (such as the one proposed by [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF]).

In the other Pe cases, RW shows a monotonic decrease in G. However, WP underestimates mixing at early times. This mismatch is attributed to the mesh distinction given that similar discrepancies are again observed. In contrast, WP overestimates mixing with respect to RW at late times. This is consistent with the WP standard deviation behavior observed previously. We suspect this overestimation is related to incomplete mixing [110,111]. In other words, concentration variability remains at scales lower than velocity discretization, whereas WP changes the entire parcel concentration in the mixing process (i.e., homogenization), which confirms the mixing results observed in the above section. Some alternatives Energies 2021, 14, 6562 13 of 18 to incomplete mixing have been proposed [107]. However, the mixing criterion of these authors is dependent on the solute state, whereas we seek to base it on the water volumes.

Conclusions

We present and test the MAWMA formulation for transport through heterogeneous porous media using the WP solution method. The formulation is an extension of WMA, which considers transport of water instead of solute concentration. Exchange of water volumes is used to reproduce mixing instead of individual species diffusion. Individual species concentrations are considered attributes of water. They may vary spatially, in which case the net solute mass exchange turns out to be proportional to the concentration gradient. However, water and solutes exchanges occur independently of concentration gradients, which is why no concentration gradient is used to calculate the mixing process.

We use the WP method to reproduce MAWMA. WP requires a velocity discretization and two transition matrices: one to reproduce advection transitions, which is Markovian in space (i.e., transitions occur after fixed spatial steps, which is consistent with a fixed heterogeneity), and one to reproduce mixing, which is Markovian in time (i.e., velocity transitions occur at a constant rate in time, which is consistent with Brownian motion). We have described how to compute these matrices from RW models. Our study shows that it is possible to obtain the water transition matrices from the classic solute transition matrices.

We use Markovian models with transition matrices computed from different time steps of the RW simulations to compare transition matrices. We showed that, unlike M c , M w is invariant in time, and the proposed calculation method is robust. Then, we test the performance of WP of advection, dispersion, and mixing with the RW simulations. We find that advection (mean displacement) and dispersion (mean spreading) are similar for the two methods. However, mixing (2nd moment of concentration) results are different, which we attribute to two reasons: (a) mesh inequivalence and (b) incomplete mixing. The structured mesh of RW is Eulerian in contrast to the unstructured Lagrangian mesh used in WP. The WP method assumes a homogeneous concentration value within the parcel volume. However, in our study, this is inappropriate because the mixing scale (i.e., the scale at which concentration can be considered constant) is much smaller than the scale required for velocity definition.

In summary, the RW concentration evaluation requires a Lagrangian mesh (such as the isochronal one proposed by [START_REF] Soler-Sagarra | Water Mixing Approach (WMA) for Reactive Transport Modeling[END_REF]). Moreover, WP needs a new evaluation of concentration in order to take into account heterogeneity inside the parcels. This new evaluation should consider water volumes, which will ensure the independence of concentration states. This will facilitate coupling with chemical reaction calculations.

Figure 1 .

 1 Figure 1. Scheme of particle transport processes through a continuum heterogeneous domain. The left image is a computed velocity field. The right-top image displays the advection path of two particles. The bottom-right image shows the diffusion possibilities of a single particle.

Figure 1 .

 1 Figure 1. Scheme of particle transport processes through a continuum heterogeneous domain. The left image is a computed velocity field. The right-top image displays the advection path of two particles. The bottom-right image shows the diffusion possibilities of a single particle.

Figure 2 .

 2 Figure 2. Scheme of water parcels (WP) method using the multi-advective water mixing approach formulation with three velocity classes: (a) initial distribution of parcels in the (x,v) domain and the initial concentration condition, (b) advective process for a single water parcel, and (c) mixing process for a single water parcel.

Figure 2 .

 2 Figure 2. Scheme of water parcels (WP) method using the multi-advective water mixing approach formulation with three velocity classes: (a) initial distribution of parcels in the (x,v) domain and the initial concentration condition, (b) advective process for a single water parcel, and (c) mixing process for a single water parcel.

and σ 2 Y

 2 is the log-permeability variance. The random field Y(x) with mean µ Y = Y(x) is generated using the Random Fields Package [100] of the R software environment [101]. Groundwater steady-state saturated flow is solved by imposing mass conservation and the Darcy equation:

Figure 3 . 2 (Figure 3 .

 323 Figure 3. Spatial distribution of the solute concentration in RW model at times (a) 5 s, (b) 20, and (c) 100 s and at case Pe = 50. The spatial proportion between the x and y axis is 1:6.

Figure 4 .

 4 Figure 4. Time (step) evolution of global mixing state of Markov models using the computed mixing transition matrices 𝑀 𝑃 , 𝑀 𝑐 , and 𝑀 𝑤 from RW simulations at different time steps. Note that the evolutions computed with the three 𝑀 𝑤 are identical and identical to the evolution resulting from 𝑀 𝑝 and 𝑀 𝑐 computed after 250 RW time steps, when particles have exhaustively sampled the velocity distribution.

Figure 4 .

 4 Figure 4. Time (step) evolution of global mixing state of Markov models using the computed mixing transition matrices M P , M c , and M w from RW simulations at different time steps. Note that the evolutions computed with the three M w are identical and identical to the evolution resulting from M p and M c computed after 250 RW time steps, when particles have exhaustively sampled the velocity distribution.

Figure 5 .

 5 Figure 5. Temporal evolution comparison of (a) mean x position and (b) solute concentration variance in x between random walk and water parcel methods at σ 2 Y = 1 heterogeneity level.

Figure 5 .

 5 Figure 5. Temporal evolution comparison of (a) mean x position and (b) solute concentration variance in x between random walk and water parcel methods at 𝜎 𝑌 2 = 1 heterogeneity level.

Figure 6 .

 6 Figure 6. Temporal evolution comparison of the global mixing state between the random walk and water parcel methods at 𝜎 𝑌 2 = 1 heterogeneity level.

Figure 6 .

 6 Figure 6. Temporal evolution comparison of the global mixing state between the random walk and water parcel methods at σ 2 Y = 1 heterogeneity level.

Table 1 .

 1 Flow and transport problem parameters and simulation details.

	Flow			Transport
	λ (m)	10	Num. time steps	100	c re f (kg/m 3 )	1
	L x	600λ	φ	0.3	∆t (s)	1
	L y	150λ	RW		Cases
	∆x, ∆y	λ/10	N Particles	2.25 × 10 6	D w m 2 /s	Peclet
	µ Y	0	WP		0	∞
	σ 2 Y	1	N v N parcels	30 1.44 × 10 5	10 -2 0.2	10 3 50

Acknowledgments: J.S.-S. is indebted to KRATOS for providing a powerful tool and for help in developing the MAWMA application.

Data Availability Statement: Data used for producing the figures can be downloaded from digital.csic.es (https://digital.csic.es/handle/10261/245368 (accessed on 15 September 2021)).

Funding: This research was funded by MICINN (Spanish government) through MEDISTRAES project II and III (grant numbers CGL2013-48869-C2-2-R and PID2019-110212RB-C22). J.S.-S. wishes to acknowledge the financial support from the AGAUR (Generalitat de Catalunya, Spain) through the "Grant for universities and research centers for the recruitment of new research personnel (FI-DGR 2013)". IDAEA-CSIC is a Centre of Excellence Severo Ochoa (Spanish Ministry of Science and Innovation, Project CEX2018-000794-S).

Conflicts of Interest:

The authors declare no conflict of interest.