Hydrodynamic Dispersion in Porous Media and the Significance of Lagrangian Time and Space Scales

Abstract

Transport in porous media is critical for many applications in the environment and in the chemical process industry. A key parameter for modeling this transport is the hydrodynamic dispersion coefficient for particles and scalars in a porous medium, which has been found to depend on properties of the medium structure, on the dispersing compound, and on the flow field characteristics. Previous studies have resulted in suggestions of different equation forms, showing the relationship between the hydrodynamic dispersion coefficient for various types of porous media in various flow regimes and the Peclet number. The Peclet number is calculated based on a Eulerian length scale, such as the diameter of the spheres in packed beds, or the pore diameter. However, the nature of hydrodynamic dispersion is Lagrangian, and it should take the molecular diffusion effects, as well as the convection effects, into account. This work shifts attention to the Lagrangian time and length scales for the definition of the Peclet number. It is focused on the dependence of the longitudinal hydrodynamic dispersion coefficient on the effective Lagrangian Peclet number by using a Lagrangian length scale and the effective molecular diffusivity. The lattice Boltzmann method (LBM) was employed to simulate flow in porous media that were constituted by packed spheres, and Lagrangian particle tracking (LPT) was used to track the movement of individual dispersing particles. It was found that the hydrodynamic dispersion coefficient linearly depends on the effective Lagrangian Peclet number for packed beds with different types of packing. This linear equation describing the dependence of the dispersion coefficient on the effective Lagrangian Peclet number is both simpler and more accurate than the one formed using the effective Eulerian Peclet number. In addition, the slope of the line is a characteristic coefficient for a given medium.

1. Introduction

Hydrodynamic dispersion through porous media has long been of great importance in many engineering fields. For instance, in the field of petroleum engineering, chemicals such as surfactants are injected into hydrocarbon reservoirs to enhance oil recovery [1]. The theory of hydrodynamic dispersion can assist engineers to estimate how far the chemicals reach, thereby designing appropriate injection points to cover the targeted zone. In environmental applications, dispersion in porous media has helped to predict how pollutants spread in the ground water and contaminate aquifers [2,3,4].

What is hydrodynamic dispersion? In general, substances or particles that disperse in a porous medium do so because of two effects: molecular diffusion and convective dispersion (or mechanical dispersion). The transport that combines both mechanisms of dispersion is known as hydrodynamic dispersion [3,5,6]. The effectiveness of solute transport is based on the hydrodynamic dispersion coefficient, and many researchers have worked to determine this coefficient for various porous geometries at different flow regions by using both experimental and numerical methods. All studies have agreed that the hydrodynamic dispersion coefficient is the sum of an effective molecular diffusivity coefficient and a mechanical dispersion coefficient. Furthermore, the coefficient for effective molecular diffusion in porous media $D_m^{\prime }$ is always smaller than that for molecular diffusion in a pure solvent $D_m$ [3,4,5,6,7,8].

There have been several suggestions for calculating the effective molecular diffusion coefficient in a porous medium. Several researchers have proposed to estimate this parameter based on the diffusive tortuosity, which is defined by the squared ratio of the average length of the path traveled by a substance as it diffuses through the porous medium over the straight-line length ${\tau }_d={\left(\frac{L_d}{L_s}\right)}^2$ [9]. The diffusive tortuosity is often defined as the molecular diffusion coefficient of a species in a free fluid relative to that in a porous medium ${\tau }_d=\frac{D_m}{D_m^{\prime }}$ [9,10,11,12,13,14]. However, others have argued that the porosity ε should be present in the equation of tortuosity as ${\tau }_d=\epsilon \frac{D_m}{D_m^{\prime }}$ [15,16,17]. From another perspective, other researchers have found that there is a link between the diffusion in a porous medium and its electrical conductivity, proposing a correlation such as $\frac{1}{F\epsilon }=\frac{D_m}{D_m^{\prime }}$, where F is the formation electrical resistivity factor [6,18,19]. At the end of the day, it appears that the most accurate method is to experimentally measure the diffusive tortuosity of a porous medium by measuring the diffusion coefficient, both in the bulk fluid and in the porous medium [20].

Regarding the hydrodynamic dispersion coefficient, a vast amount of studies for different types of porous media have ended up with different equation forms. In some studies, the longitudinal and transverse dispersion coefficient have been expressed in the general form:

$$D_L=D_m^{\prime }+{\alpha }_L{u}^n$$

(1)

$$D_T=D_m^{\prime }+{\alpha }_T{u}^n$$

(2)

where n is an empirical constant, varying from 1 to 2; u is the pore velocity; and ${\alpha }_L$ and ${\alpha }_T$ represent the effects of geometry in the longitudinal and transverse dispersion, respectively. Other researchers have worked with beds packed with spheres and found different results for different flow regimes; these are briefly summarized in

Table 1. Summary of studies on hydrodynamic dispersion through packed beds in the randomly packed spheres (RPS) configuration.

