Dispersive transport and symmetry of the dispersion tensor in porous media

The macroscopic laws controlling the advection and diffusion of solute at the scale of the porous continuum are derived in a general manner that does not place limitations on the geometry and time evolution of the pore space. Special focus is given to the definition and symmetry of the dispersion tensor that is controlling how a solute plume spreads out. We show that the dispersion tensor is not symmetric and that the asymmetry derives from the advective derivative in the pore-scale advection-diffusion equation. When flow is spatially variable across a voxel, such as in the presence of a permeability gradient, the amount of asymmetry can be large. As first shown by Auriault [J.-L. Auriault et al. Transp. Porous Med. 85, 771 (2010)] in the limit of low Péclet number, we show that at any Péclet number, the dispersion tensor $D_{ij}$ satisfies the flow-reversal symmetry $D_{ij}(+\mathbf{q})=D_{ji}(-\mathbf{q})$ where $\mathbf{q}$ is the mean flow in the voxel under analysis; however, Reynold's number must be sufficiently small that the flow is reversible when the force driving the flow changes sign. We also demonstrate these symmetries using lattice-Boltzmann simulations and discuss some subtle aspects of how to measure the dispersion tensor numerically. In particular, the numerical experiments demonstrate that the off-diagonal components of the dispersion tensor are antisymmetric which is consistent with the analytical dependence on the average flow gradients that we propose for these off-diagonal components.

Figure 1

A cubic voxel of porous material $\mathrm{\Omega }\left(\mathbf{r}\right)$ centered on point $\mathbf{r}=(r_1,r_2,r_3)$ within the larger porous system and having local coordinates $\mathbf{x}=(x_1,x_2,x_3)$ as shown.

Figure 2

The D2Q6 triangular lattice focusing on a particular node ${\mathbf{x}}_b$ sitting on the boundary surrounded by the six links $i=1,6$ that have total mass populations $N_i$ and solute mass concentrations ${\mathrm{\Delta }}_i$. The constant velocities ${\mathbf{c}}_i$ are shown as arrows and are given by ${\mathbf{c}}_1={\widehat{\mathbf{x}}}_1, {\mathbf{c}}_2={\widehat{\mathbf{x}}}_1/2+\sqrt{3}{\widehat{\mathbf{x}}}_2/2, {\mathbf{c}}_3=-{\widehat{\mathbf{x}}}_1/2+\sqrt{3}{\widehat{\mathbf{x}}}_2/2, {\mathbf{c}}_4=-{\widehat{\mathbf{x}}}_1, {\mathbf{c}}_5=-{\widehat{\mathbf{x}}}_1/2-\sqrt{3}{\widehat{\mathbf{x}}}_2/2, {\mathbf{c}}_6={\widehat{\mathbf{x}}}_1/2-\sqrt{3}{\widehat{\mathbf{x}}}_2/2$.

Figure 3

Sinusoidal forcing $\mathbf{f}={\widehat{\mathbf{x}}}_{2}sin(2\pi x_1/L)$ with $\text{Pe}=3.0$. (a) The flow field. (b) The concentration distribution near steady state when the applied concentration gradient is in the vertical direction. (c) When the concentration gradient is in the horizontal direction. Yellow denotes zero concentration, white positive, and blue negative concentration.

Figure 4

The dispersion components corresponding to the geometry and forcing in Fig. 3 but with $\text{Pe}=q_{\max }{L}_1/D_m=0.75$, 1.5, and 3.0 from bottom to top. The diffusion time constant is defined $t_{\mathrm{Diff}}=L_1^2/D_m$. The lower curve near the origin is $D_{21}(-f)$ for each Pe.

Figure 5

(a) The flow field corresponding to Eq. (106) when no solid obstacles are present. (b) The concentration field near steady state when there is an imposed concentration gradient in the vertical direction (in the same direction as the flow field). (c) When the concentration gradient is in the horizontal direction.

Figure 6

The dispersion coefficients corresponding to the geometry and forcing in Fig. 5 with $\text{Pe}=\mathrm{\Delta }{q}_2{L}_1/D_m$ ranging from 1 to 5. The wiggles in the reversed-flow ($-f$) results is a numerical artifact.

Figure 7

Upper panel: The off-diagonal components $D_{12}\left(f\right)$ and $D_{21}\left(f\right)$ as a function of Pe. The data fits are made using $D_{21}=\mathrm{Pe}\left(a+bln\mathrm{Pe}\right)$ with $a$ and $b$ given by Eq. (107). The numerical results demonstrate the anti-symmetric behavior $D_{12}\left(f\right)\approx -D_{21}\left(f\right)$. Lower panel: The difference $[D_{21}\left(f\right)-D_{12}\left(f\right)]/D_m$ [which is $\approx 2D_{21}\left(f\right)/D_m$]. Both panels correspond to the steady-state behavior observed in Fig. 6.

Figure 8

The $\alpha$(Pe) coefficient corresponding to Fig. 6 and obtained from Eq. (111).