HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via ISI Web of Science (1) Citing Articles via Google Scholar Google Scholar Articles by Jiménez-Hornero, F. J. Articles by Laguna, A. Search for Related Content PubMed Articles by Jiménez-Hornero, F. J. Articles by Laguna, A. GeoRef GeoRef Citation Agricola Articles by Jiménez-Hornero, F. J. Articles by Laguna, A. Related Collections Numerical Solutions Solute Transport Models Dispersion

SPECIAL SECTION: ZNS'03 VADOSE ZONE RESEARCH

Simulation of Tracer Dispersion in Porous Media Using Lattice Boltzmann and Random Walk Models

^a Dep. of Agronomy, University of Córdoba, P.O. Box 3048, 14080 Córdoba, Spain
^b Dep. of Agronomy, University of Córdoba and Institute of Sustainable Agriculture, CSIC, P.O. Box 3048, 14080 Córdoba, Spain
^c Dep. of Applied Physics, University of Córdoba, P.O. Box 3048, 14080 Córdoba, Spain
^* Corresponding author (jihof{at}arrakis.es)

Received 16 June 2004.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS CONCLUSIONS REFERENCES

 
Pore-scale flow and transport processes are generally difficult to simulate using conventional models. Two state-of-the-art approaches are adopted in this study of solute transport in porous media. The first approach is a lattice Boltzmann model with the Bhatnagar, Gross, and Krook simplification (BGK) to determine both the advection and diffusion components of the transport process, while the second approach uses the BGK model to calculate the velocity field in conjunction with a random walk (RW) model to characterize diffusion (BGK/RW). Both approaches yield similar results. The BGK/RW model is an attractive alternative to the BGK model when it is necessary to simulate small values of the diffusion coefficient so as to avoid instabilities in the numerical solution. Nevertheless, the BGK/RW model is less accurate than the BGK model alone when compared with the analytic solution of the well-known Taylor–Aris dispersion model.

Abbreviations: BGK, Bhatnagar, Gross, and Krook model • BGK/RW, Bhatnagar, Gross, and Krook model, with a random walk model • LBM, lattice Boltzmann model • l.u., lattice units • RW, random walk model

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS CONCLUSIONS REFERENCES

 
ONE OF THE MOST important challenges in the soil and hydrologic sciences is accurate simulation of solute transport in porous media. The wide variety of soil aggregate types and sizes, and, consequently, of pores between them, hinders a rigorous analysis of solute transport in the subsurface. While many simplified expressions exist for macroscopic transport based on simplifying assumptions, these solutions generally fail to account for the geometric complexities of the porous medium, and hence provide only first-order approximations to the solute transport process. In the past, incorporating these complexities made it necessary to use elaborate numerical models that are extremely difficult to implement.

The purpose of this work is to introduce a mesoscopic lattice BGK model that has been successfully applied to similar problems in other fields. Models of this type are derived from the lattice Boltzmann model (LBM), a numerical scheme recently developed for studying fluid dynamics problems (Chen and Doolen, 1998; Wolf-Gladrow, 2000, Chapter 5). The LBM model describes phenomena using a mesoscopic scale in which the real system is transformed into a synthetic medium consisting of simple particles interacting with each other according to relatively simple rules. Most of the reported versions of the LBM are extensions of those proposed by Chen et al. (1992) and Qian et al. (1992) on the basis of a simplification of the relaxation time as suggested by Bhatnagar, Gross, and Krook (Bhatnagar et al., 1954).

Among many other applications, BGK models have been adopted in the analysis of dispersion phenomena (Flekkøy, 1993; Flekkøy et al., 1995; Zhang et al., 2002). Application of the BGK model to solute transport in porous media allows a more accurate description of the two main processes, advection and dispersion.

Nevertheless, the use of the BGK model in solute dispersion problems is limited by the difficult simulation of prescribed values of kinematic viscosity and diffusion coefficient, especially when the values of these parameters imply relaxation times near 0.5 that produce numerical instabilities, as will be explained in the Materials and Methods section. For problems of solute transport in porous media, simulations of kinematic viscosities do not pose any problem since the corresponding relaxation times are >0.5. However, the simulated diffusion coefficients are usually small and the corresponding relaxation times are near 0.5.

One solution to this problem is the introduction of the RW model for simulating the diffusive component of the transport process. The RW model has been used previously for the description of solute transport by Ackerer (1988), Kinzelbach (1988), Dagan and Neuman (1997), and Bodin et al. (2003), among others.