Reference	Conditions	Equations of Hydrodynamic Dispersion Coefficient
Fried-Combranous [21]	$P{e}_m^E<300$ $P{e}_m^E<{10}^5$	$\frac{D_L}{D_m}=\frac{1}{{\tau }_d}+0.5P{e}_m^E{}^{1.2}$$\frac{D_L}{D_m}=\frac{1}{{\tau }_d}+\left(1.8\pm 0.4\right)P{e}_m^E$
Hiby [5,22]	$Re<100$	$\frac{D_L}{D_m}=\frac{1}{{\tau }_d}+\frac{0.65P{e}_m^E}{1+7P{e}_m^E{}^{-0.5}}$
Harleman et al. [22]	$20<P{e}_m^E<4000$	$\frac{D_L}{D_m}=0.67+0.66P{e}_m{}^{1.2}$
Koch and Brady [7]	Impermeable packed bed ${\left(1-\epsilon \right)}^2<P{e}_m^{\prime E}/2<1$ $P{e}_m^{\prime E}/2>1$	$\frac{D_L}{D_m^{\prime }}=1+\frac{3}{4}\frac{P{e}_m^{\prime E}}{2}$ $\frac{D_L}{D_m^{\prime }}=1+\frac{3}{4}\frac{P{e}_m^{\prime E}}{2}+\frac{{\pi }^2}{6}\left(1-\epsilon \right)\frac{P{e}_m^{\prime E}}{2}\ln \left(\frac{P{e}_m^{\prime E}}{2}\right)$
Delgado et al. [5]	Longitudinal dispersion
	Diffusion regime $(P{e}_m^E<0.1)$	$\frac{D_L}{D_m^{\prime }}=1$
	Predominant diffusional regime $(0.1<P{e}_m^E<4)$	$\frac{D_L}{D_m^{\prime }}=\frac{P{e}_m^{\prime E}}{0.8/P{e}_m^{\prime E}+0.4}$
	Predominant mechanical dispersion ($4<P{e}_m^E$ and $Re<10$)	$\frac{D_L}{D_m^{\prime }}=\frac{P{e}_m^{\prime E}}{\sqrt{18P{e}_m^{\prime E}{}^{-1.2}+2.35S{c}^{-0.38}}}$
	Pure mechanical dispersion ($10<Re$ and $P{e}_m^E<{10}^6$)	$\frac{D_L}{D_m^{\prime }}=\frac{P{e}_m^{\prime E}}{25S{c}^{1.14}/P{e}_m^{\prime E}+0.5 }$
	Dispersion out of Darcy domain $(P{e}_m^E>{10}^6$)	$\frac{D_L}{D_m^{\prime }}=\frac{P{e}_m^{\prime E}}{2}$
	Transverse dispersion
	Diffusion regime ($P{e}_m^E<1$)	$\frac{D_T}{D_m^{\prime }}=1$
	Predominant mechanical dispersion ($1<P{e}_m^E<1600$)	$\frac{D_T}{D_m^{\prime }}=1+\frac{1}{2.7\times {10}^{-5}Sc+12/P{e}_m^{\prime E}}$
	Pure mechanical dispersion ($1600<P{e}_m^E<{10}^6$)	$\frac{D_T}{D_m^{\prime }}=\frac{P{e}_m^{\prime E}}{\left(0.058Sc+14\right)-\left(0.058Sc+2\right)\exp \left(-500S{c}^{0.5}P{e}_m^{\prime E}\right)}$
	Dispersion out of Darcy domain ($P{e}_m^E>{10}^6$)	$\frac{D_T}{D_m^{\prime }}=\frac{P{e}_m^{\prime E}}{12}$

[5,7,21,22]. The longitudinal hydrodynamic dispersion coefficient has been found to depend on the Eulerian Peclet number $P{e}_m^E$ determined by the diameter of the spheres $d_p$ and molecular diffusion in the bulk fluid $D_m$ [5,6], as follows:

$$P{e}_m^E=u{d}_p/D_m.$$

(3)

Koch and Brady suggested equations for both the longitudinal and transverse dispersion coefficient as a function of an effective Eulerian Peclet number defined with the use of the effective diffusivity, $P{e}_m^{\prime E}$,

$$P{e}_m^{\prime E}=u{d}_p/D_m^{\prime }$$

(4)

by analytically solving the Stokes flow equations in packed beds consisting of randomly packed spheres [7].

Delgado et al. collected their own data and combined them with experimental data available in the literature to characterize the dispersion in packed beds based on the effective Eulerian Peclet number at different flow regimes [5].

It is commonly accepted that the Peclet number is estimated based on the diameter of the spheres that constitute the porous media packing, d_p. The sphere diameter is an Eulerian length scale. It has the advantage of being known a priori. However, the nature of dispersion is Lagrangian, so the hypothesis here is that a Lagrangian length scale should be used in order to combine the molecular diffusion effects and the convection effects in dispersion. In the present study, we set up model equations for the dispersion coefficient with the effective Peclet number based on either the Eulerian or the Lagrangian length scale and examine their predictive ability. The porous media that were investigated in this study included the face-centered cubic (FCC) packing of spheres and randomly packed spheres (RPS), and we then extended to the body-centered cubic (BCC) sphere arrangement for validation purposes. The images of these three packing types are shown in Figure 1. The lattice Boltzmann method (LBM) was used to simulate the flow of water at room temperature through the given porous media, and then Lagrangian particle tracking (LPT) was applied to follow particle trajectories in space and time. The transport of nanoparticles (NPs) of different sizes, corresponding to different dimensionless Schmidt numbers, was investigated. Since all prior studies have confirmed that the dispersion in the streamwise direction (flow direction) is far larger than the one in spanwise direction [3,5], this paper mainly focused on longitudinal dispersion. The contributions of the present work are (a) a comparison of Eulerian and Lagrangian scales for dispersion in porous media; (b) the introduction of a Lagrangian effective Peclet number that is found to be more appropriate for the description of dispersivity, thus unifying hydrodynamic dispersion across porous media types; and (c) the development of model equations that relate longitudinal dispersivity to the Peclet number.

