Modeling Mixing in Stratified Heterogeneous Media: The Role of Water Velocity Discretization in Phase Space Formulation

Modeling solute transport in heterogeneous porous media faces two challenges: scale dependence of dispersion and reproducing mixing separately from spreading. Both are crucial since real applications may require km scales whereas reactions, often controlled by mixing, may occur at the pore scale. Methods have been developed in response to these challenges, but none has satisfactorily characterized both processes. In this paper, we propose a formulation based on the Water Mixing Approach extended to account for velocity variability. Velocity is taken as an independent variable, so that concentration depends on time, space and velocity. Therefore, we term the formulation the Multi-Advective Water Mixing Approach. A new mixing term between velocity classes emerges in this formulation. We test it on Poiseuille’s stratified flow using the Water Parcel method. Results show high accuracy of the formulation in both dispersion and mixing. Moreover, the mixing process exhibits Markovianity in space even though it is modeled in time.


Introduction
Solute transport in homogeneous media is well reproduced by the advection-dispersion equation (ADE). However, this is not the case in real aquifers because they are heterogeneous (Le Borgne et al. 2008a, b;Gjetvaj et al. 2015;Willmann et al. 2008), which leads to a commonly observed nonequilibrium (Alcolea et al., 2008;Vogel et al. 2006). Observed transport is termed anomalous (i.e., non-Fickian). Anomalous transport is evidenced by the scale dependence of dispersivity (Lallemand-Barres and Peaudecerf 1978), or by tailing in concentration breakthrough curves (Valocchi 1985;Carrera 1993). But beyond conservative transport, accurate representation of anomalous transport is critical for simulating chemical reactions (Battiato et al. 2009;Sadhukhan et al. 2014;Scheibe et al. 2015;Soler-Sagarra et al. 2016;Tartakovsky et al. 2009). The main limitation of the ADE lies on not distinguishing between dispersion (solute spreading) and mixing (solute diffusion and dilution). The two processes are linked, but they are different (Dentz and Carrera 2007). In contrast to dispersion, mixing is a direct cause of chemical reactions (Cirpka and Valocchi 2007;Rezaei et al. 2005;De Simoni et al. 2005Soler-Sagarra et al. 2022;Tartakovsky et al. 2008;Cirpka 2002;Herrera et al. 2017). The ADE employs Fick's law (Fick 1855) to characterize both processes and must, therefore, be considered inadequate for reactive transport . A new formulation is needed to reproduce advection, dispersion and mixing (de Dreuzy et al. 2012;de Dreuzy and Carrera 2016). A large number of particle-based methods have been proposed as alternatives to the ADE (Benson and Meerschaert 2009;Bijeljic and Blunt 2006;Le Borgne et al. 2008a;Delay et al. 2005;Lester et al. 2014;Painter and Cvetkovic 2005;Russian et al. 2016;Schmidt et al. 2017;Sole-Mari et al. 2020). All of them are relevant to anomalous transport. But, while they have yielded new insights on transport, none of them considers mixing explicitly. An essential such insight is that velocity transitions after every step can be viewed as a correlated random process. This process is Markovian when transitions are made not after a fixed time step, but after particles have covered a fixed spatial distance (Le Borgne et al. 2008b). The fact that velocities may change after a fixed spatial step is consistent with a fixed heterogeneity structure. We conjecture that this is a good basis for alternative transport formulations.
The difficulty in representing mixing lies in its close relationship with spreading. Velocity variations produce stretching of lamellas, which enhances mixing by increasing the contact area between different waters (Le . The fact that velocity variations occur at all scales and that they control mixing suggest using velocity as a new dimension of the state variable (like time and space), which leads to a phase space formulation. That is, concentration at any representative volume will be given by a velocity-dependent distribution representing not so much uncertainty as actual variability (e.g., concentration at the leading edge of a plume will be larger at high-velocity paths than at low-velocity paths).
The success of Markovian formulations further suggests representing velocity variability as a Markov process. Markovian processes are typically represented by means of a transition matrix M vs p (or transition probability density for continuous representations of velocity), which expresses the probability, p, of a particle to change the velocity state v given constant steps in space phase s (De Anna et al. 2013;Kang et al. 2011Kang et al. , 2014Kang et al. , 2015Kang et al. , 2017. Transitions may occur either because of heterogeneity along a flow line, which do not produce mixing, or because of water particles diffusion between adjacent flow tubes, which is the mixing mechanism associated with plume stretching. A proper representation of mixing should distinguish these two types of transitions. In this paper, we propose a phase space formulation for transport that acknowledges velocity transitions driven by both heterogeneity along flowlines and by diffusion across. The solution method is based on the Water Mixing Approach (WMA) , because our ultimate motivation is reactive transport (RT).

Governing Equations
The ADE is the most widely used formulation for transport. The ADE expresses the solute mass balance when transport is driven by advection and dispersion.
where [-] is porosity, where the term q D c represents solute exchange driven by water dispersion and mixing. That is, q D [L 3 /L 2 /T] represents water flux exchange with respect to the mean water flux, which is accounted for in the advection term. Soler-Sagarra et al. (2022) discuss in detail the meaning of q D c . The challenge is how to write this exchange term so that it represents separately dispersion and mixing. Kang et al. (2017) proposed a phase space formulation for heterogeneous domains to find an alternative to ADE. Phase space formulations express state variables not only as dependent on time and space, but also on velocity. The formulation was originally presented for porescale models using particle probabilities, p. However, it can be easily extended to Darcy scale and written it in terms of concentrations, c = c (x, v, t) [M/L 3 ] by using basic definitions to write p = c ∕M , where M is the total solute mass. With these definitions, Kang et al. (2017) can be rewritten as where l [L] is the characteristic length, g vs [T/L] is the transition probability density of jumping from v ′ to v after a l space step and c � = c(x, v � , t) the concentration at v ′ . Thus, the formulation implies a distribution of (velocity-dependent) concentrations at every location and time. Equation (3) can be viewed as a mass balance equation, where the first term in the RHS expresses advective changes, the last term includes sink and sources, and the second and third terms represent mass losses and gains, respectively, due to velocity transitions (i.e., losses due to particles that transit from v to another velocity and the reverse gains). (1) This way, dispersion is explicitly represented in a natural way. The problem with this representation is that it does not distinguish between diffusion transitions (purple arrow in Fig. 1) and advection transitions (green arrow in Fig. 1). This is inappropriate since it does not allow treating mixing and dispersion as separate processes. To overcome this limitation, we propose to (a) restrict transitions caused by heterogeneity to advection transitions (i.e., changes in velocity along streamlines, like defined by green arrow in Fig. 1) which express dispersion explicitly, and (b) simulate diffusion separately. Furthermore, we use the WMA formulation proposed by  to define the mixing process, so as to facilitate generalization to reactive transport. It is important to notice that q D c exchanges in Eq.
(2) are now restricted to molecular diffusion exchanges, because the effect of velocity fluctuations is already included in Eq. (3)). Treating these exchanges as water exchanges may sound confusing since diffusion is commonly associated with solute fluxes, rather than water exchanges. In reality, water is exchanged by diffusion at a rate comparable to that of solutes (Harris and Woolf 1980). Furthermore, without entering into this debate (see Soler-Sagarra et al. (2022) for details), the formulation of Eq. (2) is equivalent to Fickian diffusion if q D = D w ∕L D , where D w [L 2 /T] and L D [L] are water molecular diffusion coefficient and the characteristic diffusion scale, respectively. We are writing Eq. (3) per unit volume of water, which is more convenient than per unit volume of medium (as in Eq. (2)). Therefore, we will write these exchanges as v Dm c , where the subscript emphasizes that we only include diffusive exchanges.
The new formulation is obtained by restricting Eq.
(2) as to diffusion term Equation (4) is not a complete formulation yet. Velocity transitions occur due to advection and diffusion ( Fig. 1). Since we are still restricting g vs to characterize velocity transitions along streamlines, we need a new transition term for velocity changes driven by diffusion (purple arrow in Fig. 1). Advection changes characterized by g vs are Markovian in space (Le Borgne et al. 2008b), but diffusion driven exchanges should be Markovian in time. The two transitions must be described independently. Diffusion transitions require adding two terms to (3), which yields where f vt v|v ′ [T/L/T] is the probability density of the rate of mixing transitions between velocity states. The fourth term on the RHS define the diffusion process in space domain (orange arrow in Fig. 1), while the fifth and sixth expresses the diffusive mass balance in velocity domain (purple arrow in Fig. 1). Since we are assuming diffusion to be independent of the solute, then f vt is also solute independent, so that this equation is a WMA equation. Therefore, we term this formulation Multi-Advective Water Mixing Approach (MAWMA). Note that we consider the requirements highlighted by De Dreuzy and Carrera (2016): adequate separation of advection, diffusion and dispersion.
Assessing the validity of the MAWMA, Eq. (5), could be arduous. Given that the novel processes presented are the fifth and sixth terms on the RHS, we make three simplifications for testing: (a) Flow is assumed stratified. This implies that no velocity transition occurs due to advection, which leads to neglecting the second and third terms on the RHS of the Eq. (5).
We further assume that all strata carry the same flow rate (i.e., high-velocity strata are narrower than low-velocity strata) to simplify space and velocity discretization ( Fig. 2b)  (Bolster et al. 2011;Dentz and Carrera 2007;Taylor 1953). The fourth term of the RHS may therefore be ignored. Transverse mixing has been proven to be of paramount importance when chemical reactions are involved (Werth, et al. 2006) (c) A Lagrangian formulation is adopted for advection by using material derivative d ⋅ ∕dt to minimize numerical dispersion and mixing (Batlle et al. 2002;Bell and Binning 2004;Cirpka et al. 1999;Ramasomanana et al. 2012;Soler-Sagarra et al. 2022;Zhang et al. 2007).
These three simplifications allow us to rewrite Eq. (5) as Time, space and velocity should be discretized and integrated in Eq. (6).

Methodology
This section describes the method to solve Eq. (6) and to test the solution by comparison with existing methods. We first describe how to discretize velocity and how to build an isochronal mesh. Then, we describe the RW and Isochronal Water (IW) methods, which are used for comparison purposes. We finally describe the proposed solution approach and how to compute the velocity transition matrix M vt w for each method that will be used by the Water Parcel (WP) method based on MAWMA.

Water Velocity Discretization
Velocity in the domain may be discretized adopting Eulerian or Lagrangian distribution. The Eulerian distribution yields the probability density function (pdf E ) of velocity sampled randomly in space. Discretizing pdf E with equal Eulerian probability velocity classes yields classes with the same water volume. The Lagrangian distribution pdf L is obtained by sampling the velocity equidistantly along streamlines (s-Lagrangian velocities according to Dentz et al. (2016)). Even in a stationary flux field, the pdf E and pdf L might be dynamic for a set of solute particles that are transported (Puyguiraud et al. 2019). However, in the Water Mixing concept, the water is transported instead of solute. As a result, pdf E and pdf L are stationaries in a stationary flow field. Discretizing pdf L determines velocity classes with the same flux (i.e., same injection probability, see Fig. 2c). The pdf L is related to pdf E as  where v mean is the mean of the pdf E . Discretizing velocities in classes with the same flux implies that low velocities are more probable than high velocities in the domain (Gotovac et al. 2009), which is realistic but contrasts sharply with Eulerian equi-probable discretization. Lagrangian discretization is the appropriate one in our case because Eq. (6) is written in Lagrangian form. We discretize velocities in N v classes, whose boundaries are: where cdf −1 L is the inverse of the Lagrangian cumulative distribution function. Once the boundaries are set, a representative class velocity v m is calculated as the mean of this interval as: The probability p m of the velocity class m in the domain (i.e., volumetric fraction of domain occupied by the class) is: Note that p m is inversely proportional to v m because the classes were chosen with identical Lagrangian probability, whereas the Eulerian pdf is inversely proportional to velocity [Eq. (7)].

Isochronal Mesh
The term isochronal mesh stands for a streamline-oriented mesh that advects the entire volume of water in a cell (we term it parcel) to the cell downstream within the stream tube. The stream tubes widths are proportional to the velocity classes probabilities, p m . The length ∆x m of the cells in the stream tube associated with velocity class m is calculated as: where Δt is the time step. Note that the mesh is built once the simulation time step is known unlike most meshes. Figure 2b shows a scheme of the isochronal mesh. The velocity classes are equally flux because a single stream tube is associated with each class. Moreover, all the water parcel of this isochronal mesh has the same volume because as the probability p m is inversely proportional to its velocity v m (see Fig. 2b).

Random Walk Method
A number of particles N p are injected in the domain in a time step Δt . The particles are flux weighted distributed along the domain width a. The solute mass m p associated with each particle is calculated as: where c inj is the concentration of injection. Note that the numerator in the right-hand side of Eq. (13) is the total mass injected during a time step. Each particle is advected with the To calculate concentrations, we use the isochronal mesh described in Sect. 3.2. The concentration of a cell i at time k is given by where N pi is the number of particles inside the cell and V w is its water volume.

Isochronal Water Model (WMA)
The isochronal structured water (IW) method is based on the Lagrangian form of the WMA (Eq. (2) with the advection term moved to the left-hand side so as to get the material derivative) and uses the isochronal mesh detailed in Sect. 3.2. The mesh is discretized in cells with the same water volume (parcel). Once the spatial discretization is set, Eq. (2) can be written in matrix form like: where S is the storage matrix, c k A is the vector of concentrations at all parcels after advection, θ is the time weight parameter ( 0 ≤ θ ≤ 1 ), and P [1/T] is the transition rate matrix that in Eq. (2) accounts for both mixing and spreading. However, here it only accounts for mixing (diffusion) because of the stratified flow conditions. The concentration can be expressed as: where M st w = (S − θΔtP) −1 (S + (1 − θ)ΔtP) is the transition matrix of water volumes in space s at time phase t. The vector c k A is obtained for every cell i, considering the isochronal property of the mesh which advects the entire parcel of the upstream cell i-, as: Note that if all species have the same transport parameters, Eqs. (18) and (19) can be viewed as the transport of water parcels instead of individual solutes and the concentrations would be attributes of these parcels. Soler-Sagarra et al. (2022) demonstrate the immediate extension from this scheme to a reactive transport method.
To obtain the transition matrix M vt w that WP method uses, Eqs. (17) and (18) are performed in velocity phase v instead of space phase s. The later implies that the diagonal of S contains the probability of each velocity class and P the transition rate between classes. In the stratified case, it can be viewed as any entire stream tube of the Fig. 2b is considered as a single element.

Water Parcel Model (MAWMA)
We now describe the Water Parcel (WP) method used to solve Eq. (6). Other methods might also be used. The spatial domain is discretized in parcels ( Fig. 2c) with the same water volume as the IW method (Fig. 2b). The parcel discretization covers the entire water domain, in saturated conditions (as it is the case). Like in IW, the concentration is only considered an attribute of each parcel. Each water parcel is associated with a centroid that determines its position in both x axis and velocity state (see Fig. 2c). Centroids are injected and displaced through the domain like a single solute particle.
The velocity state of each parcel is assigned randomly after injection (as discussed in Sect. 3.1, all classes are equally probable) and leads to an unstructured mesh unlike IW (see Fig. 2b, c). We integrate along the y coordinate for simplicity and for demonstration purposes. That is, we perform a dimension reduction, so that concentration in Eq. (6) depends solely on x and v. As shown in Fig. 2, this simplification might look trivial as it suggests that we are substituting the y coordinate by v. Note, however, that we assume that we do not know the vertical structure of velocity, but only its velocity distribution and transition probabilities. Figure 2c shows the parcel shape dependence on velocity state. As suggested by the tub lines of IW (see Fig. 2b), the longitudinal axis of our water parcels is proportional to their velocities, while their width is inversely proportional (see Fig. 2c). Another explanation is that the distance traveled Δx is proportional to the velocity v at the same time step Δt (explained in Sect. 3.4). As a consequence, the water parcels with low velocity tend to cram longitudinally (i.e., number of low-velocity water parcels per unit length is inversely proportional to velocity). This ensures an adequate representation of the entire distribution of velocities.
The simulation proceeds by integrating Eq. (6) in time steps. Therefore, the question is how to reproduce mixing between the parcels (i.e., the first two terms in the RHS). We use the finite volume method. The discretized form of (6) for every parcel i in velocity class l is where N mi is the number of parcels of velocity class m connected to parcel i. F vt lm is the volume of water exchanged between velocity classes l and m per unit time. k is time steps number and I l is the domain associated with the velocity class l. Finally, a ij is the fraction of this flux that will be exchanged between parcels i and j which could be either equidistributed or weighted by contact area. The latter is assumed here. As in WMA, this is an exchange process, which implies that F vt lm = F vt ml and a ij = a ji . The expression of a concentration in time step k + 1 can be obtained.
Two remarks can be made regarding this expression: (a) lm = Δt ⋅ F vt lm ∕V wl is a water mass mixing ratio, similar to the one of Soler-Sagarra et al. (2022)  Note that the mixing of parcels depends on their velocity class. However, the unstructured mesh does not ensure mass conservation in the mixing process such as ij = ji , which diminishes exchanges of water volume. This is why mass conservation is imposed after a first calculation of lambdas by defining ij = ji = max ij , ji .
Although chemical reactions are not the objective of this work, they are the ultimate goal of our research. By using the water mass mixing ratio formulation, the link with chemical processes is immediate .
The proposed WP method uses the transition matrix M vt w for just diffusion. Le Borgne et al. (2008b) demonstrated Markovianity in space using a similar transition matrix, M st p , which accounted for both advection and diffusion. Since diffusion is Markovian in time, the question can be raised if their conclusion was an artefact of small diffusion or M vt w is also Markovian in space, which needs to be tested. The computation of these matrices is quite intensive. We can only hope that they can eventually be parametrized as a function of geological understanding and flow regime. But this requires much work beyond fully ascertaining the validity of the proposed method.

Applications
The efficiency of the MAWMA formulation, solved with the WP method described in Sect. 3.5, is tested by comparison with the RW and IW methods, described in Sects. 3.3 and 3.4, respectively. Comparisons are made in terms of mixing and spreading indicators in Sect. 4.1, where we also test whether WP results are sensitive to the way matrix M vt w was computed (from RW or IW). In Sect. 4.2, the Markovianity in space is tested by comparing where v mean is the mean velocity of the distribution, y is the vertical position and a is the half distance between the planes. Owing to the horizontal symmetry of the case, we only model the half domain y = {0, a} (Fig. 2). We focus on the time evolution, especially at times earlier than the dispersion time scale τ D , which denotes the typical time for the macrodispersive spreading of the solute.
Moreover, we considered a continuous injection of solute instead of an instantaneous injection used by other authors (Bolster et al. 2011;Dentz and Carrera 2007). Simulation details are shown in Table 1. The WP was simulated using the KRATOS framework (Dadvand et al. 2010). The RW simulations were performed as proposed by Dentz and Carrera (2007), using a flux weighted injection at every time step.

Statistical Indicators of Mixing and Spreading
Mixing and spreading must be tested. Mixing can be assessed by the global mixing rate (i.e., average value of the mixing factor, ∇ t cD∇c ). However, Le Borgne et al. (2010) proposed a far simpler and more robust and stable quantification by means of the scalar dissipation rate (Pope 2000), which measures the time derivative of the concentration variance, but requires no solute flux at the boundaries. We extend in here (Appendix A) the continuous injection case, where Ω is the simulation domain, Γ is its boundary and c inj is the injection concentration. The global mixing rates computed with the three models are plotted in Fig. 3a. Only a slight mismatch is observed at the earliest times. Some oscillations are observed at early 3 V w parcel 3.3 · 10 -5 (m water ) 3 V w parcel 2 · 10 -4 (m water ) 3 times for the RW and WP solutions. These oscillations may be attributed to the fact that concentrations are calculated from the particle positions (Fig. 2a) in the IW mesh (Fig. 2b).
Despite oscillations, the same overall behavior is displayed by all the models. The global mixing rate typically displays a diminishing monotonous behavior for Dirac delta injections (Bolster et al. 2011;Le Borgne et al. 2010). In our case, we observe an initial increase, resulting in a bell shape which peaks at 0.1 D . Thus, two regimes are distinguished. At early time, mixing is enhanced in response to the increase in longitudinal dispersion, which  (Haber and Mauris 1988) and the asymptotic dispersion (Aris 1956), respectively stretches the contact area between the invading and the resident waters. This stretching phenomenon has already been observed in instantaneous injection (Le Borgne et al. 2013. However, it is the continuous injection what causes the dissipation rate to increase. After the peak, mixing decreases as the front approximates macrodispersion regime because the concentration becomes smooth. A similar behavior can be observed for finite size initial conditions (Kapoor and Kitanidis 1998).
As for dispersion, the adopted continuous injection makes it inappropriate the traditional definition (rate of growth of the second spatial moment). Therefore, we computed the spatial variance of the concentration gradient distribution instead of the concentration distribution. Other methods could have been used including a fit to the 1D analytical solution or correcting the spatial second moment with the uniform injection. We integrate vertically the concentration to obtain its correlation with the x coordinate. The (apparent) dispersion D app is where 2 ∇c is the variance of the vertically integrated concentration gradient. The analytical solution of the temporal dispersion evolution D a (Haber and Mauris 1988) and its asymptotic value D asy (Aris 1956) are also computed.
The results are plotted in Fig. 3b. Although the WP dispersion oscillates (owing to the unstructured character of the mesh) a satisfactory agreement is again observed. As in the dissipation rate, at least two different regimes of the solute distribution may be distinguished: (a) a linear increase in the variance is observed. This confirms the stretching phenomenon described above; (b) an asymptotic macrodispersion regime is attained close to 0.1 D . The transition regime roughly coincides with the scalar dissipation rate, suggesting a link between both behaviors. Indeed, spreading enhances mixing in the earlier regime. In the later regime, solute plume extension is limited since sufficient mixing occurs.

Markovianity in Space
Although solute only changes its velocity class because of the mixing process (which is Markovian in time), we can calculate the transition matrix in space M vt p (Le Borgne et al. 2008b) from the particle RW model. We believe that they are also Markovian in space, which is consistent with (Le Borgne et al. 2008b). We tested the Markovianity by comparing the transition probabilities with the ones obtained from a Markov chain model. The transition model must satisfy the Chapman-Kolmogorov equation (Risken 1996), which reads for the transition matrices M(x) of a discrete Markov chain such as with x, Δx > 0. The latter implies Although the agreement is not exact (Fig. 4), the particles satisfactorily reproduce the Markov chain results. We can therefore conclude that mixing is Markovian not only in time, but also in space. This result suggests that we can calculate transport using time steps simulations, which are more adequate to model mixing.

Conclusions
We present a new formulation for solute transport in heterogeneous cases, termed MAWMA. The formulation aims to reproduce diffusive mixing and dispersion. The formulation is an extension of WMA by making the transport state dependent on velocity as well as on time and space.
Water parcel models were employed for numerical solution of the proposed equation. Each parcel was associated with its centroid, which defines its velocity and position at a given time. We tested the velocity transition produced by mixing applying a water transition matrix in time M vt w . This differs from the solute transition matrix in space M vs w used in the correlated CTRW model.
The formulation was tested on a Poiseuille's stratified flow case. MAWMA was compared to WMA and Random Walk methods in terms of global mixing and dispersion. A good agreement was observed. Moreover, mixing shows Markovianity in space even when it is modeled with constant time steps. The results suggest that MAWMA will perform well for high heterogeneity cases using a matrix M st w for advection transitions.

Appendix A: Dissipation Rate in Continuum Injection
Several works have contributed to the use of the scalar dissipation rate to measure global mixing (Pope 2000;Le Borgne et al. 2010;Hidalgo et al. 2012;Jha et al. 2011;Nicolaides et al. 2015). To address the continuous injection case, we start with the basic definition where n is the unit vector normal to boundary Γ. We substitute the second term of the RHS of equation (A2) by using the ADE ∇ ⋅ (D∇c) = ( c∕ t) + q∇c We can employ again Green's first identity to the third term of the RHS We use the flow equation ∕ t = −∇q in the fourth term of the RHS and regroup the equation We now consider the inlet boundary condition qc inj = −(qc − D∇c) ⋅ n| Γ Note that n has a sign opposite to that of the flux because the inlet boundary faces backwards. This is why q = −q ⋅ n| Γ . Given the latter and assuming that porosity is constant, we end up with