Simulating anomalous dispersion in porous media using the unstructured lattice Boltzmann method

Flow in porous media is a significant challenge to many computational fluid dynamics methods because of the complex boundaries separating pore fluid and host medium. However, the rapid development of the lattice Boltzmann methods and experimental imaging techniques now allow us to efficiently and robustly simulate flows in the pore space of porous rocks. Here we study the flow and dispersion in the pore space of limestone samples using the unstructured, characteristic based off-lattice Boltzmann method. We use the method to investigate the anomalous dispersion of particles in the pore space. We further show that the complex pore network limits the effectivity by which pollutants in the pore space can be removed by continuous flushing. In the smallest pores, diffusive transport dominates over advective transport and therefore cycles of flushing and no flushing, respectively, might be a more efficient strategy for pollutant removal.

1. Introduction

Recent numerical simulations of flow in porous media in both two and three dimensions have revealed highly non-trivial flow patterns. Particles passively advected by the flow in the pore-space, so-called Lagrangian particles, follow trajectories with periods of almost stagnation punctuated by bursts of fast fluctuating motion [1–3]. In a Berea sandstone sample, it was shown by a three dimensional simulation that this intermittent behavior was equally significant in both longitudinal (the mean flow direction) and transverse directions [3]. The chaotic Lagrangian dynamic emerges despite the fact that the inertia of the fluid is negligible in many natural and industrial settings, i.e., the flow happens at low Reynolds numbers. Nonetheless, the flow might exhibit a complicated spatio-temporal behavior owing to the high heterogeneity of the pore space. The heterogeneity in fact gives rise to anomalous (non-Fickian) transport properties similar to those typically associated to high Reynolds number flow [4]. The anomalous transport has implications for the local dispersion of particles in a porous medium. It can alter the reaction rates of transported reactants or accelerate the small scale dispersion of pollutants.

In addition to the strong intermittent behavior, it has been found in both two and three dimensional simulations that the velocity autocorrelation, as a function of space, is short ranged while as a function of time the longitudinal velocity exhibits long time correlations. In the three dimensions, it was found [3] that the transverse velocity autocorrelation decays faster, though the absolute value is strongly correlated in time, a feature also observed in turbulent flows [5]. These remarkable properties have led to the suggestion that the Lagrangian particle velocities is given by a Markov process in space (at equidistant positions along the Lagrangian trajectories) and not in time [6, 7]. In a Markov process the particle motion can be seen as a correlated continuous time random walk (CTRW). Using a correlated CTRW, several signatures of anomalous transport behavior were accurately reproduced such as the long tails of the first passage time distribution [6], the non-linear scaling of the second centered moment of the particles longitudinal and transverse displacements [3, 6, 8], the probability distribution function of the Lagrangian velocity increments for different time lags among others [1].

Interesting results were also obtained in the study of scalar mixing in two dimensional porous media with different structural disorder [9]. Numerical simulations showed that the scalar concentration forms, at advanced time, a fractal-like structure composed of stripes of high concentration and with high lacunarity. The probability distribution function of the pair separation of advecting particles was fitted to a log-normal distribution and the averaged squared separation was found to grow with time as a power law whose exponent depends on the geometry. For weak and strong heterogeneities of the porous sample the exponents were found to be 1.8 and 2.8, respectively. Remarkably, we can see that the second exponent is very close to the Richardson dispersion law, characteristic of the separation of particles whose distance lie in the inertial range interval in turbulent flows. The authors proposed a stochastic model based on ad hoc arguments to describe the behavior of the pair separation of advecting particles l(t), which reflects the multiplicative nature of the stretching processes and accounts for the observed short range temporal correlations of the Lagrangian stretching rates.

Here we formulate a version of the lattice Boltzmann method on unstructured grids to efficiently simulate flow in complex pore geometries. We use the method to assess the Lagrangian dynamics of flows in the pore space of limestone samples computed by x-ray tomography. In order to faithfully compute the broadly distributed fluid velocities, it is required to represent the pore boundaries with a high resolution. Compared to classical lattice Boltzmann methods on regular grids, the unstructured grids offer a superior resolution and adaptivity [10]. We further analyze the impact of the complex pore geometry on pollutant removal from the pore space.

2. Materials and Methods

2.1. Porous Rock Samples

We have performed our pollutant dispersion studies on two different, natural porous samples, discretized using tetrahedral elements (the properties of the tetrahedral meshes are presented in Table 1).

TABLE 1

Table 1. Properties of the tetrahedral meshes used in the simulations and the wall-clock time spent per 5000 iterations of the unstructured lattice Boltzmann simulation (Intel^® Xeon^® 2.6 GHz).