2. Materials and Methods

2.1. Lattice Boltzmann Method

In order to simulate the flow through different types of porous media, we used an in-house code based on the LBM. In the LBM, one applies the discretized Boltzmann equation to calculate the probability of fluid particles to move and distribute in different pre-determined directions. The method is mesoscopic, where the macroscopic properties of the flow and the fluid, such as velocity and density, are accurately recovered [23,24,25]. It is important to highlight that the fluid particles in the LBM can move only in certain directions on a lattice, and these directions are defined by the selected lattice velocity model. In general, an LBM lattice is specified as DmQn, where m indicates the dimensionality (2 or 3 for a two-dimensional or a three-dimensional space, respectively) and n indicates the number of lattice directions that particles are permitted to travel. Hence, the first step in applying the LBM is to divide the computational domain geometry into structured cartesian meshes, and each node in the lattice is represented by a binary value; for instance, the nodes inside the solid phase of a porous medium are described as ‘TRUE,’ whereas the ones in the pore space are designated as ‘FALSE’ [26]. By adopting this discretization of the porous medium space, one may lose fidelity close to the solid boundaries. On the other hand, Cartesian meshing with the LBM is very convenient when simulating the flow field in complicated geometries, like the pore space in digital rock images, where the generation of unstructured grids is required for each porous medium realization without using LBM. Other techniques for simulating flow fields have been recently developed, which demonstrate high order of accuracy and can resolve flow velocities at points very close to the solid grains [27,28]. In Appendix A, we present details of the grid resolution analysis for the LBM to justify the use of the method, in addition to its other computational advantages.

The LBM is relatively simple to implement computationally, and LBM algorithms are parallelizable due to the application of the discrete Boltzmann equation to each single fluid node. Each node needs information from its neighboring nodes in order to calculate particle distribution functions, making the algorithm a good candidate for Message Passing Interface (MPI) [23,25,29]. The key term in this method is the particle distribution function f, which represents the fraction of particles at a fluid node position at a certain time. Basically, there are three main steps for solving the discretized Boltzmann equation, i.e., the streaming, the collision, and the forcing steps, as follows:

$$\underset{Streaming}{\underbrace{f_i\left(\overrightarrow{x}+{\overrightarrow{e}}_i\Delta t,t+\Delta t\right)=f_i\left(\overrightarrow{x},t\right)}}+\underset{Collision}{\underbrace{{\Omega }_i\left(\overrightarrow{x},t\right)}}\pm \underset{Forcing}{\underbrace{f{f}_i}}.$$

(5)

In the streaming step, the particle distribution function $f_i$ at the position $\overrightarrow{x}$ at time t moves along the direction i with velocity ${\overrightarrow{e}}_i$ to the new position $(\overrightarrow{x}+{\overrightarrow{e}}_i\Delta t$) during one time step $\Delta t$. The collisions occurring during streaming steps are demonstrated by the collision operator ${\Omega }_i$ [2] that represents the change of $f_i$ due to the collision among particles. In the present work, the Bhatnagar, Gross, and Krook (BGK) approximation for the single-relaxation time was used to estimate the collision operator. The BGK has been a simple and reliable relaxation time approximation, and it is given by:

$${\Omega }_i\left(\overrightarrow{x},t\right)=-\frac{1}{{\tau }_r}\left(f_i-f_i^{eq}\right).$$

(6)

The term ${\tau }_r$ is the relaxation time towards the equilibrium distribution function $f_i^{eq}$, and it relates to the fluid viscosity $\mathsf{\nu }$ as follows:

$$\mathsf{\nu }=\frac{1}{3}\left({\tau }_r-\frac{1}{2}\right).$$

(7)

The particle equilibrium distribution function $f_i^{eq}$ is calculated by:

$$f_i^{eq}\left(\overrightarrow{x}\right)=w_i\rho \left[1+\frac{{\overrightarrow{e}}_i.\overrightarrow{U}}{c^2}+\frac{{({\overrightarrow{e}}_i.\overrightarrow{U})}^2}{2c^4}-\frac{{\overrightarrow{U}}^2}{2c^2}\right]$$

(8)

where $w_i$ is the lattice specific weighting factor for the lattice direction i, $\rho$ is the density, $\overrightarrow{U}$ is the macroscopic velocity, and c is the speed of sound. The last component added in the Boltzmann equation is the forcing factor to account for the pressure drop during streaming steps [23,25,29].

As mentioned above, the main purpose of using the LBM is to calculate the macroscopic properties such as the fluid density $\rho$ and macroscopic velocity $\overrightarrow{U}$; this requirement can be achieved by applying the conservation equations of mass and momentum, as follows:

$$\rho ={\displaystyle \sum }_{i=0}^n{f}_i$$

