Upscaling Mixing in Highly Heterogeneous Porous Media via a Spatial Markov Model

Abstract

In this work, we develop a novel Lagrangian model able to predict solute mixing in heterogeneous porous media. The Spatial Markov model has previously been used to predict effective mean conservative transport in flows through heterogeneous porous media. In predicting effective measures of mixing on larger scales, knowledge of only the mean transport is insufficient. Mixing is a small scale process driven by diffusion and the deformation of a plume by a non-uniform flow. In order to capture these small scale processes that are associated with mixing, the upscaled Spatial Markov model must be extended in such a way that it can adequately represent fluctuations in concentration. To address this problem, we develop downscaling procedures within the upscaled model to predict measures of mixing and dilution of a solute moving through an idealized heterogeneous porous medium. The upscaled model results are compared to measurements from a fully resolved simulation and found to be in good agreement.

1. Introduction

Mixing is the process that brings dissolved chemical species together. Thus, accurately accounting for mixing is particularly important in the context of correctly predicting mixing-driven chemical reactions [1,2,3,4]. Many existing reactive transport modeling approaches assume perfect mixing at some scale, while in reality incomplete mixing of reactive species occurs below this scale. The effects of this have been observed in theory [5,6,7,8], numerical modeling [5,6,7,9,10], laboratory experiments [11,12,13,14,15], and field studies [16,17,18]. Incomplete mixing typically results in an overestimation of the amount of reaction that will occur, presenting the need to artificially, and often non-physically, alter the effective reaction rate used in the model to better match observations [19]. This has led to the development of upscaled models that aim to account for subscale concentration fluctuations that limit reaction e.g., [20,21,22]. The correct prediction of reaction rates has many practical implications in the context of porous media and aquifers, e.g., the prediction of contaminant migration [23,24], the remediation of contaminated groundwater [25,26], and the fundamental prediction of naturally occurring geochemical reactions that shape the subsurface below us. For these reasons, the development of models that accurately describe mixing behaviors are necessary.

In general, real flows through geologic porous media are non-uniform with heterogeneity present over a wide range of scales [2]. For many such systems, it is simply not computationally feasible to resolve the full flow heterogeneity across all of these scales [2]. This presents the need for upscaled models that capture the effects of these small scale processes without resolving them. In highly heterogeneous porous media, the complex flow organization is characterized by variations over several orders of magnitude in the hydraulic conductivity field and strong channeling of high velocity regions [27,28]. This results in larger Lagrangian velocity correlations in space for particles in higher velocity zones than in lower velocities [27]. It can also lead to strong Lagrangian accelerations as particles transition between high and low velocity regions [28]. To accurately model mixing in highly heterogeneous systems, these flow features must be taken into account.

Mixing generally happens more quickly in non-uniform systems as a result of the stretching and compressing experienced by a fluid due to flow heterogeneity [29,30,31]. In advection-dominated transport settings, the underlying flow structure will have a large impact on transport and mixing in a system. The Péclet number (Pe), which quantifies the ratio between advective and diffusive processes, plays a key role in determining the enhancement of mixing due to the non-uniform flow. In high Pe systems with strong velocity contrasts, the flow will deform a solute plume and significantly impact mixing. This leads to very different behaviors if compared with those emerging in homogeneous media. This has a clear impact on reaction dynamics as investigated in detail by [32]. Flow topology itself can also have a strong impact on mixing. This is exemplified by recent studies which observed that regions of chaotic mixing can occur in Darcy-scale groundwater flows under transient forcing conditions [33] and examined the role of twisting streamlines on the enhancement of mixing in three-dimensional anisotropic porous media [34].

There are many ways to quantify mixing, ranging from local to global scales. Many studies consider global measures of mixing such as the scalar dissipation rate [1,35,36] and the dilution index [30,35,37]. These global metrics consider the amount of variation in the concentration field over the entire flow domain as it evolves in time. At local scales, metrics relating to plume deformation induced by velocity gradients have been shown to be useful. These include, for example, the Okubo-Weiss parameter [30,38], the largest eigenvalue of the Cauchy-Green strain tensor [38], the finite-time Lyapunov exponent [38], and the trace of the local strain matrix squared [39]. Recently, strong mixing regions identified by these local metrics of flow deformation have been linked to higher amounts of reaction [40]. The Okubo-Weiss parameter [30] and the trace of the local strain matrix squared [39] have also been explicitly linked to the rate of evolution of the the dilution index over time, demonstrating a correlation between global mixing behavior and local flow deformation.