In this study a combination of the BGK and RW models (referred to as the BGK/RW model) is proposed as an alternative method to the BGK model when this approach produces numerical instabilities for situations with small values of the diffusion coefficient (high Péclet numbers). The combined BGK/RW hence uses the BGK model for the advective component of the solute transport process, and the RW model for the dispersive part, as indicated by Maier et al. (1998)(2000). The purpose of this work is to compare both methods for simulating solute transport processes in porous media.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS CONCLUSIONS REFERENCES

 
The Lattice BGK Model
The BGK model is an extension of the LBM often used for solving the Navier–Stokes equations (Rothman and Zaleski, 1997, Chapter 6.3; Chen and Doolen, 1998; Chopard and Droz, 1998, Chapter 3.5; Wolf-Gladrow, 2000, Chapter 4.2), among other applications. The model can also be used for solution of the diffusion equation (Flekkøy, 1993). The BGK model consists of an ensemble of particles interacting with each other while conserving mass and momentum (Rothman and Zaleski, 1997, Chapter 6.1; Chopard and Droz, 1998, Chapter 3.5). The domain of the problem is a regular grid whose nodes are linked to their neighbors through a vicinity relationship, chosen according to the complexity of the process. A good choice for two-dimensional advection processes is an array of eight nodes surrounding a central node, (the d2q9 model), while less nodes are required for diffusion (e.g., four neighbors in the d2q4 model). Figure 1 shows both types of vicinity models.

View larger version (8K): [in this window] [in a new window]   Fig. 1. Neighborhood relationships used in the BGK model for the description of advection (d2q9), and diffusion (d2q4), where d is the dimension number and q the number of adopted neighbors.

[1]

where w_i = (r/t)·c_j, is the velocity of a particle in link i of node r node at time t, r is the length of the lattice spacing, t is the time elapsed during one time-step, c_i are the coordinates of the neighbor node linked to through link i, and q is the number of neighbors being considered.

The fundamental mesoscopic equation of the BGK model is (Qian et al., 1992; Chen et al., 1992):

[2]

where f⁽⁰⁾ is the local equilibrium term and is the relaxation time, using concepts proposed by Bhatnagar et al. (1954). The general form of f_i⁽⁰⁾ in our work is

[3]

where the Einstein summation convention has been adopted, and where k_ß is the component of any vector k in dimension ß. The t_p are weighing factors for the rest particles, t, and those moving along the vertical and horizontal, t₁, and in the diagonal directions, t₂. Finally, c_s is a parameter (known as the velocity of sound), which depends on the chosen neighborhood relationship. Table 1 gives parameters values for the d2q4 and d2q9 lattices (Succi, 2001, Chapter 5.2).

The relaxation time facilitates simple identification of two parameters: the kinematic viscosity, , (and consequently the Reynolds number, Re) and the diffusion coefficient D, which determines the Péclet number Pe:

[4]

where and _D are the relaxation times that adjust the kinematic viscosity and the diffusion coefficient, respectively, and c_s and c_sD are the velocities of sound in the selected neighborhood models for the description of advection and diffusion.

The stability of the proposed BGK model depends on the value of the two relaxation times. If one of these is <0.5, the kinematic viscosity and the diffusion coefficient would be negative. Since solute transport in porous media has a low Re value, the advective component does not experience any stability problems. Nevertheless, the diffusion coefficient is frequently small, thus leading to a high Pe value that precludes the simulation of the diffusive component of solute transport with the BGK model.

From Mesoscopic to Macroscopic Scale
Table 2 shows the conversion rules between variables in the BGK model and their physical entities (Succi, 2001, Appendix D). The scale factors r and t are expressed in SI units.

[5]

where x(t + t) is the position of the particle, initially at x(t). The advective component, v[x(t)], is the fluid velocity determined from the velocity field using the BGK model. The diffusive component, x_D(Dt), corresponds to the vector whose module is constant and equal to d, where d is the number of dimensions. The random direction of x_D(Dt) follows a Gaussian probability distribution function with a mean value of zero and a variance of 2Dt. For the two-dimensional case, d = 2,

[6]

because of the fact that the solution of the diffusion equation,

[7]

in a d-dimensional space is the Gaussian probability distribution function

[8]

with variance ² = 2Dt.

If the tracer particles do not interact between each other, and if they are suspended in the homogeneous carrier fluid, the trajectories described are random walks that occur in stages that are themselves random events. In this case it is possible to use the central limit theorem that permits the estimation of the average displacement as

[9]

and the square average displacement

[10]

