Effects of incomplete mixing on reactive transport in flows through heterogeneous porous media

The phenomenon of incomplete mixing reduces bulk effective reaction rates in reactive transport. Many existing models do not account for these effects, resulting in the overestimation of reaction rates in laboratory and field settings. To date, most studies on incomplete mixing have focused on diffusive systems; here, we extend these to explore the role that flow heterogeneity has on incomplete mixing. To do this, we examine reactive transport using a Lagrangian reactive particle tracking algorithm in two-dimensional idealized heterogeneous porous media. Contingent on the nondimensional Peclet and Damköhler numbers in the system, it was found that near well-mixed behavior could be observed at late times in the heterogeneous flow field simulations. We look at three common flow deformation metrics that describe the enhancement of mixing in the flow due to velocity gradients: the Okubo-Weiss parameter $\left(\theta \right)$, the largest eigenvalue of the Cauchy-Green strain tensor $\left({\lambda }_C\right)$, and the finite-time Lyapunov exponent $\left(\mathrm{\Lambda }\right)$. Strong mixing regions in the heterogeneous flow field identified by these metrics were found to correspond to regions with higher numbers of reactions, but the infrequency of these regions compared to the large numbers of reactions occurring elsewhere in the domain imply that these strong mixing regions are insufficient in explaining the observed near well-mixed behavior. Since it was found that reactive transport in these heterogeneous flows could overcome the effects of incomplete mixing, we also search for a closure for the mean concentration. The conservative quantity $\overline{u^2}$, where $u=C_A-C_B$, was found to predict the late time scaling of the mean concentration, i.e., $\overline{C_i}\sim \overline{u^2}$.

Figure 1

(Top) The natural logarithm of the permeability field, $ln\left(\kappa \right)$. (Middle) Horizontal velocity $v_x$ and (Bottom) vertical velocity $v_y$. Lengths have been normalized to the domain extent in the $y$ direction.

Figure 2

The black line corresponds to the well-mixed solution for the mean concentration and scales like $t^{-1}$ at late times. The red line is the result for the mean concentration in a purely diffusive system and scales as $t^{-1/2}$ at late times. The blue line corresponds to the mean concentration from our simulation in an idealized heterogeneous porous medium with Pe = 1200 and Da = 14 400. At late times, our model results for this particular heterogeneous field scales as if it were almost well-mixed.

Figure 3

The A and B particle locations for a system with Pe = 1200 and Da = 14 400 are shown at an early time $\left(t^{*}=12\right)$, early intermediate time $\left(t^{*}=54.1\right)$, late intermediate time $\left(t^{*}=155.4\right)$, and late time $\left(t^{*}=630.6\right)$. It is observed that the islands grow in size over time, demonstrating the presence of incomplete mixing in the system.

Figure 4

The mean concentration vs. time is shown for the simulation results of nine cases with different Peclet and Damköhler numbers. The red lines are the well-mixed solution defined in Eq. (7) and the blue lines are the simulation results for heterogeneous flow field cases at the specified Peclet and Damköhler numbers.

Figure 5

(a) The two-dimensional histogram of reaction locations for 100 realizations of the Pe = 1200 and Da = 14 400 case. Every time a reaction occurs, the location of the reaction is recorded. (b) The corresponding fields of the Okubo-Weiss parameter $\left(\theta \right)$, (c) ${\lambda }_C$, and (d) the finite-time Lyapunov exponent $\left(\mathrm{\Lambda }\right)$ are shown.

Figure 6

A scatterplot of the local values of $\theta ,{\lambda }_C$, and $\mathrm{\Lambda }$ and the number of reactions at the same locations for the Pe = 1200 and Da = 14 400 case.

Figure 7

The average $\theta ,{\lambda }_C$, and $\mathrm{\Lambda }$ values at locations where the amount of reactions that have occurred is above the threshold $\tau$ indicated on the x-axis for the Pe = 1200 and Da = 14 400 case.

Figure 8

The evolution of $\overline{C_A^{\prime }{C}_A^{\prime }},\overline{C_A^{\prime }{C}_B^{\prime }}$, and $\overline{u^2}$ in time for a system with Pe = 1200 and Da = 14 400. At late times, $\overline{u^2}$ scales like $\overline{C_A^{\prime }{C}_A^{\prime }}$. The quantities $\overline{C_A^{\prime }{C}_A^{\prime }},\overline{C_A^{\prime }{C}_B^{\prime }}$, and $\overline{u^2}$ have been nondimensionalized by $C_0^2$.

Figure 9

The evolution of $\sqrt{-\overline{C_A^{\prime }{C}_B^{\prime }}\phantom{{}^{l^l}}}$, and $\overline{C_A^{\prime }{C}_A^{\prime }}$ in time for a system with Pe = 1200 and Da = 14 400. At late times, $\overline{C_A^{\prime }{C}_A^{\prime }}$ scales like $\sqrt{-\overline{C_A^{\prime }{C}_B^{\prime }}\phantom{{}^{l^1}}}$.

Figure 10

The evolution of $\overline{C_i}$ and $\overline{u^2}$ in time. $\overline{C_i}$ is given by the blue line and the black line is $\overline{u^2}$. At late times, $\overline{C_i}\sim \overline{u^2}$.