Many approaches exist for modeling effective upscaled transport in highly heterogeneous porous media, e.g., multi-rate mass transfer models [41], fractional advection-dispersion equations [10], and continuous time random walks [42]. In this paper, we will focus on one particular model, the Spatial Markov model (SMM). The SMM has been used to predict effective mean transport in a wide variety of flowing hydrologic systems, including highly heterogeneous geologic media [28,43], fracture media [44,45,46,47,48,49,50], and complex pore scale systems [51]. It has also been extended to predict first-order and bi-molecular reactions in relatively simple idealized pore scale systems [52,53]. The SMM belongs to the broader family of continuous time random walk (CTRW) models. In its discrete implementation, it imposes correlation between successive particle steps. This is unlike many other approaches that assume statistical independence between subsequent jumps. By doing so, one can typically capture upscaled behaviors over smaller scales than an uncorrelated model allows. It has been shown that accounting for such correlation is particularly important in advection-dominated, steady flows [54,55].

The SMM, as we will apply it, is a reduced-dimensionality 1-d model that can describe the location of a solute plume over time and thus predict mean concentrations. While this is useful for predicting common depth-averaged observables such as breakthrough curves, having information on only the mean transport is insufficient to predict mixing. Since mixing is a small scale process associated with diffusion and the deformation of a plume by a non-uniform flow, fluctuations in concentration must also be captured in a model, so as to accurately predict large scale mixing dynamics. Therefore, improvements and adjustments must be made to the traditional SMM in order to successfully upscale mixing. Recently, the SMM was extended to predict mixing using a trajectory-based method in a periodic flow domain [56], where the periodic nature of the geometry and flow were exploited to downscale the full concentration field from mean transport predictions. Here, we propose a conceptually similar approach by examining transport through a highly heterogeneous porous medium. Given the lack of periodicity, the downscaling approach is not as intuitive or obvious. For this reason, we consider a variety of approaches to predict effective mixing using the SMM. To do this, we develop several downscaling methods within the upscaled SMM to downscale solute particle locations and reconstruct representative concentration fields that can be used to quantify mixing.

2. System Setup

The velocity field $\mathbf{v}\left(\mathbf{x}\right)$ considered in this work is representative of an idealized two-dimensional heterogeneous porous medium. To obtain the flow field, a log-normal random Gaussian-correlated permeability field $\kappa (x,y)$ (dimensions of $\left[L^2\right]$) is generated with zero mean and variance ${\sigma }_{ln\kappa }^2=9$. The permeability is discretized in squared grid cells of unit length and we fix correlation length $\lambda =2$. The natural log of the permeability field used in this study is depicted in Figure 1a. The distribution, variance, and correlation length of the permeability field were chosen to align with the field used in the work of Le Borgne et al. [28,43], which first introduced the SMM and applied it to successfully predict mean transport characteristics, including the rate of spreading of the plume as quantified by the second centered moment as well as breakthrough curves. After generating our permeability field, we solve Darcy’s law coupled with incompressibility

$$\mathbf{v}\left(\mathbf{x}\right)=\frac{\kappa \left(\mathbf{x}\right)}{\mu }\nabla p,\nabla ·\mathbf{v}=0,$$

(1)

to obtain our velocity fields, where v is the velocity $\left[L{T}^{-1}\right]$, $p\left[M{L}^{-1}{T}^{-2}\right]$ is pressure, and $\mu \left[M{L}^{-1}{T}^{-1}\right]$ is viscosity. System (1) is solved using a finite volume method [57]. At the boundaries we impose permeameter like conditions: no flux across the lateral boundaries and constant head at the upstream and downstream boundaries. The mean velocity can then be rescaled to any desired value, so for convenience we choose $\overline{v_x}=1$. Figure 1b shows the natural log of the absolute value of the velocity $v=\sqrt{v_x^2+v_y^2}$, where $v_x$ and $v_y$ are the velocity fields in the horizontal and vertical directions, respectively. The flow domain was generated to be large enough so that when we simulate transport the solute will not interact with the boundaries. Domain lengths in the x and y directions equal to $L_x=8000$ and $L_y=2000$, respectively, were deemed sufficient for this purpose.

The governing equation for transport is the advection dispersion equation

$$\frac{\partial C(\mathbf{x},t)}{\partial t}=D{\nabla }^{2}C(\mathbf{x},t)-\mathbf{v}·\nabla C(\mathbf{x},t),$$

(2)

where C is the solute concentration $\left[M{L}^{-2}\right]$ and D is the constant isotropic dispersion coefficient $\left[L^2{T}^{-1}\right]$. Our dispersion coefficient is selected to be a constant for simplicity and to align with the work of Le Borgne et al. [28,43], although it should be noted that using a velocity-dependent dispersion coefficient can have an impact on mixing [58,59]. Rather than solving Equation (2) directly, we model it with an equivalent particle tracking random walk approach. Our plume is discretized into a large number of particles. Each particle’s mass is given by $m_p=m_{tot}/N_p$, where $m_{tot}$ is the total amount of mass of the solute, $N_p$ is the number of particles, and $m_p$ is the mass carried by each particle. $m_{tot}=1$, $N_p={10}^5$, and $m_p={10}^{-5}$ are chosen for this study. In this fully resolved model, the particles move in time according to Langevin equation

