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

: 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 deﬁned 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. 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.

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 [28][29][30][31][32][33][34][35]. 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 [36]. Thus, a transition matrix M vs is needed [18,[37][38][39][40]. M vs stands for the 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, 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 [15,49]) 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). [41]) 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 [15,49]) 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). [41]) 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 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, [49,97] proposed expressing diffusive transitions in velocity (purple arrow in Figure 1) as 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 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:

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 [49,97]. 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 (7). For simplicity, the discretization procedure is the one used by [49,97]. 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 [98]).
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 [15,49]. 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 [98]). 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. 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. 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.
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 [18,24,35,37,39,41] 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 is the time step, Nli is the number of parcels with velocity l spatially connected to i, F is the water volumetric flux diffused, 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.
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 [18,24,35,37,39,41] 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 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 [15], which leads to 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, 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. [15,49] 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 [99]. 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 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: 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 (11) 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., [15]. 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, 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 [35] 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 [98]) and velocity class concentration is expressed as 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 Combining (12) and (13) Sc k+1 = M p Sc k Then, the concentrations for the next time step are as follows 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 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 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 where M V is the volume exchange matrix computed initially from Equation (17). Finally, the water transition matrix is expressed as 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 [48]. 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) 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 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 The WP method should reproduce mean advection, dispersion, mixing, and "be flexible enough to be applicable to real problems" [27]. 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 [15,49]. Still, a number of developments are needed to address the real cases with a level of maturity comparable to stochastic methods [32,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 [15,49]. Still, a num of developments are needed to address the real cases with a level of maturity compara to stochastic methods [32,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 where m (k) x is the k-th order moment of the solute distribution in space 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 Global mixing G [105] is defined as Note that we can also define global mixing G in terms of the velocity domain such as 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 (12). The initial solute probability distribution for any velocity class i is 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 (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: , , 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 (12). 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 and 6. 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 [15]).
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 [15]).
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 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 [15]). 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.