(9)

$$\rho \overrightarrow{U}={\displaystyle \sum }_{i=0}^n{f}_i{\overrightarrow{e}}_i$$

(10)

where n is the number of lattice directions that particles are permitted to move, including the rest position. In our work, D3Q15 was selected as the lattice velocity model with 15 allowable lattice directions. Moreover, the periodic boundary condition was used in three directions x, y, and z; and at the solid–fluid interfaces, the no-slip boundary condition with the bounce-back technique was applied [24,30].

2.2. LBM Code Validation and Verification

In order to verify the LBM code, we calculated the permeability of periodic simple cubic (SC), FCC, or BCC porous medium and then compared our results with those of two other groups. The size of the simulation box used for each geometry was 100 × 100 × 100 μm, and the domain consisted of 201 nodes in each x, y, and z direction. The LBM was used to run for different values of the forcing factor (the pressure drop over the domain length), resulting in different velocity fields. Next, Darcy’s empirical law, showing the relation between the pressure drop $\frac{\Delta P}{L}$ through the flow and the superficial fluid velocity ${\mathrm{U}}_s$ was applied to determine the permeability $k$ of different packing arrays [31]. Darcy’s law is expressed as:

$${\mathrm{U}}_s=-\frac{k}{\mu }\frac{\Delta P}{L}$$

(11)

where k is the medium permeability and μ is the dynamic viscosity of the fluid. The k value was then divided by the sphere diameter squared to obtain a dimensionless permeability. Our results for the dimensionless permeability of different geometries, including SC, FCC, and BCC, compared well to the ones of Chapman and Higdon [32] that were obtained by solving the unsteady Stokes equations for the microscopic flow. Moreover, Eshghinejadfard et al. [33] used their own LBM code to calculate the permeability of FCC and BCC and showed similar results to ours. The details of the comparison are shown in

Table 2. Comparison of the results for dimensionless permeability with prior results. The reported error is the relative absolute error compared to references.

Geometry	Porosity	Dimensionless Permeability (k/dp2)			Relative Error (%)Compared to
		Present Results	Reference (1) [32]	Reference (2) [33]	Reference (1) [32]	Reference (2) [33]
FCC	0.26	1.75 × 10−4	1.74 × 10−4	1.70 × 10−4	0.76%	2.96%
BCC	0.32	5.10 × 10−4	5.02 × 10−4	5.10 × 10−4	1.49%	0.10%
SC	0.48	2.61 × 10−3	2.53 × 10−3		3.10%

.

Moreover, the LBM code was verified for laminar flow through different sphere packing types (SC, BCC, and FCC) of fixed-bed columns. This flow is usually modeled by the Blake–Kozeny (BK) equation [31], which shows the correlation between the pressure drop over the length with superficial velocity:

$${\mathrm{U}}_s=\frac{\Delta P}{L}\frac{d_p^2}{150\mu }\frac{{\epsilon }^3}{{\left(1-\epsilon \right)}^2}.$$

(12)

The results are presented in

Table 3. The lattice Boltzmann method (LBM) code verification by Blake–Kozeny equation. The reported error is the relative absolute error for the calculated superficial velocity.

Parameter	SC	BCC	FCC	RPS
Porosity, $\epsilon$	0.48	0.32	0.26	0.37
Spheres’ diameter, dp (cm)	1.00 × 10−3	8.66 × 10−3	7.07 × 10−3	1.41 × 10−3
Pressure drop, $\Delta P/L$ (g/cm2 s2)	3.66 × 103	1.67 × 104	5.93 × 104	5.33 × 105
Pore velocity (cm/s)	2.00 × 10−1	2.00 × 10−1	2.00 × 10−1	2.00 × 10−1
Superficial velocity, $U_s$ (cm/s)	9.53 × 10−2	6.40 × 10−2	5.19 × 10−2	7.42 × 10−2
Blake–Kozeny velocity (cm/s)	9.62 × 10−2	5.92 × 10−2	6.31 × 10−2	9.08 × 10−2
Relative error	0.91%	8.08%	17.70%	18.26%

, revealing that SC followed the BK correlation with a less than 1% error, while other types of packing deviated from the BK equation with errors of 8.08% and 17.7% for BCC and FCC, respectively. In prior work out of our laboratory, this code was also validated against flow cases with known analytical solutions, such as flow in a slit, in a pipe, and flow around an infinite array of spheres [29].

2.3. Lagrangian Particle Tracking

LPT is an effective method to track the trajectories of particles in a Lagrangian framework while they move in a flow field. However, there are several assumptions when this method is applied. Firstly, the nanoparticles (NPs) are considered as passive, so the presence of NPs does not modify the flow field. Secondly, the NP concentration is low, so the particle–particle interactions are not taken into consideration. Thirdly, the NP–wall interactions and the adsorption of NPs on the wall of the packing spheres are neglected. In short, there are just two kinds of transport contributing to the motion of the NPs—the convection by the flow field and the molecular diffusion of NPs. Therefore, the equation of motion of the NPs can be expressed by:

$${\overrightarrow{X}}_{t+\Delta t}={\overrightarrow{X}}_t+\Delta t.{\overrightarrow{U}}_t+\Delta {\overrightarrow{X}}_m$$

(13)