$$\mathbf{x}\left(t_{n+1}\right)=\mathbf{x}\left(t_n\right)+\mathbf{v}\left(t_n\right)\Delta t+\sqrt{2D\Delta t}\mathbf{\xi }$$

(3)

where $\mathbf{\xi }\sim N(0,1)$, $\Delta t$ is the time step, and $\mathbf{v}$ is the particle velocity interpolated from the velocity field shown in Figure 1b. At every time step, the particle moves with an advective step and a random jump reflecting diffusion with zero mean and a variance equal to $2D\Delta t$.

We consider an advection-dominated flow with Péclet number equal to $Pe=\frac{\overline{v_x}\lambda }{D}=200$, where $\overline{v_x}=1$ is the mean velocity in the horizontal direction of the flow, $\lambda =2$ is the correlation length of the permeability field, and $D={10}^{-2}$ is the constant dispersion coefficient. Our simulation time step is $\Delta t={10}^{-3}$. For our initial condition, we consider a flux-weighted line injection of a solute. The line injection is placed at $x=0.4L_x$ between $y_{min}=0.25L_y$ and $y_{max}=0.75L_y$, as indicated on Figure 1. This location was chosen so that the solute would not interact with the flow boundaries over the times that we are interested in measuring mixing.

To quantify mixing in this system, we consider two commonly used metrics of global mixing. First, we measure the integral of the solute concentration squared, given by

$$M\left(t\right)=\int C^2(\mathbf{x},t)d\mathbf{x}.$$

(4)

$M\left(t\right)$ describes how fluctuations in concentration evolve in time over the entire domain. The time derivative of this quantity multiplied by $-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ is known as the scalar dissipation rate and is a well-known metric used to quantify the rate of mixing [1,35,36]. Our second metric of mixing is the dilution index, defined as

$$E\left(t\right)=exp\left(-\int C(\mathbf{x},t)\ln [C(\mathbf{x},t)]d\mathbf{x}\right).$$

(5)

The dilution index describes the degree of solute mass uniformity (i.e., dilution) in a system and is proportional to the volume occupied by a solute [37]. Under perfect mixing when $C(\mathbf{x},t)$ is uniform over the entire domain under consideration, the integral of the squared concentration $M\left(t\right)$ will be at a minimum, the scalar dissipation rate will be equal to zero, and the dilution index $E\left(t\right)$ will be at a maximum.

In order to compute $M\left(t\right)$ and $E\left(t\right)$, the concentration field must first be calculated. To obtain the concentration field, we map particles onto the same spatial grid as the velocity field. Then, the concentration in the ith grid cell is given by $C_i=N_i{m}_p/A_{grid}$, where $N_i$ is the number of particles in the grid cell and $A_{grid}$ is the area of the grid cell. Figure 2 and Figure 3 show the calculation of $M\left(t\right)$ and $E\left(t\right)$ in time from the fully resolved simulation using Equations (4) and (5), respectively (solid blue line). In the development of our upscaled model, our goal is to match the simulation results of this fully resolved simulation as well as possible.

3. Upscaled Model: Spatial Markov model

Next we describe the SMM and suite of modifications to the traditional method that were made for the purpose of predicting mixing. We provide a detailed overview of the SMM, but for more details refer the interested reader to the detailed descriptions provided in [54,60].

3.1. Model Parameterization

As in previous implementations of the SMM, we parameterize the model by running high resolution particle tracking simulations over two representative cells in our flow field using the random walk method described in the previous section using Equation (3). We initially have our flux weighted line injection located at $x=0.4L_x$ and for our model parameterization we track the evolution of each of these particles in time until they have reached the downstream location $x=0.4L_x+2L_{cell}$. A cell length of $L_{cell}=24\lambda$ was deemed appropriate for our system. Further discussion on this choice of cell length can be found in Section 3.2.1.

During the parameterization step, we record several features of a particle’s trajectory: (i) the particles’ y position at the inlet of the first cell ($y_0$), (ii) the particles’ y position at the inlet of the second cell ($y_1$), (iii) the time it takes for the particles to travel through the first cell (${\tau }_1$), and (iv) the time it takes for the particle to travel through the second cell (${\tau }_2$). This process is illustrated in Figure 4. After the joint distribution $f(y_0,y_1,{\tau }_1,{\tau }_2)$ is obtained, it is used to inform the upscaled model. Typically, only ${\tau }_1$ and ${\tau }_2$ are recorded in the parameterization of the Spatial Markov model, but the additional information on $y_0$ and $y_1$ is required here for our downscaling procedure, which aims to predict mixing.