Sample A (Figure 1) is chalk from off-shore drilling in the North Sea region (Hod chalk #15). Sample B (Figure 2) came from an outcrop of bryozoan chalk in Rødvig, Denmark. The digital 3D images were obtained by x-ray nanotomography [11] measured at beamline ID22 at the European Synchrotron Radiation Facility, France. The reconstructed volume of Sample A had a voxel size of 25 nm and an optical resolution about 150 nm. Details about the data collection and reconstruction can be found in Müter et al. [12]. The reconstructed images were corrected for ring artifacts [13] before segmentation by a dual filtering and Otsu thresholding procedure [14]. For the LBM calculations we used a subvolume of 100³ voxels, which gives a side length of 2.5 microns. Sample B was reconstructed and segmented in the same way as Sample A. The only difference is that it was reconstructed with a voxel size of 50 nm leading to an optical resolution of about 300 nm. The LBM calculations were performed on a subvolume of 400³ voxels (side length of 20 microns). In Sample A, the diameter of the pore throats corresponds to at least 6 voxel sizes, while the typical pore diameter is 30–40 voxels. In Sample B, the smallest pore throats are about 5 voxels, however, the flow is dominated by a large, well connected pore system, where the pore diameter is 50 voxels and above.

FIGURE 1

Figure 1. Left: The surface plot of the investigated volume of Sample A, colored according to the surface velocities. The values on the solid boundary are at least three orders of magnitude lower than the values on the outlet (top surface of the mesh). Right: The streamlines of the steady-state velocity field, generated for τ = 0.02, Δρ = 2.5 · 10⁻⁶. Both images were generated using ParaView 4.3.1. [32].

FIGURE 2

Figure 2. Left: The surface plot of the investigated volume of Sample B, colored according to the surface velocities. The values on the solid boundary are at least three orders of magnitude lower than the values on the outlet (top surface of the mesh). Right: The streamlines of the steady-state velocity field, generated for τ = 0.1, Δρ = 2.5 · 10⁻⁶. Both images were generated using ParaView 4.3.1. [32].

In order to construct the tetrahedral meshes, the initial triangular surface meshes were produced using the marching cubes isosurface as implemented in ImageJ [15]. The resolution of those surface meshes is the same as in the corresponding tomographic data (1 voxel accounting for 1-2 triangles). A tetrahedral volume mesh was constructed from the surface mesh using Tetgen [16].

2.2. Unstructured Lattice Boltzmann Method

The appeal of the lattice Boltzmann method to the geophysics research lies in its ability to model and simulate several types of phenomena related to reactive fluid flow in porous rocks: single and multiphase flows [17], unsteady flows [18], flows in complex geometries [19, 20] at a wide range of Reynolds numbers, as well as chemical processes, such as dissolution and precipitation [21]. For our experiments, we have used a variant of the finite element, off-lattice Boltzmann method [22], which uses the characteristic based integration, first proposed by Zienkiewicz and Codina [23]. We chose this method over the finite volume approach [10, 24] due to its improved stability. The characteristic based, off-lattice Boltzmann method has been specifically designed to remove the coupling between the positions of the computational grid's nodes and the discrete set of velocities, while retaining the stability of the standard, regular grid-based approach. The new scheme is derived directly from the spatio-temporal continuous formulation of the lattice Boltzmann equation

$\begin{array}{cc}\frac{\partial f_i}{\partial t}+c_i·\nabla f_i={\Omega }_i(f(x,t)),& (1)\end{array}$

where f_i(x, t) is the particle velocity distribution along the i-th discrete direction, c_i is the associated discrete velocity, and Ω_i is the collision operator, accounting for the rate of change of f_i due to collision. We have employed the two relaxation time collision (TRT) operator [19, 25], which is a special case of a multi relaxation time collision operator (MRT, [26]), most commonly used in the lattice Boltzmann studies of porous flow due to its improved treatment of the solid boundary conditions in comparison to the standard, single relaxation time Bhatnagar-Gross-Krook (BGK) operator [19].

2.2.1. Integration Along Characteristics

Integrating Equation (1) along the characteristic curves and approximating the right-hand side yields

$\begin{array}{cc}\begin{array}{l}{f}_i^{n+1}(x+c_i\delta t)=f_i^n(x)+\delta t(\theta {\Omega }_i^{n+1}(f^{n+1}(x+c_i\delta t))\hfill \\ +(1-\theta ){\Omega }_i^n(f^n(x))),\hfill \end{array}\hfill & (2)\hfill \end{array}$

where n denotes the time step (iteration number), 0 ≤ θ ≤ 1 refers to the choice of the time integration method for the collision operator (0 for fully explicit, 1 for fully implicit, 0.5 for Crank-Nicolson), fⁿ_i(x) = f_i(x, tⁿ), and f^{n + 1}_i(x) = f_i(x, tⁿ+ δt). In the general case of a MRT collision operator Ω_i(f(x, t)) = ∑_jA_ij(f^eq_j − f_j) [26] the implicit character of Equation (2) can be removed by introducing new variable

$\begin{array}{cc}{g}_i(x,t)=f_i(x,t)-\delta t\theta {\Omega }_i(f(x,t)),& (3)\end{array}$

and expressing the collision operator in terms of g [22]

$\begin{array}{cc}\Omega (f(x,t))={(I+\delta t\theta A)}^{-1}\Omega (g(x,t)).& (4)\end{array}$

As the collision operator preserves both mass and momentum, the moments of g_i equal corresponding moments of f_i, ∑_ig_i = ρ, ∑_ig_ic_i = ρu, hence also g^eq_i = f^eq_i. Finally, Equation (2) becomes

$\begin{array}{cc}{g}_i^{n+1}(x+c_i\delta t)=g_i^n(x)+\delta t{\displaystyle \sum _j{\widehat{A}}_{ij}}(g_j^{eq,n}(x)-g_j^n(x)),& (5)\end{array}$

where the modified collision matrix reads $\widehat{A}$ = (I + δtθA)⁻¹A. By substituting A by A_BGK = τ⁻¹I, which corresponds to the single relaxation time BGK collision operator, one readily obtains $\widehat{A}$_BGK = (τ + θδt)⁻¹I. The resulting shift in the effective relaxation time $\widehat{\tau }$ = τ + θδt accounts for the unconditional stability of the Crank-Nicolson method and the implicit Euler method, since the condition δt < 2$\widehat{\tau }$ is always fulfilled as long as θ ≥ $\frac{1}{2}$ and τ > 0. The standard Courant-Friedrichs-Lewy restriction on δt still applies.

2.2.2. Finite Element Scheme

In order to solve Equation (5) numerically on an unstructured grid, one has to be able to express the value of gⁿ⁺¹_i at a given point x* in terms of the value of gⁿ_i (and its derivatives) evaluated one time step earlier, at the same point. In order to enhance stability, this is typically done by the second order Taylor expansion of gⁿ_i(x* − c_iδt) [22, 23]

$\begin{array}{cc}\begin{array}{l}{g}_i^n(x^{\ast }-c_i\delta t)=g_i^n(x^{\ast })-\delta t{c}_{ir}{\partial }_r{g}_i^n(x^{\ast })\hfill \\ +\frac{\delta t^2}{2}{c}_{ir}{c}_{is}{\partial }_s{\partial }_r{g}_i^n(x^{\ast })+{𝒪}(\delta t^3).\hfill \end{array}\hfill & (6)\hfill \end{array}$

After applying the same approximation to the equilibrium distribution function, Equation (5) takes the following form

$\begin{array}{cc}\begin{array}{l}{g}_i^{n+1}=g_i^n-\delta t\left(c_{ir}{\partial }_r{g}_i^n-{\displaystyle \sum _j{\widehat{A}}_{ij}}(g_i^{eq,n}-g_i^n)\right)\hfill \\ +\frac{\delta t^2}{2}{c}_{ir}{\partial }_r\left(c_{is}{\partial }_s{g}_i^n-2{\displaystyle \sum _j{\widehat{A}}_{ij}}(g_i^{eq,n}-g_i^n)\right)\hfill \\ +\frac{\delta t^3}{2}{c}_{ir}{c}_{is}{\partial }_s{\partial }_r\left({\displaystyle \sum _j{\widehat{A}}_{ij}}(g_j^{eq,n}-g_j^n)\right),\hfill \end{array}\hfill & (7)\hfill \end{array}$

which is now suitable for discretization in space using the Galerkin finite element method. Here, the spatial decomposition using linear, tetrahedral elements has been applied. The modified particle distribution functions are sampled at the nodes (vertices) of the tetrahedral mesh, and interpolated at other points

$\begin{array}{cc}{\tilde{g}}_i^n(x)=\varphi {(x)}^T{\tilde{g}}_i^n& (8)\end{array}$

where ${\tilde{g}}_i^n$ is the approximate solution, ${\tilde{g}}_i^n=[{\tilde{g}}_i^n(v_1),{{\tilde{g}}_i^n(v_2),\dots ,{\tilde{g}}_i^n(v_{N_V})]}^T$ is the vector storing the values of ${\tilde{g}}_i^n$ at all vertices in the mesh, whose positions are v₁, v₂, …, v_NV. Furthermore, ϕ(x)^T = [ϕ¹(x), ϕ²(x), …, ϕ^N_V(x)], where ϕ^k(x) is a piecewise-linear shape function (also known as a hat or tent function) associated with the vertex k, i.e., ϕ^k(v_j) = δ_jk, and ϕ^k is linear over each tetrahedral element. By applying the Bubnov-Galerkin method, we finally obtain the discrete, weak form of Equation (7)

$\begin{array}{cc}\begin{array}{l}M({\tilde{g}}_i^{n+1}-{\tilde{g}}_i^n)=\left(-\delta t{C}_i-\delta t^2{D}_i\right){\tilde{g}}_i^n+\left(\delta tM-\delta t^2{C}_i-\delta t^3{D}_i\right)\\ {\displaystyle \sum _j{\widehat{A}}_{ij}}({\tilde{g}}_j^{eq,n}-{\tilde{g}}_j^n),\end{array}& (9)\end{array}$

where matrices M, C_i, D_i ∈ ℝ^{N_V × N_V} are defined as

$\begin{array}{cc}\begin{array}{l}M={\displaystyle {\int }_{{𝒟}}\varphi {\varphi }^{T}dV},C_i={\displaystyle {\int }_{{𝒟}}\varphi c_{ir}{\partial }_r{\varphi }^{T}dV},\hfill \\ D_i=\frac{1}{2}{\displaystyle {\int }_{{𝒟}}{\partial }_s}\varphi c_{is}{c}_{ir}{\partial }_r{\varphi }^{T}dV.\hfill \end{array}\hfill & (10)\hfill \end{array}$

2.2.3. Boundary Conditions

The quality of any lattice Boltzmann simulation hinges on the choice of the boundary conditions used both on the solid and open boundaries. In this work, we have used a recent formulation of the solid boundary conditions for high Reynolds number flows, first proposed by Chikatamarla and Karlin [27]. Similarly to the popular bounce back method, the first step is to determine the subset D of the discrete directions that are corresponding to the “missing” populations (i.e., point away from the solid phase and into the fluid phase). The procedure for finding such subset in an unstructured grid setting is described in detail in Misztal et al. [10]. By setting the target values of the density ρ_tg and velocity u_tg at the boundary node, the local populations are recovered using the following rules

$\begin{array}{cc}{f}_i\leftarrow \{\begin{array}{l}2f_i^{eq}({\rho }_{tg},u_{tg})-f_i^{eq}(\rho ,u),i\in \overline{D}\\ f_i+f_i^{eq}({\rho }_{tg},u_{tg})-f_i^{eq}(\rho ,u),i\cancel{\in }\overline{D}\end{array},& (11)\end{array}$

where ρ = ∑_i∈Df^eq_i(ρ_tg, u_tg) + ∑_i∉Df_i, and ρu = ∑_i∈Dc_if^eq_i(ρ_tg, u_tg) + ∑_i∉Dc_if_i. In our simulations we assumed no-slip boundary condition on the static, solid walls, corresponding to u_tg = 0. The value of ρ_tg can be estimated either by using the bounce back method or, if the former approach fails (depending on the local geometry of the solid boundary), by extrapolating the values of density from the neighboring bulk nodes. We have observed that this choice of the solid boundary condition formulation yields significantly more stable simulations of flows with lower values of the kinetic viscosity and higher Reynolds number, compared to the bounce back method.

The boundary conditions at the inlet and at the outlet can also be implemented using the same procedure. We have chosen to use pressure boundary conditions, specified by the constant pressure values: p_in at the inlet, and p_out at the outlet; p_in > p_out. Since p = c²_sρ in the lattice Boltzmann method, this is equivalent to setting ρ_tg = p_in/c²_s at the inlet and ρ_tg = p_out/c²_s at the outlet (like most authors, we use c²_s = 1/3). The values of u_tg at the open boundary nodes can be approximated using the averaged values of the velocity at the neighboring bulk nodes.

3. Results

The unstructured lattice Boltzmann method allows us to robustly compute single phase flow fields in arbitrary, complex channel networks for a wide range of flow parameters. In this section, we apply this method to study the flow in each of the natural porous samples described in Section 2.1. We then use the precomputed flow fields to study pollutant dispersion in the samples.

3.1. Steady-State Flow Computation

Our first goal is to compute the steady state flow in each porous rock sample due to applying a constant pressure difference between the inlet and the outlet planes. The particle distribution functions are initialized with the equilibrium distributions corresponding to ρ = 1, u = 0 in the entire mesh. For a given time step δt, relaxation time τ, and constant inlet and outlet densities ρ_in = 1 + Δρ and ρ_out = 1, we run the lattice Boltzmann simulation until the steady state is reached. In order to detect when steady-state has been reached, we compute the inlet flow rate Q_in by evaluating the following integral numerically at every time step tⁿ

$\begin{array}{cc}{Q}_{in}={\displaystyle {\int }_{A_{in}}u·dA}\approx {\displaystyle \sum _{T_k=\langle pqr\rangle }\frac{u_p+u_q+u_r}{3}}·n_k{A}_k,& (12)\end{array}$

where the sum on the right-hand side is over each inlet triangle T_k; u_p, u_q, u_r are its nodal velocities, n_k is the unit-length normal vector to T_k, and A_k is its area. The outlet flow rate Q_out is evaluated the same way. The simulation ends when |Q_in(tⁿ⁺¹) − Q_in(tⁿ)| < ϵ|Q_in(tⁿ)| and |Q_out(tⁿ⁺¹) − Q_out(tⁿ)| < ϵ|Q_out(tⁿ)|, where ϵ = 10⁻⁵.

Examples of the obtained velocity fields in each of the porous samples are shown in Figures 1, 2. In Figure 3 we present the typical convergence plots for Q_in and Q_out in our simulations. The initial oscillations are caused by the discontinuous initial conditions, as well as the formulation of the open boundary conditions; setting ρ_in = 1 + Δρ gives rise to a shock wave, which first travels to the outlet, and gets reflected off the outlet and inlet planes, while its amplitude decays exponentially. Once the pressure wave vanishes, the values of Q_in and Q_out continue to converge exponentially to the steady state values Q^∞_in and Q^∞_out. We expect these values to be equal, however, since the expression (12) is only first-order accurate, we observe a small difference between these two values, which is particularly evident in case of the more coarsely meshed Sample A. We have performed steady-state flow computation for Δρ ranging from 2.5 · 10⁻⁷ to 2.5 · 10⁻³, and relaxation times τ = 0.01 − 0.1, and we have obtained stable results obeying the Darcy's law, which states that under steady-state flow conditions, the flow rate Q through a given cross section is proportional to the pressure drop Δp = c²_sΔρ that drives the fluid.

FIGURE 3

Figure 3. Convergence plot for the steady-state computation in Sample A (on the left, τ = 0.02, Δρ = 10⁻⁶) and in Sample B (on the right, τ = 0.1, Δρ = 2.5 · 10⁻⁶). The exponential fits are computed in R, using robust nonlinear regression, after weighting out the initial, oscillatory phase (flow rates and time are given in lattice Boltzmann units).

These simulation results shown in Figures 1, 2 correspond to a Reynolds number approximately equal to 0.01. The velocity in the figures is expressed in lattice Boltzmann units u_LB. By dimensional analysis we can express the velocity in another arbitrary system of units according to the following transformation

$\begin{array}{cc}{u}_{SU}=u_{LB}\frac{l_{LB}}{l_{SU}}\frac{{\nu }_{SU}}{{\nu }_{LB}},& (13)\end{array}$

where the l_LB and ν_LB are the system size and kinetic viscosity in lattice Boltzmann units, l_SU and ν_SU are the same quantities expressed in another system of units. If the fluid is water, substituting the real-world parameters into Equation (13) yields the following scaling formulas for the real-world velocities in the SI system of units: u_SI = u_LB · 6000$\frac{\text{m}}{\text{s}}$ in Sample A, and u_SI = u_LB · 600$\frac{\text{m}}{\text{s}}$ in Sample B.

3.2. Dispersion in Porous Samples

Having computed the steady state flow field in each of the porous rock samples, we obtain the Lagrangian trajectories of passive tracers, which satisfy

$\begin{array}{cc}{\dot{x}}_k(t)=u(x_k(t)),& (14)\end{array}$

in the absence of molecular diffusion. The solution of the above equation is fixed by the initial position x_k(0). From these Lagrangian trajectories we study the statistics of the separation d_kj(t) = x_k(t) − x_j(t) between two fluid particles, which initially are close to each other. The separation of particles obeys the differential equation

$\begin{array}{cc}{\dot{d}}_{kj}(t)=u(x_k(t))-u(x_j(t)).& (15)\end{array}$

The constraints imposed by the porous heterogeneous structure may cause non-smooth velocity fields due to the intrinsic branching of channel connectivity, leading to broad distribution of velocities and absence of decorrelation in time. For non-smooth velocity fields we might have, for instance, that the velocity difference between two points scales as |u(x_k(t)) − u(x_j(t))| ∝ |x_k(t) − x_j(t)|^β = d^β_kj, with β < 1, as in the inertial range in turbulent flows. Multiplying the left and right hand sides terms in Equation (15) by d_kj(t) and substitution of the latter scaling lead to

$\begin{array}{cc}{\dot{d}}^2=2d_{kj}·\left[u(x_k(t))-u(x_j(t))\right]\propto d^{1+\beta },& (16)\end{array}$

whose solution for finite time gives the behavior

$\begin{array}{cc}{\left(\frac{d(t)}{d_0}\right)}^2\propto t^{\frac{2}{1-\beta }}.& (17)\end{array}$

As seen in Figure 4, we find from a flow simulation in Sample A that β ≈ 1/3, which coincides with Richardson dispersion law. Of course, such spatial scaling should not be associated to the classical Richardson-Kolmogorov cascade picture. Instead, the anomalous dispersion originates from the high degree of heterogeneity in Sample A and the broad distributions of pore sizes. This is further confirmed by the simulations performed in Sample B, here we see a particle separation in time, which is consistent with ballistic motion of particles in a more simple geometry. This is due to the fact, that the flow in Sample B is dominated by one large pore and the heterogeneity of the much smaller neighboring pores is only having a negligible impact on the averaging (see Figure 4).

FIGURE 4

Figure 4. The black line shows the change with time of the average of the squared distance of particles, which initially are closer than 1% of the system size (the linear size). In the left panel, we show data from Sample A. The green dashed line is proportional to t³ and is consistent with the Richardson dispersion law. Despite the limited scaling region, we see an anomalous dispersion, where particles separates faster than in a ballistic regime. The envelope surrounding the black line shows the standard deviation of individual particle pairs. In the right panel, we show similar data for Sample B. The dashed line is here proportional to t² and is therefore consistent with a ballistic separation of particles.

3.3. Evacuation of Particles through Advection and Diffusion

In our final experiment, we apply the previously obtained steady state velocity field u(x) in a study of pollutant removal from Sample A. In particular, we perform a comparative, qualitative study of particle evacuation rate from the porous sample, using two different flushing strategies: through continuous steady flow and through the use of rectangular pulses of equal width and amplitude, increasing period and decreasing duty cycle (see Figure 5).

FIGURE 5

Figure 5. Continuous (blue) and rectangular pulse envelope (red). The pulse width T_on is constant in all our particle evacuation simulations, and we run the simulations with three different pulse separation times T_off. The duty cycle is defined as d = T_on/T_cycle. The pulse separation time T_off is on the order of 10% of the diffusive time, defined as R²/D, where R is the characteristic pore size, and D is the diffusivity.

As the first step of this numerical experiment, we randomly seed a large number of particles N₀ (in our experiments N₀ = 250,000), uniformly distributed in the pore space of the sample. We assume that the particles are passively advected by the fluid, but also undergo molecular diffusion. The motion of such passive tracers is described by the following Langevin equation

$\begin{array}{cc}{\dot{x}}_k(t)=\overline{u}(x_k(t),t)+{\eta }_k(t),& (18)\end{array}$

where x_k(t) is the position of the kth particle at time t, and u(x_k(t), t) = u(x_k(t))ψ(t) is the fluid's velocity evaluated at particle's position at time t; ψ(t) is the pulse envelope function. Here we assume that the frequency of the pulses is smaller than the inverse of the characteristic time in which viscous effect are propagated inside the porous sample. The latter can be estimated as R²/ν, where R is on the order of the largest channel width. The final term, ${\eta }_k(t)=\sqrt{2D}\xi (x_k(t),t),$, is the random “kick” function, where D is the diffusivity and ξ_k(t) is a vector-valued Gaussian noise with mean amplitude equal to unity, whose correlation is given by 〈ξ_k(t) ξ_k′(t′)〉 = δ_kk′δ(t − t′). The above equation for passive tracers can be obtained in the limit of vanishing Stokes time from the following system of Langevin equations describing the Brownian motion of inertial particles in a moving fluid [28]

$\begin{array}{cc}\begin{array}{l}{\dot{x}}_k(t)=v_k(t),{\dot{v}}_k(t)=\frac{1}{\tau }(-v_k(t)+\overline{u}(x_k(t),t)\\ +{\eta }_k(t)),\end{array}& (19)\end{array}$

in which τ is the Stokes time. We integrate Equation (18) for each particle using the forward Euler method. Additionally we assume elastic collisions between the particles and the solid wall, and that the particles which cross the outlet plane are removed from the system. Then, we track the number N of particles remaining in the pore space as a function of the total volume of fluid flushed through the system

$\begin{array}{cc}V(t)={\displaystyle {\int }_0^t\left({\displaystyle {\int }_{A_{in}}\overline{u}(x,t)·dA}\right)d}t.& (20)\end{array}$

The results of our experiments for Sample A are presented in Figure 6. The parameters were chosen so that the effective Péclet number Pe ≈ 1600. Evidently, continuous flushing is the least efficient in terms of the amount of resources used (water, energy) in order to remove a given fraction of pollutant from the sample. Increasing the separation between the pressure pulses yields increasingly more efficient use of resources. This is due to the fact that the water flow through the sample quickly decreases the particle concentration in the main channel network, however, some pores are inaccessible to the flow, i.e., u = 0 in such spaces. The only way for the particles to be released from these pores is through diffusion, and, as it is a slower process, increasing the pulse separation increases the time for diffusion to act on. In contrast, as the particle concentration in the main channel network decreases, the efficiency of continuous flushing quickly drops.

FIGURE 6

Figure 6. The fraction of the pollutant left inside Sample A as a function of the total volume of water that passed through the sample (given as a fraction of the pore space volume), using continuous flushing (duty cycle d = 1) and increasingly separated, constant width, rectangular pressure pulses.

4. Discussion

4.1. Evaluation of the ULBM

Our choice to use the unstructured lattice Boltzmann method, as opposed to the regular grid based approaches, was dictated by the higher flexibility and accuracy of the unstructured meshes at representing complex solid boundaries. Unlike regular grids, unstructured meshes allow for locating the nodes precisely at the boundary, yielding the boundary approximation error on the order of h², as opposed to h in case of regular grids, where h is the grid spacing (or local edge length, in case of unstructured meshes). This error can be decreased further, without significantly increasing the problem size, by employing adaptive meshes, i.e., using higher precision (lower h) to resolve fine features, while keeping the mesh coarser in the regions where the boundary is almost flat, or where the flow field is approximately linear. Furthermore, unstructured meshes do not introduce extra memory overhead related to the storage of the solid phase elements, which is particularly beneficial in case of samples with low porosity. Overall, the number of volumetric elements required to represent a porous structure with precision h on the boundary scales like h⁻² log h⁻¹, as opposed to h⁻³ in case of basic regular grids. Certain techniques, such as the immersed boundary method [29] and hierarchical voxel structures [30], allow regular grid based approaches reach the effective storage efficiency similar to that of unstructured meshes. However, they increase the complexity of the fluid simulation method, and the lack of coupling between the nodal positions and the precise location of the solid boundary persists. As a consequence, boundary features finer than the minimum grid spacing cannot be represented correctly and, in the worst case scenario, the topology of the channel network might be ambiguous or inaccurate. In contrast, unstructured grids yield much simpler formulation of the boundary conditions, allow for a faithful representation of arbitrarily fine surface details, and reproduce the channel network topology exactly (although, in our case, the accuracy of the representation is still limited by the resolution of the tomographic data).

The benefits of the ULBM come at a cost of increased number of computations per grid node. At each time step, at a given node, regular grid based LBM requires a single evaluation of the collision term, one assignment (copy) operation and a single addition per each discrete direction. ULBM also requires just a single evaluation of the collision term per node, at every time step, however, the streaming and collision steps require averaging the values of f_i and Ω_i from all neighboring nodes, yielding 2(v + 1) additions and multiplications per each discrete direction, where v is the valency of the node. In the meshes used in our simulations, the average valency 〈v〉 ≈ 12 (Table 1). Finally, lattice Boltzmann methods compute the full temporal evolution of the flow and, when applied to the problem of finding the steady-state flow with constant pressure boundary condition, it converges asymptotically to the solution. Consequently, in this particular application, it can be outperformed by implicit, finite-element Stokes solvers, which typically require fewer sparse matrix multiplications (due to a lower number of variables, and a higher convergence rate). However, implicit Stokes solvers are unsuitable for simulations of unsteady flows, such as: multiphase flows, turbulent flows, flows influenced by moving solid boundaries or by oscillating pressure differences, all of which do not pose a great difficulty for the LBM. Furthermore, the LBM's convergence rate could be improved by employing more refined open boundary conditions, for example by specifying the velocity profile at the inlet, and using Grad's approximation [31] at the outlet.

5. Concluding Remarks

The complex pore geometry of most porous media is a significant obstacle for efficient numerical simulation of flow. We have here discussed the unstructured lattice Boltzmann method for fluid flow. We have used the method on a couple of rock samples and have in the sample with a more complex geometry observed an anomalous dispersion of Lagrangian particles without any particle diffusion introduced. One sample with a less complex pore space showed no sign of anomalous dispersion. We further showed that in the more complex parts of the pore space, diffusive transport might dominate over advective transport and therefore the evacuation of a pollutant might be, from a resource point of view, more effective if cycles of flushing water through the sample are used. Further simulations however would be needed to establish the most efficient strategy for pollutant removal.

Funding

Provided through the grant Earth Patterns from the Villum Foundation and by the Innovation Fund Denmark and Maersk Oil and Gas A/S through the P³ project. The Danish National Research Council (via Danscatt) provided support for the experimental work. AH acknowledges funding from the European Unions Seventh Framework Programme for research, technological development and demonstration under grant agreement no 316889. DM acknowledges funding from the People Programme (Marie Curie Actions) of the European Unions Seventh Framework Programme FP7/2007-2013/under REA grant no 297921.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank B. Vinter and J. Avery for help with improving the performance of the ULBM implementation, and S. Pedersen, K. N. Dalby and H. Suhonen for help with the experimental work at beamline ID22 at the European Synchrotron (ESRF).

References

1. de Anna P, Le Borgne T, Dentz M, Tartakovsky AM, Bolster D, Davy P. Flow intermittency, dispersion, and correlated continuous time random alks in porous media. Phys Rev Lett. (2013) 110:184502. doi: 10.1103/PhysRevLett.110.184502

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Lester DR, Metcalfe G, Trefry MG. Is chaotic advection inherent to porous media flow? Phys Rev Lett. (2013) 111:174101. doi: 10.1103/PhysRevLett.111.174101

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Kang PK, de Anna P, Nunes JP, Bijeljic B, Blunt MJ, Juanes R. Pore-scale intermittent velocity structure underpinning anomalous transport through 3-D porous media. Geophys Res Lett. (2014) 41:6184–90. doi: 10.1002/2014GL061475

CrossRef Full Text | Google Scholar

4. Benson DA, Wheatcraft SW, Meerschaert MM. Application of a fractional advection-dispersion equation. Water Res Res. (2000) 36:1403–12. doi: 10.1029/2000WR900031

CrossRef Full Text | Google Scholar

5. Mordant N, Delour J, Léveque E, Arnéodo A, Pinton JF. Long Time correlations in lagrangian dynamics: a key to intermittency in turbulence. Phys Rev Lett. (2002) 89:254502. doi: 10.1103/PhysRevLett.89.254502

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Le Borgne T, Dentz M, Carrera J. Lagrangian statistical model for transport in highly heterogeneous velocity fields. Phys Rev Lett. (2008) 101:090601. doi: 10.1103/PhysRevLett.101.090601

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Le Borgne T, Dentz M, Carrera J. Spatial Markov processes for modeling Lagrangian particle dynamics in heterogeneous porous media. Phys Rev E. (2008) 78:026308. doi: 10.1103/PhysRevE.78.026308

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Lester DR, Metcalfe G, Trefry MG. Anomalous transport and chaotic advection in homogeneous porous media. Phys Rev E. (2014) 90:063012. doi: 10.1103/PhysRevE.90.063012

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Le Borgne T, Dentz M, Villermaux E. Stretching, coalescence, and mixing in porous media. Phys Rev Lett. (2013) 110:204501. doi: 10.1103/PhysRevLett.110.204501

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Misztal MK, Hernandez-Garcia A, Matin R, Sørensen HO, Mathiesen J. Detailed analysis of the lattice Boltzmann method on unstructured grids. J Comput Phys. (2015) 297:316–39. doi: 10.1016/j.jcp.2015.05.019

CrossRef Full Text | Google Scholar

11. Cloetens P, Ludwig W, Baruchel J, Van Dyck D, Van Landuyt J, Guigay J, et al. Holotomography: quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays. Appl Phys Lett. (1999) 75:2912–4. doi: 10.1063/1.125225

CrossRef Full Text | Google Scholar

12. Müter D, Sørensen H, Jha D, Harti R, Dalby K, Suhonen H, et al. Resolution dependence of petrophysical parameters derived from X-ray tomography of chalk. Appl Phys Lett. (2014) 105:043108. doi: 10.1063/1.4891965

CrossRef Full Text | Google Scholar

13. Jha D, Sørensen HO, Dobberschütz S, Feidenhans R, Stipp SLS. Adaptive center determination for effective suppression of ring artifacts in tomography images. Appl Phys Lett. (2014) 105:143107. doi: 10.1063/1.4897441

CrossRef Full Text | Google Scholar

14. Müter D, Pedersen S, Sørensen HO, Feidenhans'l R, Stipp SLS. Improved segmentation of X-ray tomography data from porous rocks using a dual filtering approach. Comput Geosci. (2012) 49:131–9. doi: 10.1016/j.cageo.2012.06.024

CrossRef Full Text | Google Scholar

15. Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of image analysis. Nat Meth. (2012) 9:671–5. doi: 10.1038/nmeth.2089

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Si H. TetGen, a delaunay-based quality tetrahedral mesh generator. ACM Trans Math Softw. (2015) 41:11:1–11:36. doi: 10.1145/2629697

CrossRef Full Text | Google Scholar

17. Huang H, Wang L, Lu XY. Evaluation of three lattice Boltzmann models for multiphase flows in porous media. Comput Math Appl. (2011) 61:3606–17. doi: 10.1016/j.camwa.2010.06.034

CrossRef Full Text | Google Scholar

18. Pazdniakou A, Adler PM. Dynamic permeability of porous media by the lattice Boltzmann method. Adv Water Res. (2013) 62, Pt B:292–302. doi: 10.1016/j.advwatres.2013.06.001

CrossRef Full Text | Google Scholar

19. Pan C, Luo LS, Miller CT. An evaluation of lattice Boltzmann schemes for porous medium flow simulation. Comput Fluids. (2006) 35:898–909. doi: 10.1016/j.compfluid.2005.03.008

CrossRef Full Text | Google Scholar

20. Chang C, Liu CH, Lin CA. Boundary conditions for lattice boltzmann simulations with complex geometry flows. Comput Math Appl. (2009) 58:940–9. doi: 10.1016/j.camwa.2009.02.016

CrossRef Full Text | Google Scholar

21. Kang Q, Zhang D, Chen S. Simulation of dissolution and precipitation in porous media. J Geophys Res. (2003) 108:2505. doi: 10.1029/2003jb002504

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Bardow A, Karlin IV, Gusev AA. General characteristic-based algorithm for off-lattice Boltzmann simulations. Europhys Lett. (2006) 75:434–40. doi: 10.1209/epl/i2006-10138-1

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zienkiewicz OC, Codina R. A general algorithm for compressible and incompressible flow–Part I. the split, characteristic-based scheme. Int J Numer Meth Fluids. (1995) 20:869–85. doi: 10.1002/fld.1650200812

CrossRef Full Text | Google Scholar

24. Rossi N, Ubertini S, Bella G, Succi S. Unstructured lattice Boltzmann method in three dimensions. Int J Num Methods Fluids. (2005) 49:619–33. doi: 10.1002/fld.1018

CrossRef Full Text | Google Scholar

25. Talon L, Bauer D, Gland N, Youssef S, Auradou H, Ginzburg I. Assessment of the two relaxation time Lattice-Boltzmann scheme to simulate Stokes flow in porous media. Water Resour Res. (2012) 48:W04526. doi: 10.1029/2011WR011385

CrossRef Full Text | Google Scholar

26. d'Humieres D, Ginzburg I, Krafczyk M, Lallemand P, Luo LS. Multiple-relaxation-time lattice boltzmann models in three dimensions. Philos Trans Math Phys Eng Sci. (2002) 360:437–51. doi: 10.2307/3066323

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Chikatamarla SS, Karlin IV. Entropic lattice Boltzmann method for turbulent flow simulations: Boundary conditions. Phys A Statis Mech Appl. (2013) 392:1925–30. doi: 10.1016/j.physa.2012.12.034

CrossRef Full Text | Google Scholar

28. Cardy J, Falkovich G, Gawedzki K, Nazarenko S, Zaboronski OV. Non-equilibrium Statistical Mechanics and Turbulence. Cambridge: Cambridge University Press (2008). doi: 10.1017/CBO9780511812149

CrossRef Full Text | Google Scholar

29. Chen DJ, Lin KH, Lin CA. Immersed boundary method based lattice boltzmann method to simulate 2d and 3d complex geometry flows. Int J Modern Phys C. (2007) 18:585–94. doi: 10.1142/S0129183107010826

CrossRef Full Text | Google Scholar

30. Stiebler M, Tölke J, Krafczyk M. Advection-diffusion lattice Boltzmann scheme for hierarchical grids. Comput Math Appl. (2008) 55:1576–84. doi: 10.1016/j.camwa.2007.08.024

CrossRef Full Text | Google Scholar

31. Chikatamarla SS, Ansumali S, Karlin IV. Grad's approximation for missing data in lattice Boltzmann simulations. Europhys Lett. (2006) 74:215. doi: 10.1209/epl/i2005-10535-x

CrossRef Full Text | Google Scholar

32. Ayachit U. The ParaView Guide: A Parallel Visualization Application. Kitware (2015).