where ${\overrightarrow{X}}_t$ is the position of a NP at time t and ${\overrightarrow{X}}_{t+\Delta t}$ is the new position of that NP in the next time step ($t+\Delta t)$ [24,26]. The velocity of the nanoparticle ${\overrightarrow{U}}_t$, which leads to the convective transport, can be obtained from the LBM flow field using a trilinear interpolation scheme to calculate the velocity of the NPs,${\overrightarrow{U}}_t$, at its position between the grid nodes [2,24,30]. The second motion $\Delta {\overrightarrow{X}}_m$, the molecular diffusion, of the particles is estimated by Einstein’s theory for Brownian motion. After moving because of convection, each particle undertakes a Brownian motion jump that is calculated according to a random number drawn from a normal distribution with a zero mean (because the Brownian motion is equally probable to be in the positive or the negative direction) and a standard deviation $\sigma$ determined by the NP diffusivity $D_m$. The standard deviation in each space direction is found as follows:

$$\sigma =\sqrt{2D_m\Delta t}=\sqrt{\frac{2\mathsf{\nu }\Delta t}{Sc}}$$

(14)

where the Schmidt number Sc is the ratio of the fluid kinematic viscosity $\mathsf{\nu }$ and the NP molecular diffusivity [31]. At this point, by estimating the travelling distance of all NPs caused by convection and Brownian motion, the position of each particle can be tracked at every single time step. In the case that the NPs move into the solid packing after a time step, they are bounced back to their current position in the fluid phase. To be more specific, the new position of each particle in every single time step was checked, and if it was inside the grains, the new position was not updated. The time step in LPT was chosen so that NPs never moved more than one-half of the LBM mesh unit. We determined the timestep by ensuring that any LPT particle that moved with the maximum fluid velocity in every time step plus the maximum Brownian motion in the same direction did not move more than half of the lattice unit. This algorithm has been validated with Taylor–Aris diffusion theory and has shown excellent results [24].

2.4. Velocity Autocorrelation Function

A Lagrangian length scale is determined based on the relevant Lagrangian time scale and the average velocity of all NPs, which is equal to the average pore velocity of $l^L={\tau }^L{u}_p$. Taylor defined the Lagrangian time scale as [34]:

$${\tau }^L={{\displaystyle \int }}_0^{\infty }{R}^L\left(\tau \right)d\tau$$

(15)

where $R^L\left(\tau \right)$ is the velocity autocorrelation function (VACF). While Taylor’s work was focused on turbulent velocity fluctuations along the trajectories of fluid particles, the analogous VACF for porous media would be based on the velocity fluctuations of the NP velocities relative to the average NP velocity at each time step. Furthermore, since the NPs move with both Brownian motion and convection, the VACF should be the material autocorrelation function, as proposed by Saffman [35]. In that work, it was recognized that a scalar marker does not follow the same path as a fluid particle, since it can move off a streamline because of Brownian diffusion. The VACF would need to be defined with the velocity of a diffusing particle, instead of the velocity of a fluid particle, as follows:

$$R^L\left(t,t_0\right)=\frac{\overline{V_j^{\prime }\left(t\right)V_j^{\prime }\left(t_0\right)}}{{\overline{V_j^{\prime 2}\left(t\right)}}^{1/2}{\overline{V_j^{\prime 2}\left(t_0\right)}}^{1/2}}$$

(16)

where $t_0$ is the initial time, and $V_j^{\prime }\left(t\right)$ is the velocity fluctuation of NP j, $V_j^{\prime }\left(t\right)=V_j\left(t\right)-{\overline{V}}_j\left(t\right)$. This equation has been used for dispersion in turbulent flow to estimate the Lagrangian time scale [36,37,38,39,40]. At small times after the NP release, the VACF is almost equal to unity, because the initial velocity is correlated with itself; however, over time, the VACF drops to zero. This happens because the random molecular movement of the particles takes them away from their original velocity streamlines, and the velocity of NPs as time advances has no correlation with the initial velocity of the NPs [34,36,37]. Hence, the Lagrangian time scale represents the time that the particles need to “forget” the velocity of the location from where they came, and it is appropriate to be used to cover the effect of molecular diffusion in the dispersion.

2.5. Scope of Work

Numerical experiments were conducted for two types of porous media, with infinite arrays of spheres packed in the FCC and RPS configurations. The computational domain size was the same for both types of porous media (100 × 100 × 100 μm), but the diameter of the spheres was different. The FCC domain was the minimum periodic cell for the FCC configuration, which contained four spheres with a diameter of 70.76 μm, whereas the simulation box for randomly packed spheres consisted of 432 spheres of a 14.06 μm diameter each. In the RPS configuration, 432 spheres were packed in a cubic domain by the Lubachevsky–Stillinger simulation algorithm [30,41,42]. The images of two packing types are shown in Figure 1.