3.2. Model Mean Longitudinal Transport

The SMM is a random walk particle-based model. As for the fully resolved model described previously, the plume is discretized into a large number of particles of equal mass $m_p$. Unlike the fully resolved model, a particle takes a uniform step in space of size $L_{cell}$ and moves forward in time a random amount $\tau$ at each step in the upscaled model. In the SMM, during each model step, particles move forward in time and space according to Langevin equation

$$\begin{array}{cc}\hfill x^{(n+1)}& =x^{\left(n\right)}+L_{cell}\hfill \\ \hfill t^{(n+1)}& =t^{\left(n\right)}+{\tau }^{(n+1)}n=0,1,2,\dots \hfill \end{array}$$

(6)

where ${\tau }^{(n+1)}$ is sampled from

$$f\left(\tau \right)=\left\{\begin{array}{cc}f\left({\tau }_1\right)\hfill & \mathrm{if}n=0\hfill \\ f\left({\tau }^{(n+1)}\right|{\tau }^{\left(n\right)})\hfill & \mathrm{if}n=1,2,\dots .\hfill \end{array}\right.$$

(7)

The distribution of ${\tau }_1$ ($f\left({\tau }_1\right)$) is measured directly during the parameterization step and the distribution of ${\tau }_2$ given ${\tau }_1$ ($f\left({\tau }_2\right|{\tau }_1)$) is modeled by a transition matrix. To obtain the transition matrix, $f\left({\tau }_1\right)$ is separated into 20 equiprobable bins and the cutoff times associated with each bin are recorded. This bin number is chosen based on a convergence test and is in line with previous estimates for the required number of bins [60]. Bin 1 contains the particles with the fastest travel times and Bin 20 contains particles with the slowest travel times. A particle with travel time ${\tau }_p$ is in Bin i if $t_{c,i}\le {\tau }_p<t_{c,i+1}$, where $t_{c,i}$ is the cutoff time for Bin i, $t_{c,1}=0$, and $t_{c,21}$ is greater than the maximum value of ${\tau }_1$ and ${\tau }_2$ for all particles. Then, the transition matrix is defined by

$$T_{i,j}=P({\tau }_2\in \mathrm{Bin}j|{\tau }_1\in \mathrm{Bin}i)\approx f\left({\tau }_2\right|{\tau }_1).$$

(8)

It is assumed that $f\left({\tau }^{(n+1)}\right|{\tau }^{\left(n\right)})=f\left({\tau }_2\right|{\tau }_1)$. Thus, each block of the transition matrix, $T_{i,j}$, describes the probability that a particle will have a travel time in Bin j given that its travel time was in Bin i in the previous step.

Using Equation (6) in the SMM, we know when a particle is at a location $N{L}_{cell}$, where N is an integer. In other words, at a specified time $t^{\ast }$, the SMM provides information about which cell each particle is in, how long it has been there, and how long it will remain there, but the exact location of a given particle beyond this is unknown. In order to predict mixing by estimating $M\left(t\right)$ or $E\left(t\right)$ as defined in Equations (4) and (5), we need to have an estimate of the spatial distribution of the concentration field at a given moment in time. To obtain such information from the SMM, we must make an educated guess as to where a particle will be located within a cell. This means that we must develop a downscaling procedure within the upscaled model to predict the location of each particle. Note that we are not attempting to accurately reproduce the detailed concentration fields observed in the fully resolved model, but are instead trying to develop a downscaling procedure that provides an accurate estimate of the integral of the squared concentration, M, and the dilution index, E, in time—i.e., we are trying to generate a numerical closure via a downscaling process. In Section 4, we describe a variety of different downscaling methods considered here for the purpose of predicting mixing. While we actually implemented several more than those described here, we focus on this limited number to highlight specific features.

3.2.1. Choice of Cell Length

To determine the appropriate cell length, breakthrough curves (BTCs) were measured using the upscaled SMM at four downstream locations at distances 48λ, 96λ, 144λ, and 192λ from the inlet and compared to BTCs measured with the fully resolved model. The SMM was run using a variety of cell lengths in order to determine a cell length that would be appropriate for use in our expanded SMM to measure mixing. The BTCs measured using the SMM with cell lengths of $L_{cell}$ = 6λ, 12λ, and 24λ (dotted lines) are shown on Figure 5 along with the corresponding BTCs measured with the fully resolved model (solid lines). The upscaled model should be able accurately capture the peaks and tails of the BTCs. By examination of Figure 5, it is clear that the upscaled model with $L_{cell}=6\lambda$ is unable to capture either the peaks or the tails of the BTCs and the SMM results seem to deteriorate as the particles move downstream. For this reason, we must choose a larger cell length than $6\lambda$ for the SMM. From Figure 5, the SMM with $L_{cell}=12\lambda$ appears to capture the peaks and tails of the BTCs relatively well, but $L_{cell}=24\lambda$ gives better results. For this reason, a cell length of $L_{cell}=24\lambda$ was deemed appropriate for our system. This result implies that a length scale of 24$\lambda$ is required to capture the transport dynamics in the heterogeneous system we consider. While limited to the specific level of heterogeneity investigated here, this result quantifies the characteristic scale necessary to identify the parameters of our effective transport model with respect to the one associated with the conductivity field. This scale is likely associated with connected structures (e.g., channels) emerging in the velocity field and that are relevant for the characterization of the travel time of successive jumps. This estimate is fully consistent with the detailed work of [61].

4. Downscaling Procedure to Predict Mixing

In this section we discuss the downscaling procedures developed in order to construct full concentration fields within the upscaled model. The goal here is not to accurately approximate the fully resolved concentration field, but rather reconstruct some effective concentration field that can be used to predict our desired global measures of mixing. To obtain this information, we must develop a procedure to estimate the x and y locations of each particle at any given time.

4.1. Downscaling the Streamwise Location x

As mentioned previously, at any time $t^{\ast }$ we know which cell each particle is in, how long it has been in that cell, and how much longer it will be there. Consider a particle that is in Cell $n+1$ at time $t^{\ast }$. The particle will travel through that cell over an amount of time ${\tau }^{(n+1)}$ and has been in that cell for a time equal to $t^{\ast }-t^{\left(n\right)}$, where $t^{\left(n\right)}$ is the total time traveled by the particle upon entering Cell $n+1$. In order to make an educated guess for the x location of the particle within the cell at time $t^{\ast }$, we assign each particle a mean longitudinal velocity through the cell equal to $\overline{v^{(n+1)}}=L_{cell}/{\tau }^{(n+1)}$. For our downscaling procedure in x, we assume that the particle moves straight across Cell $n+1$ with a uniform velocity equal to $\overline{v^{(n+1)}}$. Thus, we linearly interpolate along this path to determine the downscaled x position, $x^{\ast }$, i.e.,

$$x^{\ast }=(t^{\ast }-t^{\left(n\right)})\overline{v^{(n+1)}}+x^{\left(n\right)}.$$

(9)

This is the same choice for predicting x locations used in [53] and this process is illustrated schematically in Figure 6.

Figure 7 shows the histogram of particle x positions for both the fully resolved and upscaled models at various points in time throughout the simulation. This figure shows that the downscaling method described here does a reasonable job of predicting the x positions of particles relative to the small scale simulation. While it is certainly not perfect, we deem this sufficient and choose it as the method we will use to predict particle x locations for all cases. A feature that can be observed is that longitudinal particles’ locations are better matched as time advances (compare, for instance, results in Figure 7a–c to Figure 7d–f). This is consistent with the fact that our method considers a constant longitudinal velocity across a single step and therefore ignores subscale fluctuations that appear to average out across multiple steps of a single particle. Finding a good method to guess the y location of a given particles requires more effort, and we propose several different methods in the following sections.

4.2. Downscaling the Spanwise Location y

4.2.1. Method 0

Before developing more sophisticated downscaling procedures to predict the particle y locations, we start with probably the most naive approach; that is, we assume we know nothing about a particle’s y location and allow it to take any value with uniform probability between $y_{min}$ and $y_{max}$. Thus, the downscaled x and y positions are given by

$$\begin{array}{cc}\hfill x^{\ast }& =(t^{\ast }-t^{\left(n\right)})\overline{v^{(n+1)}}+x^{\left(n\right)}\hfill \\ \hfill y^{\ast }& =(y_{max}-y_{min}){\eta }^{\ast }+y_{min}\hfill \end{array}$$

(10)

where ${\eta }^{\ast }\in U(0,1)$. This method corresponds to a case where information about the probable y positions of particles is either unavailable or unimportant.

Figure 2 shows the calculation of M vs. time and Figure 3 shows E vs. time for both the fully resolved model and the upscaled model using Method 0. From Figure 2 and Figure 3, it is clear that Method 0 causes our upscaled model to strongly overpredict mixing. This result is unsurprising, as the placement of particles in the y direction based on a uniform distribution is equivalent to complete mixing of the solute in the y-direction. In reality, particles traveling on similar streamlines will tend to congregate with one another, leading to large transversal fluctuations in the concentration field. For this reason, a uniform distribution is not suitable here. More knowledge of the system is needed in order to adequately represent the particle y locations in the downscaling procedure.

4.2.2. Method 1

As observed in the previous section, it is crucial to incorporate more information from the system in the prediction of the particle positions in the y direction. Recall that in our parameterization step we recorded the joint distribution $f(y_0,y_1,{\tau }_1,{\tau }_2)$. For Method 1, we assign each particle in the upscaled model a y location equal to its initial y position, $y_0$, at all times. Therefore, the downscaled x and y positions are given by

$$\begin{array}{cc}\hfill x^{\ast }& =(t^{\ast }-t^{\left(n\right)})\overline{v^{(n+1)}}+x^{\left(n\right)}\hfill \\ \hfill y^{\ast }& =y_0.\hfill \end{array}$$

(11)

This process is illustrated in Figure 8. With this method, we are proposing that the particle’s initial position in y may be a sufficient approximation at all times in our upscaled model. From a statistical perspective this aligns with the idea that this might be its most likely location.

From Figure 2 and Figure 3, it is observed that there is a large improvement in the prediction of $M\left(t\right)$ and $E\left(t\right)$ using Method 1 relative to Method 0. This method does well at earlier times when particles have not moved very far from their initial position, but starts to under-predict mixing at ∼t = 10 as the particles have moved farther downstream and their initial y position is no longer as good of an approximation. This under-prediction of mixing can be attributed to the fact that particles would not move along a straight streamline in reality. This method does not allow for any form of transverse spreading, which is something that should be taken into account. Moving forward, we aim to improve this result by incorporating more information from the system.

4.2.3. Method 2

We know in reality that particles are not traveling along straight horizontal streamlines through our heterogeneous flow. The next step in improving our downscaling procedure would be to incorporate this idea and attempt to simulate a more realistic particle trajectory in our upscaled model. We adjust the Langevin equation to include a y component in our Spatial Markov model, i.e.,

$$\begin{array}{cc}\hfill x^{(n+1)}& =x^{\left(n\right)}+L_{cell}\hfill \\ \hfill y^{(n+1)}& =y^{\left(n\right)}+\Delta y\hfill \\ \hfill t^{(n+1)}& =t^{\left(n\right)}+{\tau }^{(n+1)}n=0,1,2,\dots .\hfill \end{array}$$

(12)

Here, $\Delta y=y_1-y_0$ comes from $f(y_0,y_1,{\tau }_1,{\tau }_2)$, which was determined during the parameterization step of the model. Each $\tau$ value has an associated $\Delta y$ value; $\tau$ is drawn as usual using the transition matrix and its corresponding $\Delta y$ value is selected for each particle at every step. Then in the downscaling procedure, we linearly interpolate between $y^{\left(n\right)}$ and $y^{(n+1)}$ over the cell length $L_{cell}$ to the point $x^{\ast }$, i.e.,

$$\begin{array}{cc}\hfill x^{\ast }& =(t^{\ast }-t^{\left(n\right)})\overline{v^{(n+1)}}+x^{\left(n\right)}\hfill \\ \hfill y^{\ast }& =\frac{y^{(n+1)}-y^{\left(n\right)}}{L_{cell}}(x^{\ast }-x^{\left(n\right)})+y^{\left(n\right)}.\hfill \end{array}$$

(13)

This method is illustrated schematically in Figure 9.

The results of Method 2 in the calculation of $M\left(t\right)$ and $E\left(t\right)$ compared with the fully resolved model results are shown on Figure 2 and Figure 3, respectively. This method appears to show an improvement over Method 1, except at early times (before time = 10) when most of the particles are still in the first cell. It is worth noting that at all times throughout our simulation this method is over-predicting mixing. This can be attributed to the fact that Method 2 is not capable of incorporating the idea that particles traveling near each other are likely to have similar trajectories. In this version of the upscaled model, each particle draws a new ($\tau$,$\Delta y$) pair when it enters a new cell. The method is unable to account for the fact that particles traveling near each other will likely have similar $\tau$ and $\Delta y$ values. While $\tau$ values are correlated by the implementation of the transition matrix in the SMM, particles with similar travel times can have very different values of $\Delta y$ associated with them due to the highly heterogeneous nature of the flow field. This results in the observed over-mixing for Method 2 in Figure 2 and Figure 3.

4.2.4. Method 3

In this final method, we will add some complexity to the parameterization step of the model. We previously measured the joint distribution $f(y_0,y_1,{\tau }_1,{\tau }_2)$ in the model parameterization. Now, we expand this step by recording information on the particle trajectories within the first cell. To measure these trajectories, we record the y positions for each particle as they travel across the first cell at $N_t$ evenly spaced points in the x-direction. Therefore, in our updated parameterization step we record the joint distribution $f(y_0,y_{t,1},y_{t,2},\dots ,y_{t,N_t}=y_1,{\tau }_1,{\tau }_2)$, where $y_{t,\#}$ are the particle y positions measured at $N_t$ points along their trajectory. $N_t=20$ was determined to be a sufficient number of sampling positions for capturing the particle trajectories within the cell.

In order to incorporate this additional particle trajectory information into our downscaling procedure, we separate the particles into zones based on their initial y position, $y_0$. These zones were defined by separating the $y_0$ particle positions into a number $N_z$ of equiprobable regions. Since our solute line injection defined by the particle $y_0$ positions is flux-weighted, these equiprobable zones will not be uniformly spaced in the y direction. $N_z=20$ was selected here, so that each zone will contain trajectory information from $\raisebox{1ex}{$N_p$}\!\left/ \!\raisebox{-1ex}{$N_z$}\right.=5000$ particles. A broad range of zone numbers were tested and the results were not significantly affected, indicating the robustness of this method. The histogram of $y_0$ values is shown in Figure 10a, where each color represents a different zone number.

After defining our zones, we then create a set of y values, $Y_j$, for each zone number j. The set $Y_j$ contains all of the y positions measured in the trajectories for every particle that had $y_0$ in Zone j. This means that $Y_j$ contains every y value from each particle trajectory defined by $y_0,y_{t,1},y_{t,2},\dots ,y_{t,N_t}=y_1$ that had its initial y position, $y_0$, located in zone number j. In this study, $Y_j$ is a set containing $N_t\times \raisebox{1ex}{$N_p$}\!\left/ \!\raisebox{-1ex}{$N_z$}\right.=20\times 5000$y values. Figure 10b shows the histogram of the set of y values $Y_j$ for each zone number j. The adjustment of the parameterization step is illustrated in Figure 11.

After defining our zones and determining the set of y values $Y_j$ associated with each zone number j, we will now use this information in our downscaling procedure. In this method, we will determine the downscaled y position of each particle at some time $t^{\ast }$ by randomly selecting a y value from the set of possible y values, $Y_j$, associated with each particle’s zone number j. The particle’s zone number depends exclusively on its initial y position, $y_0$, and does not change at any point throughout the simulation. For example, a particle that had $y_0$ in Zone 2 will always draw from the set of possible y values associated with that zone, $Y_2$. Thus, the downscaled x and y positions at time $t^{\ast }$ are given by

$$\begin{array}{cc}\hfill x^{\ast }& =(t^{\ast }-t^{\left(n\right)})\overline{v^{(n+1)}}+x^{\left(n\right)}\hfill \\ \hfill y^{\ast }& =y_j^{(n+1)},\hfill \end{array}$$

(14)

where $y_j^{(n+1)}$ is a y value selected randomly with uniform probability at each Spatial Markov step from the set $Y_j$ associated with the particle’s zone number j.

This method incorporates the idea that particles will in general travel along paths that have y values similar to their initial y position. They will not likely be making large jumps in the y direction, but they will also not simply maintain their initial y position as we assumed in Method 1. By selecting from a set of y values based on each particle’s initial y position, we are assuming that particles over time will have the same distribution in y as they did in the first cell during the model parameterization. The y positions of particles are assigned through a regionalization of the transversal dimension y into zones, preventing them from making excessively large jumps in the spanwise direction. The results of Method 3 in the measurement of $M\left(t\right)$ and $E\left(t\right)$ compared to the fully resolved model and the other upscaled models are shown on Figure 2 and Figure 3, respectively. As observed in these figures, Method 3 performs better than the other upscaled models in estimating $M\left(t\right)$ and $E\left(t\right)$.

5. Discussion

Now, we compare the different upscaled methods both qualitatively and quantitatively. From Figure 2 and Figure 3, it is clear that Methods 1, 2, and 3 all perform very well relative to Method 0 in predicting $M\left(t\right)$ and $E\left(t\right)$ as compared to the fully resolved random walk algorithm. Here, we will discuss the successes of each of these methods as well as their shortcomings.

Figure 12 shows the particle locations at three different times throughout the simulations for the fully resolved random walk model and the upscaled methods 0, 1, 2, and 3. As already noted above, our goal was not necessarily to reproduce the concentration fields of the fully resolved simulation, but to capture the key features that enable an accurate estimate of $M\left(t\right)$ and $E\left(t\right)$, such as the peaks in concentration and the spreading of the plume. Many closures for predicting reactive transport involve integrated nonlinear terms like these e.g., [62,63] and so accurately estimating them is an essential first step to accurate upscaling of reactive transport. As observed in Figure 12, Method 0 clearly over-predicts mixing at all times. This is consistent with the results of Figure 2 and Figure 3. In Figure 12, Method 1 qualitatively appears similar to our fully resolved simulations at $t=1$ and 10, but it appears that the concentration peaks may be too high and there is not enough spreading of the plume at $t=100$. This is in agreement with the under-prediction of mixing observed for Method 1 after $t=10$ in Figure 2 and Figure 3. Methods 2 and 3 also appear qualitatively similar to the fully resolved simulation at $t=1$ and 10, but may have slightly too much spreading. It is also observed in Figure 12 that there is some discontinuous structure for Method 3 at $t=10$. This corresponds to the point in time where the particles have begun to enter the second SMM cell. At each step in the SMM, each particle selects a new y value from the set $Y_j$ associated with that particle’s zone number j. This means that as a particle enters a new SMM cell it selects a new y value that is somewhere between the minimum and maximum values of $Y_j$, resulting in the observed structure. At $t=100$, Method 2 clearly has too much spreading, but may still be capturing the peaks in concentration. Figure 2 and Figure 3 are in agreement with these observations as Method 2 slightly over-predicts mixing at all times after $t=0.03$. At $t=100$, the plume of Method 3 has a different shape than the fully resolved simulation, but the amount of spreading and regions of high concentration are comparable. Figure 2 and Figure 3 show that Method 3 also very slightly over-predicts mixing at all times except very early times, but in general show results that are very close to that of the fully resolved simulation. Upon visual inspection of Figure 2 and Figure 3, Method 3 appears to give the best prediction of $M\left(t\right)$ and $E\left(t\right)$.

After qualitatively examining each of these methods, we now want to quantify which method performs the best. To do this, we use the mean absolute error on a log scale, defined as

$$ϵ\left(\beta \right)=\frac{{\sum }_{i=1}^{N_{pts}}|{\log }_{10}{\beta }_{1,i}-{\log }_{10}{\beta }_{2,i}|}{N_{pts}}$$

(15)

where $\beta$ is either $M\left(t\right)$ or $E\left(t\right)$, ${\beta }_1$ is the mixing metric calculated in the fully resolved model, ${\beta }_2$ is the mixing metric calculated in the upscaled model, and $N_{pts}$ is the number of points in time at which $M\left(t\right)$ and $E\left(t\right)$ are measured. We choose to calculate the mean absolute error on a log scale so as not to penalize the smaller values of $M\left(t\right)$ and $E\left(t\right)$. The results of this error calculation are shown in

Table 1. The mean absolute error of ${\log }_{10}\beta$, where $\beta$ is either $M\left(t\right)$ or $E\left(t\right)$, for each of the upscaled methods as defined by Equation (15).

	ϵ(M)	ϵ(E)
Method 0	0.5895	0.4888
Method 1	0.1275	0.0893
Method 2	0.1121	0.1217
Method 3	0.0352	0.0646

. As expected, Method 0 has a significantly larger error than Methods 1, 2, and 3. Methods 1 and 2 had very similar error values, and Method 3 had the smallest error in the measurement of both $M\left(t\right)$ and $E\left(t\right)$. Based on this error metric, Method 3 is the option that performs the best. This is consistent with the qualitative examination as well.

6. Conclusions

In this work, we have extended the Spatial Markov model to predict effective mixing in flows through idealized two-dimensional heterogeneous porous media. To do this, we developed several downscaling methods to predict the particle locations within the upscaled model and thereby approximating solute distribution within the domain. From these predicted particle locations, we were able to generate concentration fields that could then be used to calculate the evolution of the solute mixing, as quantified by integral of the squared concentration, $M\left(t\right)$, and the dilution index, $E\left(t\right)$, in time. The results of each of these downscaling methods developed for the SMM were compared to a fully resolved random walk simulation. In order to be able to predict $M\left(t\right)$ and $E\left(t\right)$ accurately with our upscaled model, it was crucial to capture the same amount of variation in the concentration field as in the fully resolved model. Each of these upscaled methods require nearly the same computation time and are on the order of 1000 times faster to run than the fully resolved method. Although Method 3 does add more complexity to the upscaled model than the other methods presented here, it is as computationally efficient as the other methods and provides the best prediction of $E\left(t\right)$ and $M\left(t\right)$. For this reason, we recommend Method 3 as an extension of the Spatial Markov model to predict mixing in uniform mean flow simulations in heterogeneous porous media.

Although we were able to predict $M\left(t\right)$ and $E\left(t\right)$ within an acceptable margin of error using our upscaled model, there is still room for improvement. Our results show that a key feature to be encoded in the upscaled approaches is that particles traveling near each other will have similar trajectories. Method 3 provides good results in terms of matching average mixing behavior through a rough regionalization of transversal particle positions. Still, we observe a slight over-prediction of mixing and the methodology is not able to capture the topological structure of the solute plume. We envision that the method can be further improved in future works, particularly with the aim of generalizing the approach to broader transport settings beyond one-dimensional transport.