where t₀ is the initial time and expresses the average operator.

When both advection and diffusion processes are considered together, the dispersion coefficient D* replaces the diffusion coefficient D. The dispersion coefficient in the x direction and time t is computed as

[11]

where

[12]

with

[13]

where x(t) is the tracer particle coordinate in the x direction. If the derivative d_x²/dt is constant, Fick's Law is a valid model for the description of the dispersion process.

The distance between two consecutive nodes, r, is the same distance as used in the BGK model for the velocity field. The time step must be such that it minimizes the number of displacements of the particle while keeping enough precision in the simulation. As a general rule, its value may be chosen from (Maier et al., 2000):

[14]

where v_max is the maximum value of the velocity field.

An important aspect of simulation is the interaction of solute with the pore wall boundaries. Although the bounce-back rule may be applied, in many cases this rule may cause either an accumulation or a reduction in the number of solute particles near the boundaries, depending on the t values. Salles et al. (1993) suggested a simple rule to avoid these problems that has been adopted here. According to the Salles et al. (1993) rule, if one particle finds a solid boundary at any site, the particle stops and remains there as a starting point for the next time step with a random direction.

The larger the particle number, the better the precision of the method, as demonstrated by Maier et al. (1998). Nevertheless, the computational effort increases with the number of solute particles.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS CONCLUSIONS REFERENCES

 
The porous medium used in this work is the random fractal lattice developed by Rappoldt and Crawford (1999), in which a cell acts as a soil matrix with a probability p = 0.7. The porosity of the medium is given by

[15]

where k is the number of recursion levels (in our case k = 4). For r = t = 1, the length and width of the problem domain are nx = ny = 156 lattice units (l.u.). The first 50 columns were taken free of obstacles to stabilize flow before entering the porous medium (Fig. 2 and 3) .

View larger version (40K): [in this window] [in a new window]   Fig. 2. Dispersion of a passive tracer in a medium with porosity = 0.76, simulated with the BGK model. The flow was forced from left to right with the parameters Re = 39.7 and Pe = 110.2. The dark and white hues correspond to regions with higher and lower tracer concentration, respectively. The concentration values are normalized with respect to their initial value. The numbers in the figure corresponds to the computed time steps.

View larger version (43K): [in this window] [in a new window]   Fig. 3. Dispersion of a passive tracer in a medium with porosity = 0.76, simulated with the BGK/RW model. The flow was forced from left to right with the parameters Re = 39.7 and Pe = 110.2. The dark and white hues correspond to regions with higher and lower tracer concentration, respectively. The concentration values are normalized with respect to their initial value. The numbers in the figure corresponds to the computed time steps.

Figure 2 shows the evolution of the tracer concentration normalized with respect to its initial concentration. Dark shaded and clear areas correspond to regions with higher and lower tracer concentrations, respectively. The corresponding number of time steps is given also. Notice that the tracer remains for a longer period in low-permeability areas of the medium where flow has difficulty passing.

The BGK/RW Model
In this case the dispersion of the solute was simulated with 9000 particles initially located between Columns 28 and 32, in the same porous medium used as before and with the same velocity field and diffusion coefficient. Zero-gradient condition was applied along the outlet (right boundary), and a periodic condition was used along the upper and lower boundaries of the domain. Figure 3 shows the temporal evolution of the normalized tracer concentration.

While Fig. 2 and 3 are similar, the latter shows that the tracer remains for a longer time in regions where dispersion is reduced due to the presence of smaller pores.

Comparison between the Two Models
The breakthrough curves for both models were computed as the amount of tracer passing through the outlet of the domain, L = nx, normalized with respect to the initial concentration and the normalized time T* = tU/L, where t is time and U the mean flow velocity. As Fig. 4 indicates, the breakthrough curve of the BGK/RW model lags behind the corresponding curve obtained in the BGK model.

View larger version (12K): [in this window] [in a new window]   Fig. 4. Breakthrough curves computed with BGK (full line) and the BGK/RW models (dashed line) for simulating the dispersion of a passive tracer in a medium with porosity = 0.76 using the parameters Re = 39.7 and Pe = 110.2.

[16]