The number and the size of the spheres in RPS were chosen to ensure that the simulation domain was representative of the porous medium (i.e., that the domain was at least as large as a representative elementary volume for the packed bed). This means that if the size of the RPS simulation box was larger and contained more spheres, the porous media properties (e.g., porosity, tortuosity, and permeability) would remain unchanged. The method used to determine the number of spheres needed to obtain a representative RPS volume included the generation of domains with more spheres keeping the same porosity. After completing flow runs, the average fluid velocity was calculated within a small cube in the center of the domain, and this calculation was repeated by increasing the calculations from the center to the full domain. In Figure 2, we plot the average fluid velocity obtained with the LBM in the streamwise (x direction) and spanwise directions (y and z directions) as a function of increasing the size of the cube from the center to the full length of the cubic domain. With the selected simulation box for RPS configuration, the average velocities in all directions were stable with the increment in the domain size. Since the pressure drop was the same, this means that the permeability of the porous medium did not change with domain size, as seen from Darcy’s law. Similar results were achieved when we changed the direction of the flow in the same RPS configuration, setting the flow direction as the y direction and then the z direction. Therefore, the created RPS computational domain was representative of an isotropic medium with RPS packing.

The simulation boxes for FCC and RPS were discretized into LBM lattice nodes. The number of lattice nodes greatly affected the accuracy of the results, as did the computational time. Details of the grid independence analysis that was conducted are presented in Appendix A. Finally, the number of nodes was chosen to balance accuracy with computational resources so that the FCC domain was meshed with a 201 × 201 × 201 mesh, while a 401 × 401 × 401 mesh was applied for the RPS domain. Next, the LBM code in conjunction with a D3Q15 lattice velocity model was implemented to obtain the flow field of water at room temperature through porous media, with the pore velocities 0, 50, 100, 200, 500, 1000, 2000, 3000, and 4000 μm/s. Next, 50,000 NPs were released in the flow domain at randomly and uniformly selected positions in the pore space, and the velocity and trajectory of each of these particles were tracked using the LPT algorithm, as already discussed. For each velocity field, different LPT runs with different values of the particle Schmidt number (Sc = 100, 1000, 7080, and 10,000) were performed, providing data for Sc covering two orders of magnitude. Moreover, the range of the Schmidt number was chosen to cover particles with diameters from 0.22 to 5 nm, and Sc = 7080 was the Schmidt number of Janus particles (d_p = 3.6 nm) in water, which has been investigated in prior work to lower oil–water interfacial tension [1]. Thus, thirty-six LBM/LPT simulations were conducted for each one of the two porous media configurations.

The resulting data were used to estimate the VACF and the Lagrangian timescale, as mentioned previously, with initial velocity obtained at the beginning of the LPT run ($t_o=0)$. After that, the hydrodynamic dispersion coefficient was determined based on the positional variance method [24,43]:

$$D=\underset{t\to \infty }{\lim }\frac{1}{2}·\frac{d{\sigma }^2}{dt}$$

(17)

where the position variance ${\sigma }^2$ is calculated in each direction as:

$${\sigma }^2=\overline{{\left(X_j-\overline{X_j}\right)}^2}.$$

(18)

The dispersion coefficient, calculated for the flow velocity equal to 0 was considered as the effective molecular dispersion $D_m^{\prime }$.

3. Results

3.1. Velocity Autocorrelation Function

Figure 3 is a presentation of the VACF of the streamwise velocity in different cases, when water at room temperature flowed through a porous bed packed with an infinite array of spheres arranged in the FCC geometry. The VACF was calculated based on the initial velocity at the beginning of the LPT run (t_o = 0), which started once the flow field had reached a steady state. Figure 3 consists of four graphs—Figure 3a–d—that represent the VACF of four different types of diffusion NPs with Schmidt numbers of 100, 1000, 7080, and 10,000, respectively. For each value of Sc, the velocity autocorrelation functions for various values of average pore velocity, ranging from 50 to 4000 μm/s, are illustrated. Similarly, Figure 4 is a demonstration of the results for the geometry of randomly packed spheres and at the same simulation conditions as Figure 3. In Figure 5 the effect of molecular diffusion on the Lagrangian timescale is highlighted. In other words, the dependence of the VACF on Sc, for both the FCC Figure 5a and RPS geometries Figure 5b, at the same average pore velocity of 1000 μm/s, is highlighted.

3.2. Hydrodynamic Dispersion Coefficient

All the data from numerical experiments related to the hydrodynamic dispersion coefficient and Lagrangian timescale for flow of NPs through FCC- and RPS-packed beds are summarized in Appendix B. There were, in total, thirty-six experiments carried out for each porous bed configuration. Four types of NPs with Sc = 100, 1000, 7080, and 1000 were released in the flow field with an average pore velocity of 0, 50, 100, 200, 500, 1000, 2000, 3000, and 4000 μm/s. The value of the dispersion coefficient in the static flow field (u = 0) was used as the effective diffusivity $D_m^{\prime }$ in order to estimate the ratio of dispersion over diffusion ($D_L/D_m^{\prime })$. In addition, the effective Peclet number was obtained from the Eulerian length scale (diameter of the packed spheres, see Equation (4)) or the Lagrangian length scale ($l^L=u{\tau }^L)$, as follows:

$$P{e}_m^{\prime L}=\frac{u{l}^L}{D_m^{\prime }}=\frac{u^2{\tau }^L}{D_m^{\prime }}.$$

(19)

Figure 6 is a plot of the ratio of the hydrodynamic dispersion coefficient over the effective molecular diffusion coefficient for particles in the FCC-packed bed as a function of the effective Peclet number. Figure 6a, b displays the dependence of the coefficient ratio on the effective Eulerian Peclet number and the effective Lagrangian Peclet number, respectively. The examined ratio is linearly correlated with the effective Lagrangian Peclet number. The same analysis was done for the randomly-packed bed; the results are reported in Figure 7.