The simulations were performed with a separation between plates of ny = 30 l.u. and a flow velocity of 0.04 l.u./time step. As before, r = t = 1. The kinematic viscosity of the carrier fluid was 0.25 l.u.²/time step ( = 1.25) and the diffusion coefficient 0.25 l.u.²/time step (_D = 1 for the BGK model). The corresponding characteristic time to reach the Taylor–Aris dispersion condition was t_c = 3600 time steps. To ensure this condition, the length of the problem domain was nx = 1200 l.u., being the number of time steps of the simulations larger than 8 times t_c. The Reynolds and Péclet numbers were computed using the separation between plates as the characteristic length obtaining for both the value of 4.8. Figure 5 shows the resulting plots of the variance (²) vs. the time (t) for both BGK and BGK/RW. The D*/D ratios, calculated with Eq. [11], were compared later with the theoretical value of 1.37753 obtained from Eq. [16]. The BGK model yielded a value of 1.40364 whereas the BGK/RW gave 1.41026, with errors of 1.860 and 2.376%, respectively.

View larger version (22K): [in this window] [in a new window]   Fig. 5. Plots of variance (2) vs. time (t) obtained with the BGK (solid line) and BGK/RW models (dashed line) for Taylor–Aris dispersion during flow with Re = Pe = 4.8. The last branch of each curve may be fitted to the indicated values for d2/dt.

[17]

where x is the distance at which the concentration of tracer is computed, and erfc is the error function. Results are shown in Fig. 6 and 7 for the BGK and BGK/RW approaches, respectively. The goodness of fit of the theoretical to the experimental breakthrough curves was estimated with three methods (Legates and McCabe, 1999): (i) the coefficient of determination R² given by

[18]

(ii) the efficiency or Nash–Sutcliffe efficiency E given by

[19]

(iii) and the agreement index d given by

[20]

View larger version (11K): [in this window] [in a new window]   Fig. 6. Breakthrough curve computed with the BGK model (solid line), for Taylor–Aris dispersion of a passive tracer during flow with the Re = Pe = 4.8, and comparison with the theoretical solution (dashed line).

View larger version (11K): [in this window] [in a new window]   Fig. 7. Breakthrough curve computed with the BGK/RW model (solid line), for Taylor–Aris dispersion of a passive tracer during flow with the Re = Pe = 4.8, and comparison with the theoretical solution (dashed line).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS CONCLUSIONS REFERENCES

When it is necessary to simulate values of the diffusion coefficient that correspond to relaxation times near 0.5, the BGK model becomes numerically unstable. The Péclet number for this type of flow is relatively high. According to the expression Pe = UL/D, there are two ways to obtain high values for Pe without modifying the diffusion coefficient: (i) increasing the velocity or (ii) the characteristic length of the flow. Both options are not advisable solutions in the case of the BGK model since U has to be lower than the velocity of sound to ensure numerical stability of the model, whereas a larger L implies an increase in the computational effort.

The use of the BGK/RW model is a feasible alternative for overcoming the BGK drawback before mentioned. For the BGK/RW approach, the main constraint is the number of particles to use to obtain the desired accuracy in the results without excessively increasing the run-time of the simulations.

Implementation of the proposed models is not difficult. The adaptation of the d2q4 neighborhood to the d2q9 code is almost trivial for the BGK model, although parallel programming techniques must be used to obtain reasonable run-times with this model. In the case of the BGK/RW approximation it is necessary to write two different codes for the BGK and RW models.

	ACKNOWLEDGMENTS

 
The support of the Spanish CICYT, through Project CAO1-1-C4-1, is thankfully acknowledged by the authors. F. J. Jiménez-Hornero is supported by Consejeria de Innovacion, Ciencia y Empresa de la Junta de Andalucia (Ayudas para facilitar el Retorno de Investigadores a Centros de Investigacion y Universidades de Andalucia).

	REFERENCES

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS CONCLUSIONS REFERENCES

 

• Ackerer, Ph. 1988. Random-walk method to simulate pollutant transport in alluvial aquifers or fracture rocks. p. 475–486. In E. Custodio et al. (ed.) Groundwater flow and quality modelling. NATO ASI Ser., Ser. C Math and Phys. Sci., 224:475–486, D. Reidel. Norwell, MA.
• Aris, R. 1956. On the dispersion of a solute in a fluid flowing through a tube. Proc. R. Soc. A 235:67–77.
• Bhatnagar, P., E.P. Gross, and M.K. Krook. 1954. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. B 94:511–525.[CrossRef]
• Bodin, J., G. Porel, and F. Delay. 2003. Simulation of solute transport in discrete fracture networks using the time domain random walk method. Earth Planet. Sci. Lett. 208:297–304.[CrossRef]
• Chen, S., and G. Doolen. 1998. Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech. 30:329–364.[CrossRef][ISI]
• Chen, S., Z. Wang, X. Shan, and G. Doolen. 1992. Lattice Boltzmann computational fluid dynamics in three dimensions. J. Stat. Phys. 68:379–400.[CrossRef]
• Chopard, B., and M. Droz. 1998. Cellular automata modeling of physical systems. Cambridge University Press, Cambridge, UK.
• Dagan, G. 1989. Flow and transport in porous formations. Springer-Verlag, New York.
• Dagan, G., and S.P. Neuman (ed.) 1997. Subsurface flow and transport. A stochastic approach. Cambridge University Press, New York.
• Flekkøy, E.G. 1993. Lattice Bhatnagar–Gross–Krook models for miscible fluids. Phys. Rev. E 47:4247–4257.[CrossRef]
• Flekkøy, E.G., U. Oxaal, J. Feder, and T. Jøssang. 1995. Hydrodynamic dispersion at stagnation points: Simulations and experiments. Phys. Rev. E 52:4952–4962.[CrossRef]
• Kinzelbach, W. 1988. The random walk method in pollutant transport simulation. p. 227–246. In E. Custodio et al. (ed.) Groundwater flow and quality modelling. NATO ASI Ser., Ser. C Math and Phys. Sci., 224:227–246. D. Reidel. Norwell, MA.
• Kutilek, M., and D.R. Nielsen. 1994. Soil hydrology. Catena Verlag, Cremlingen-Destedt, Germany.
• Legates, D.R., and G.J. McCabe. 1999. Evaluating the use of "goodness-of-fit" measures in hydrologic and hydroclimatic model validation. Water Resour. Res. 1:233–241.
• Maier, R.S., D.M. Kroll, R.S. Bernard, S.E. Howington, J.F. Peters, and H.T. Davis. 2000. Pore-scale simulation of dispersion. Phys. Fluids 12:2065–2079.[CrossRef]
• Maier, R.S., D.M. Kroll, H.T. Davis, and R.S. Bernard. 1998. Pore-scale flow and dispersion. Int. J. Modern Phys. C 9:1523–1533.[CrossRef]
• Perea-Reeves, S.J., and H.W. Stockman. 1997. A lattice-gas study of dispersion in alveolated channels. Chem. Eng. Sci. 19:3277–3286.
• Qian, Y.H., D. D‘Humieres, and P. Lallemand. 1992. Lattice BGK models for Navier-Stokes equation. Europhys. Lett. 17:479–484.
• Rappoldt, C., and J.W. Crawford. 1999. The distribution of anoxic volume in a fractal model of soil. Geoderma 88:329–347.
• Rothman, D.H., and S. Zaleski. 1997. Lattice-gas cellular automata. Simple models of complex hydrodynamics. Cambridge University Press, Cambridge, UK.
• Salles, J., J.F. Thovert, R. Delannay, L. Prevors, J.L. Auriault, and P.M. Adler. 1993. Taylor dispersion in porous media. Determination of the dispersion tensor. Phys. Fluids A 5:2348–2376.[CrossRef]
• Stockman, H.W. 1999. A 3D lattice Boltzmann code for modeling flow and multi-component dispersion. Tech. Rep. SAND99-0162. Sandia Natl. Lab., Alburquerque, NM.
• Succi, S. 2001. The lattice Boltzmann equation for fluid dynamics and beyond. Numerical mathematics and scientific computation. Oxford University Press, Oxford, UK.
• Taylor, G. 1953. Dispersion of soluble matter in solvent flowing slowly through a tube. Proc. R. Soc. A 219:186–203.
• van Genuchten, M.Th., and W.J. Alves. 1982. Analytical solutions of the one dimensional convective dispersive solute transport equation. USDA Tech. Bull. 1661. U.S. Gov. Print. Office, Washington, DC.
• Wolf-Gladrow, D.A. 2000. Lattice-gas automata and lattice Boltzmann models. Lecture Notes in Mathematics. Springer-Verlag, Berlin.
• Zhang, X., A.G. Bengough, J.W. Crawford, and I.M. Young. 2002. A lattice BGK model for advection and anisotropic dispersion equation. Adv. Water Resour. 25:1–8.

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via ISI Web of Science (1) Citing Articles via Google Scholar Google Scholar Articles by Jiménez-Hornero, F. J. Articles by Laguna, A. Search for Related Content PubMed Articles by Jiménez-Hornero, F. J. Articles by Laguna, A. GeoRef GeoRef Citation Agricola Articles by Jiménez-Hornero, F. J. Articles by Laguna, A. Related Collections Numerical Solutions Solute Transport Models Dispersion