4. Discussion

4.1. Velocity Autocorrelation Function and Lagrangian Timescale

Figure 3, Figure 4 and Figure 5 show that the velocity autocorrelation coefficients started from unity and dropped to zero in different patterns. The VACF dropped depending on the Schmidt number, average velocity, and geometry. Therefore, the Lagrangian timescale also changed with these parameters because it was calculated as the area formed by the velocity correlation and the two coordinate axes x and y (see Equations (15) and (16)). Figure 3 and Figure 4 demonstrate that the Lagrangian timescale decreased when the velocity increased. This means that particles in high pore velocity cases tend to drop off the flow lines quicker at a constant molecular diffusion coefficient. At the same time after the NP release, the NPs in high pore velocity simulations traveled farther in the streamwise direction and had the opportunity to experience a longer tortuous path than NPs in lower pore velocity cases. Moreover, from Figure 5, it is firmly seen that the higher the Schmidt number was, the larger the Lagrangian timescale was. This observation can be explained considering that, when the Schmidt number was high, the random movement because of molecular diffusion was small relative to the movement due to convection. Hence, it took a longer time for the NPs to completely stop following the flow streamlines and for the velocity autocorrelation coefficient to drop to zero. In addition, when the same type of fluid and NPs flowed through the porous media with the same average velocity, the characteristic time scale in the case of the FCC-packed bed was much larger than that in the RPS bed. The reason for this is that the FCC configuration was less tortuous than the RPS configuration. When the particles traveled in the porous media, not only random jumps of particles by diffusion but also their movement due to interaction with the spheres in the packed bed, made them move away from their streamline trajectories. Hence, the particles, moving in more tortuous paths, interacted with the wall more often. This could explain why they quickly decorrelated with their initial velocity.

There are two important messages about the Lagrangian scale that should be highlighted. First, Figure 3 and Figure 4 show that with the same geometry, the Lagrangian timescale changed with molecular diffusivity (or Schmidt number) and flow velocity. Hence, the relevant length scale also varied with Schmidt number and flow velocity. In contrast, for the same geometry, the Eulerian length scale was unchanged, because it only related to the geometrical configuration of the porous medium. Second, by comparing the change of Lagrangian timescale in Figure 3 and Figure 4 with that in Figure 5, it can be seen that the effects of the Schmidt number (molecular diffusion) were more prominent than the effects of flow velocity and geometry.

4.2. Hydrodynamic Dispersion Coefficient

By examining Figure 6 and Figure 7, it is seen that the ratio of the dispersion coefficient and the diffusion coefficient is a function of Pe. Figure 6a and Figure 7a indicate that the ratio does not linearly depend on the Eulerian Peclet number $P{e}_m^{\prime E}$. A linear equation is a good fit to the correlation of the dispersion and diffusion ratio with Lagrangian Peclet number $P{e}_m^{\prime L}$ for both types of packing (as seen in Figure 6b and Figure 7b). This finding is applicable for an Eulerian Peclet number below 4035 for FCC and 860 for RPS, as well as a Lagrangian Peclet number up to 7500 for FCC and 613 for RPS.

As the relationship between the dispersion coefficient and Eulerian Peclet number is not clear, to correlate the dispersion and diffusion ratio with the Eulerian Peclet number, we suggested two equations that give a minimum value for the mean relative error for this case. Following Taylor–Aris dispersion through a tube [44,45], in the first correlation scheme, a quadratic equation was applied. However, a quadratic equation could not represent well the whole data set for both types of packing. Hence, two equations, including one for low Pe and one for high Pe, were used for each type of packing. In the second correlation scheme, the form of Equation (1) was utilized, which also showed a good fit in this situation. While there was uncertainty in selecting the appropriate equation in the case of Eulerian Pe, only the linear equation was found to be suitable to correlate the dispersion and diffusion ratio with the Lagrangian Peclet number. This was named “Correlation Scheme 3”. The three schemes are well-presented in

Table 4. Three correlation schemes for hydrodynamic dispersion of particle transport in an FCC-packed bed.

Correlation Scheme	Correlation Equation	Mean Relative Error
1	$P{e}_m^{\prime E}<700 \frac{D_L}{D_m^{\prime }}=0.000343{\left(P{e}_m^{\prime E}\right)}^2+1$ $P{e}_m^{\prime E}>700 \frac{D_L}{D_m^{\prime }}=0.000168 {(P{e}_m^{\prime E})}^2+0.112P{e}_m^{\prime E}+1$	6.3%
2	$\frac{D_L}{D_m^{\prime }}=1+0.00126{\left(P{e}_m^{\prime E}\right)}^{1.774}$	3.3%
3	$\frac{D_L}{D_m^{\prime }}=1+0.423P{e}_m^{\prime L}$	3.3%

for FCC and

Table 5. Three correlation schemes for hydrodynamic dispersion of particles transport in RPS-packed bed.

Correlation Scheme	Correlation Equation	Mean Relative Error
1	$P{e}_m^{\prime E}<50 \frac{D_L}{D_m^{\prime }}=0.00853{(P{e}_m^{\prime E})}^2+1$ $P{e}_m^{\prime E}>50 \frac{D_L}{D_m^{\prime }}=0.000365 {(P{e}_m^{\prime E})}^2+0.462P{e}_m^{\prime E}+1$	13.5%
2	$\frac{D_L}{D_m^{\prime }}=1+0.114{(P{e}_m^{\prime E})}^{1.283}$	15.5%
3	$\frac{D_L}{D_m^{\prime }}=1+1.096P{e}_m^{\prime L}$	3.1%

for RPS.

Among these three correlation schemes, the third correlation with effective the Lagrangian Peclet number appeared to be superior for three reasons. It had a higher accuracy in comparison with the first and second schemes, based on an effective Eulerian Peclet number, as seen by the smaller mean relative error of Scheme 3, especially for the case of RPS. Moreover, while there was uncertainty in selecting the form of the correlation equation with effective Eulerian Pe, the form of the equation with effective Lagrangian Pe was uniform no matter what mode of particle transport was dominant (molecular diffusion or convection). There were at least two unknown parameters in the correlation equation with Eulerian Pe, and it was difficult to determine to which properties they related. On the contrary, the equation of the effective Lagrangian Pe contained only one unknown coefficient, A, as given by:

$$\frac{D_L}{D_m^{\prime }}=1+A P{e}_m^{\prime L}$$

(20)

or

$$D_L=D_m^{\prime }+A{u}_p{l}^L=D_m^{\prime }+A \frac{l^L{}^2}{{\tau }^L}.$$

(21)

Though this is a simple equation, it still can express the nature of hydrodynamic dispersion. According to this equation, the coefficient $D_L$ is the sum of two terms. The first term is the effective molecular diffusion coefficient $D_m^{\prime }$, and the second term $A l^L{}^2/{\tau }^L$ is proportional to the characteristic length scale squared over the characteristic timescale. Importantly, with this kind of form, the second term represents the coefficient of mechanical dispersion. Therefore, this correlation equation matches with the theory that hydrodynamic dispersion is a combination of molecular diffusion and convective transport. This confirmed that the Lagrangian scale is a proper approach to represent the theory of hydrodynamic dispersion. This correlation type was further confirmed with the BCC-packed bed, as shown in Figure 8, in order to conclude that it is valid for different packing types. In this equation, only the parameter A changes with various geometry structures; therefore, it is reasonable to argue that A represents the properties of porous media. The values of A for three porous media are briefly summarized in

Table 6. Summary of properties of different packing types of porous media together with the A value.

Geometry	Porosity	Permeability k (cm2)	Diameter of Spheres dp (cm)	Darcy Number Da = k/dp2	A
FCC	0.26	8.75 × 10−9	7.07 × 10−3	1.75 × 10−4	0.423
BCC	0.32	3.82 × 10−8	8.66 × 10−3	5.10 × 10−4	0.933
RPS-432 spheres	0.37	1.39 × 10−9	1.41 × 10−3	7.04 × 10−4	1.096

, together with their properties such as porosity, permeability, and Darcy number.

The linear connection between dispersion, the diffusion ratio, and the effective Lagrangian Pe in this study was only confirmed for packed beds. Further work should focus on the question of whether this finding holds for other types of real porous media found in the subsurface for environmental or oil recovery applications, such as sandstones or limestones. Moreover, the link between A and the properties of a porous medium should be further investigated.

5. Conclusions

This work examined the physics of hydrodynamic dispersion in porous media through a Lagrangian point of view. Hydrodynamic dispersion is one of those porous media-related topics that are typically associated with a lot of uncertainty and a lack of fundamental understanding, even though several applications are critically dependent on predicting dispersion in porous media. The conventional approach has been the Eulerian approach, often through macroscopic tools of analysis. Dispersion, however, is an effect that is driven by particle transport phenomena that are below the scale of consideration for macro-scale models. The present study was focused on providing a fresh approach and on identifying the relevant parameters, dimensionless numbers, and quantities. Given the currently available computational techniques that allow for in-depth calculations devoid of empiricisms, such as the use a lattice Boltzmann models and specialized particle tracking algorithms that provide paths of individual particles from a Lagrangian perspective, it is possible to analyze hydrodynamic dispersion in a physically sound framework and to probe its relationship to the effective Peclet number from both the Eulerian and Lagrangian viewpoints.

The results led to the definition of a Peclet number that is based on more naturally relevant scales rather than using Eulerian macroscopic quantities. The length scale used in defining an effective Eulerian Peclet number was the diameter of the spheres making up the porous media and was constant for a certain porous medium configuration. In contrast, the effective Lagrangian Peclet number was calculated from the Lagrangian length and time scales, which varied with molecular diffusion (i.e., the Schmidt number), the fluid velocity, and the pore geometry. This new definition has a strong impact on simplifying and unifying the correlation between the ratio of dispersion over diffusion with the Peclet number in packed beds, as seen in Equation (20). This unified correlation reveals a linear dependence on the effective Lagrangian Peclet number for different packing types, with the clear appearance of a geometrical coefficient A as the slope of the line. We conclude that the Lagrangian scales are the appropriate scales for studying hydrodynamic dispersion, as one might have predicted because both molecular diffusion effects and flow are taken into consideration