Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 846854, 10 pages

http://dx.doi.org/10.1155/2013/846854

Research Article

Modeling the CO₂ Sequestration Convection Problem Using the Lattice Boltzmann Method

Hassen M. Ouakad

Department of Mechanical Engineering, King Fahd University of Petroleum and Minerals, Dhahran 31261, Saudi Arabia

Received 10 May 2013; Revised 3 October 2013; Accepted 16 October 2013

Academic Editor: Jian Guo Zhou

Copyright © 2013 Hassen M. Ouakad. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This paper presents an investigation of the density-driven problem that rises during the CO₂ sequestration into saline aquifer. The lattice Boltzmann method (LBM) is implemented in a way to solve this mixing problem (the brine problem along with the solute transport). The CO₂-brine interface was located at the top of the considered domain. Different Rayleigh numbers were used in order to investigate this problem. When Rayleigh number is low, we got steady-state concentration contours describing a Rayleigh-Bénard type of convection. Moreover, when the Rayleigh number was selected to be big enough, we observe that the system is less stable and a convective fingering is initiated. This instability is caused by a higher density difference between the brine and the sequestrated CO₂. Note here that the turbulence is not taken into account in the study. After the onset this convective instability, the brine with a high CO₂ concentration migrates down into the porous medium. This study is based on a statistical LBM theory without assuming periodicity in any directions and without considering any type of disturbances in order to turnon the instability behavior.

1. Introduction

Carbon dioxide storage in deep saline aquifers is considered a possible option to bring greenhouse gas emissions under control. The final goal of this process is to dissolve the CO₂ into the sub-surface saline water. Injected CO₂ dissolves into formation brines from above, increasing brine density and creating an unstable hydrodynamic state favorable for natural convection. After injection, there exist sort of gravity forces that tend to redistribute the injected CO₂. This mechanism is known as a convection-advection process. The understanding of this mechanism and the associated mechanisms affecting this mixing may have an impact on the long-term sequestration process in deep saline aquifers. Several groups [1–18] used linear stability theory and high-order numerical integrations to study how this unstable convection affects the rate of dissolution of CO₂ whereas others [19–32] used the Lattice Boltzmann Method (LBM) to solve fluid flow problems in different media. Finally, laboratory studies [33–35] have confirmed qualitative and quantitative aspects of CO₂ solute convection process and found similarity with previous simulated results.

Lindeberg and Bergmo [1] presented numerical simulations of the long-term fate of CO₂ in an aquifer. They proposed that the gravitational instability of this convective mixing problem can be studied in the context of a finite domain, bounded at the top by a constant concentration boundary, using the Boussinesq approximation. Hassanzadeh et al. [2, 3] developed a 2D numerical model to study the diffusive and convection mixing of CO₂ in geologic storage. Their study revealed two important time scales involved in this process: the time to start the instability and the time to achieve ultimate dissolution. They showed that depending on the system Rayleigh number and the heterogeneity, convective mixing can greatly accelerate the dissolution of CO₂ in an aquifer.

Some studies [4–8] have also been performed investigating the linear and global stability analysis of the time-dependent density-driven convection in deep saline aquifers for long-term storage of CO₂. They investigated the onset time for convection, the preferred wavelength for the growth of convective fingers, and growth rates. Hassanzadeh et al. [9, 10] presented scaling analysis of the convective mixing of CO₂ in saline aquifers based on direct numerical simulations. They subjected the top layer to a sudden rise in CO₂ concentration from the top. They ascertained the stability of the system by integrating numerically a Galerkin-based reduced-order model. They showed results revealing increase in the onset time of convection with the transverse dispersion. Rapaka et al. [11] used the idea of nonmodal stability analysis to compute the maximum amplification of perturbations in the density-driven convection resulting from dissolution of CO₂ in brine in the underlying medium, optimized over the entire space of initial perturbations. They presented 3D spectral calculations of the problem governing equations. Their results showed good agreement compared to those obtained from spectral calculations. The same group [12] extended the analysis they adopted in [11] to the important cases of anisotropic and layered porous media with a permeability variation in the vertical direction and while using a modal analysis based on the so-called Galerkin method. They showed that, unlike the solution of the initial value problem, results obtained using nonmodal analysis used in [11] is insensitive to the choice of bottom boundary condition.

Ghesmat et al. [13, 14] studied the velocity dependent dispersion influence on the mechanism of dissolution of CO₂ into brine by fully nonlinear numerical simulations. Their investigation revealed that this dispersion increases the dissolution of CO₂ as its strength increases and that it also has an impact on the fingering pattern.

Some groups [15, 16] developed parallel simulator, based on the TOUGH2 code, for large-scale and long-term CO₂ geologic sequestration in saline aquifers. They investigated the convective mixing induced by a small increase in brine density due to dissolution of CO₂. Their numerical model involves more than 1 million grid blocks, which is currently impossible to run without using a parallel simulator [15]. In the same optic, Lu and Lichtner [17] investigated the convective instability of CO₂ by taking advantage of a 3D parallel computing using the code PFLOTRAN. The onset of this instability was resolved with high-resolution numerical simulations to investigate the rate of plume dissolution caused by fingering phenomena. The calculated rate of dissolution of CO₂ into the brine was found to be highly dependent on grid resolution. Farajzadeh et al. [18] investigated the transient density-driven natural convection problem when the top of the porous medium initially saturated with a liquid is exposed to a CO₂-rich gaseous phase. They solved this problem numerically using the mass and momentum conservation laws and diffusion of CO₂ into the liquid. Their results showed that the effect of natural convection increases with increasing Rayleigh number, and that the nonlinear behavior depends strongly on this nondimensional number.

Recently, the lattice Boltzmann method (LBM) has been proposed as an alternative numerical method for simulating fluid flows as well as physics in fluids. The method is principally successful in fluid flow applications involving nonlinearities as well as complicated domains and boundaries. Unlike the methods used in the above literature review (conventional numerical schemes based on discretizations of macroscopic continuum equations), the LBM is based on coupled microscopic models and mesoscopic kinetic equations. Some studies [19–21] used LBM to model thermal convection problems, while others [22–25] used it to model fluid flow, heat-solute transport problems [26], and solute transport in porous media [27–29].

Chen and Doolen [22] and Aidun and Clausen [23] presented a review of the use of the LBM as a parallel and efficient algorithm for simulating single-phase and multiphase fluid flows, for incorporating additional physical complexities, and for modeling complicated boundary conditions and multiphase interfaces. Both groups concluded that this method represents a mature and efficient substitute for Navier-Stokes solvers in many flow problems. Dong et al. [24] investigated the viscous fingering phenomenon of two immiscible fluids in a channel by applying the LBM. They claimed that the LBM can be viewed as a promising tool for investigating fluid behaviour and other immiscible displacement problems by providing good understanding of the mechanisms of viscous fingering. Thorne and Sukop [25] used the LBM to solve the classic solute-induced buoyancy Elder problem [26, 27]. Their work represented first steps toward using lattice Boltzmann models for application to groundwater problems (salt-water intrusion) in coastal regions. Bardsley et al. [28] adopted the Lattice Boltzmann modeling to simulate large-scale density-dependent ground water flow and heat/solute transport systems. Anwar and Sukop [29] presented LBM for modeling the flow and transport of solute in heterogeneous porous media and compared their results with available analytical solutions. Recently, Boek [30] described the ongoing development of lattice-Boltzmann (LB) computer simulation codes to study flow in porous media at the pore scale using X-ray Micro Tomography (XMT) pore space images. He reviewed the development of his codes to study flow in 2D micromodels as well as in 3D. Chen and Zhang [31] used a pore-scale LBM to simulate the density-driven convection problem in a porous medium in order to study the possibility of CO₂ sequestration in deep saline aquifers. Add to that, recent developments of the LBM for large-scale haemodynamic applications was discussed in [32].

From the aforementioned review, one can note that a robust numerical method is still needed to tackle the buoyancy-driven flow of CO₂ during its geological storage into saline aquifer. It is well-known fact that numerical scheme, based on finite-difference method for example, suffers from convergence as well as numerical diffusion problems. Add to that, method based on the TOUGH2 codes are time-consuming since they need to run more lots of grid blocks equations. In addition, many groups tend to use assumptions to be able to solve such problems by means of the previous mentioned numerical method (such as assuming periodicity in certain directions, or inject some noise signals to turn on the nonlinearity of the system).

The organization of this paper is as follows. First, a Lattice Boltzmann scheme is proposed to tackle the coupled convection-advection problem of CO₂ into brine. Then, a convergence is ascertained in order to get the suitable lattice to be used for this problem. Finally, CO₂ concentration as well as steam-function contours are shown for different Rayleigh numbers and iteration times.

2. Problem Formulation

In this section, we formulate the problem governing the flow of injected CO₂ into saline aquifer (brine), Figure 1.

Figure 1: Schematic of the CO₂-brine convective mixing process.

After injection, gravity forces will tend to accumulate the injected CO₂ at the top of the aquifer. This is modeled by assuming that the CO₂ will form a horizontal layer in top of the brine layer [8]. We consider here a 2D, homogeneous and isotropic porous layer of a finite thickness. We assume that the model is saturated with the brine. Add to that, symmetry condition at the lateral boundaries and no-flow in the bottom boundary at all times is been set. We assume also that CO₂ is injected at the top, which is assumed to be under thermodynamic equilibrium with the brine at the upper boundary interface. Therefore, we propose also to use a finite domain, bounded at the top by a constant concentration boundary and under the Boussinesq approximation as was assumed by [1]. The equations describing the Boussinesq-type of flow in a horizontal porous layer are given as follow [9]:

where is the Darcy velocity, is the concentration of the dissolved CO₂, is the permeability, is the diffusion coefficient, is the porosity, is the volumetric expansion factor, and is the viscosity.

For convenience, we introduce the following nondimensional variables:

where

and where is the Density difference between of the CO₂ saturated into brine and the fresh brine.

Substituting (4)-(5) into (1)–(3) and dropping the hats, the nondimensional equations are written as

where is the Rayleigh number which is the only nondimensional parameter present in the problem and is given by:

Using the vorticity formulation for (1) to eliminate the pressure. Hence, we take the curl of (1) and get

where is the vorticity and is the streamfunction related to the velocity components as follow:

Equation (8) is the so-called Poisson’s equation with two unknown functions and .

The boundary conditions are given as follow:

3. The Lattice Boltzmann Method

3.1. Overview

The lattice Boltzmann method is based on the use of a particle velocity distribution function, [35], which quantifies the probability to monitor a fluid particle with velocity at a lattice location and time . The subscripts “” of the distribution functions indicate the different discrete lattice directions of the considered lattice, Figure 2. In Figure 2, we show two examples of lattice Boltzmann: the one with five directions called D2Q5 and the one with nine directions called D2Q9.

(a)

(b)

Figure 2: Schematic figure showing the (a) D2Q5 and (b) D2Q9 lattices, with some nonzero vectors of the particle velocity distribution function at a node.

The fluid particle velocity distribution functions are defined for particles moving simultaneously along each lattice. The occurring velocity vectors depend on the number of sub-lattices. The fluid particles can collide with each other as they move under applied external forces.

In the LBM, the temporal evolution of the particle velocity distribution function satisfies the following equation:

where is the lattice time step. The index stands for the base vectors of the underlying lattice type. The left hand term of (5) is the advection term which represents the free propagation of the fluid particles. The right term of (5) represents the collision operator. The term represents the new distribution function after advection and rearrangement. The function represents the local equilibrium particle distribution, which depends only on the locally conserved mass and momentum density.

When considering an additional external source of momentum (e.g., body forces such as occurring in pressure gradients or gravitational fields), (5) is written as

where is the applied external force.

In the above equations, (11) and (12), the relaxation time is a parameter which characterizes the constitutive behavior of the considered fluid at a microscopic level. It is related to the macroscopic kinematic viscosity of the simulated fluid as

The equilibrium distribution for an incompressible fluid, , which approximates the Maxwell-Boltzmann equilibrium distribution up to a second order Taylor series, can be written as

where represent weight factors, is the mass-density, and is the velocity vector given as follow:

3.2. Solute Transport

The transport of solute or concentration is described with the advection-diffusion equation (see (2)), where is the concentration of the solute. The lattice Boltzmann equation can be used to describe the solute advection-diffusion process into another fluid. Hence, we introduce another probability distribution function, , for the solute concentration as follow:

where

and the relaxation time is a parameter which characterizes the solute and is related to the macroscopic diffusion coefficient as follows:

3.3. LBM Steps

The main steps in the LBM are the following:

(1)definition of the boundary conditions;(2)initialization of initial values for density, velocity, and concentration;(3)calculation of the equilibrium distribution functions with these given values;(4)propagation of the particle portions to the next neighbor: streaming step;(5)collision step (see (12) and (16));(6)calculation of the new density, velocity and concentration distributions (see (15) and (17));(7)after this last step, we need to iterate on time by increasing it by one unit and we start again with the calculation of the equilibrium distribution (Step  3).

3.4. Boundary Conditions

For the definition of the boundary conditions in LBMs, we need to translate them into the microscopic domain. In fact, we usually have some macroscopic information such asthe no-slip condition on walls for fluid flows. Therefore, in the LB algorithm, one has to translate the macroscopic information into the microscopic distribution functions and should be aware that the choice of certain boundary conditions is of primary importance since it affects quantitatively (accuracy) and qualitatively (stability) of the simulation results.

As described above in the problem formulation as well as the LBM overview and solute transport sections, we encounter two mixed problems: the fluid flow and the solute flow. For the fluid flow, a no-slip boundary condition was used at the top and bottom boundaries, and a periodic boundary condition was used at the left and right boundaries. For the solute transport, the concentration of the top boundary was fixed, and a zero concentration gradient boundary condition was employed at the left, right, and bottom boundaries. Hence, we have here two types of boundary conditions (see (10)): constant boundary (zero velocity and constant concentration at the top), and constant gradient (zero concentration gradient).

Note here that we used a full-way bounce-back boundary condition in order to implement the no-slip (zero-velocity) walls. In the full-way bounce-back boundary condition, each fluid particle distribution function is replaced by the value of the one with a velocity vector pointing in opposite direction. This identity completes the collision step, and the streaming step is executed directly after.

For the constant boundary condition type, we will describe here the case of constant concentration at the top by considering the boundary node of Figure 3 which is aligned with the direction. We notice that the distribution functions , , , and for the D2Q5 lattice and , , , , , and for the D2Q9 one are inside the boundary (known), whereas for the D2Q5 and , , and for the D2Q9 are outside the boundary (unknown). Hence, for both cases, after streaming the functions laying inside the boundary will be known functions. Hence, we need one equation for the first lattice case and three equations for the second lattice case.

Figure 3: Schematic figure showing the D2Q5 and the D2Q9 lattices at the top boundary. The grey area represents the outside of the considered system. The known particle velocity distribution functions from inside the system are (, , , and ) for the D2Q5 and (, , , , , and ) for the D2Q9.

For the D2Q5 case, we use (17), that is,

For the D2Q9 case, we use also (17) along with the bounce-back conditions in the longitudinal directions of the lattice [32, 36], that is,

Finally, for the zero concentration gradient (Neumann-type boundary condition), we use the following identity [37] to satisfy this condition in the rest of the boundaries (dashed boundaries in Figure 3):

where is the boundary location.

3.5. LBM Nondimensional Parameters

We should note here that there are two nondimensional parameters that are controlling the studied density-driven flow and solute transport problem in the LBM process. These parameters are the Rayleigh and Prandtl numbers defined respectively as follows:

where is the lattice height given in lattice unit (LU).

3.6. Porous Media Step

There exist in the literature several groups [25, 38, 39] that attempted to model the porous media through the LBM. One of these groups [39] stated that their model is much more efficient and can handle this problem without suffering from any kind of convergence issues. In fact, to model the porous media in the LBM, Walsh et al. [39] proposed to add a new step along with the streaming and the collision steps. We note here by the given particle velocity distribution function and by the distribution function calculated after the collision step. The new calculated distribution function that account for the porous media step is given as follow:

where is a constant for a homogenous medium and called the volumetric solids fraction parameter used to damp the evolution of each particle momentum in the considered porous media.

4. Results

4.1. Rayleigh-Benard Type of Convection

We start here our simulations by considering a nonporous media ( = 0). Figure 4 show a simulation of a Rayleigh-Bénard type of convection. Here, we set = 4000, and = 1 and we present the steady state concentration contours as well as the stream-function and velocity distributions. It can be seen that we ended-up having a Rayleigh-Bénard type of contours and this results is in perfect agreement with what was presented in [31].

(a)

(b)

(c)

(d)

Figure 4: Lattice Boltzmann simulated velocity, concentration, and stream-function contours with = 4000, = 1, and = 0 after 20000 time iteration.

4.2. Effect of Rayleigh Number

In this subsection, numerical simulations are performed to find the onset of convection in large-scale fields with some given properties. Numerical simulations were performed by changing the molecular diffusivity, model thickness, and porous medium permeability, resulting in a wide range of Rayleigh numbers. Results for the onset of convection are illustrated in Figures 5, 6, and 7 as a function of Rayleigh number. We have noticed that with the continuous increase of the Rayleigh number , we noticed that the simulated concentration contours became more and more unpredictable until we are get a convective-mixing type of behavior. The concentration contours start to contract and we saw initiation of concentration fingers that migrates downward in a convective-mixing way, Figures 5–7. The theoretical stability analyses performed by several researchers [4, 8, 10] suggest that Rayleigh should be larger than certain critical values, which depends on aquifer’s characteristics. Hence, for CO₂ storage applications, the geometrical, kinematical, and material properties should be well selected for a healthy injection process.

(a)

(b)

(c)

(d)

Figure 5: Lattice Boltzmann simulated velocity, concentration and stream-function contours with = 2 * 10⁵, = 1, and = 0 after 2000 time iteration.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Lattice Boltzmann simulated concentration contours with = 2 * 10⁴, = 1, and = 0 for various Rayleigh number.

Figure 7: Lattice Boltzmann simulated concentration contours with = 2 * 10⁵, = 1, and = 0 for various time iteration.

4.3. Effect of the Porous Media Step

Finally, we include here in the above described LBM model, the porous media step through the volumetric solids fraction parameter . The main approach used is the so-called GLB model that consider the partial-bounce back term in the LBM model and simply take a form of   *  , where represents the net fluid particles to be bounced back due to the existence of solid phase. Four GLB models, which have been proposed in many previously published works [25, 38, 39], the models take the same form of the partial-bounce back term, and they treat a solid boundary node as a special case of a gray node where = 1, that is, partial-bounce back term in the LBM reduces to the standard bounce-back but with a delay of one time step. As a result, no additional treatment is required for treating the boundary condition. Solving simultaneously the LBM equations, (14) and (16), along with the porous step equation, (24), we notice that the porous media step damp the whole mixing system and we end-up changing the whole mixing behavior from transient with concentration contours with form of fingers, Figure 8(a), to steady-state parallel concentration contours, Figure 8(b). We are currently still investigating the effect of inclusion of this step in the whole process of the LBM.

(a)

(b)

Figure 8: Lattice Boltzmann simulated concentration contours with = 1 * 10⁵, = 1, and (a) = 0, and (b) = 0.9 after 5000 time iteration.

5. Conclusions

In this paper, an investigation into the density-driven problem that rises during the CO₂ sequestration into saline aquifer was presented. The Boussinesq-type of flow has been solved using a Lattice Boltzmann method. We have shown possible scenario of a Rayleigh-Bénard type of convection for low Rayleigh number. Then, a convective-mixing loop was obtained for high Rayleigh number. To account for the porous media effect, we proposed to use a new step in the LBM process. The results indicated that the inclusion of this step may damp the whole mixing system and may affect the whole convective-mixing process.

References

E. Lindeberg and P. Bergmo, “The long-term fate of CO₂ injected into an aquifer,” in Greenhouse Gas Control Technologies, J. Gale and Y. Kaya, Eds., vol. 1, pp. 489–495, Elsevier, 2003.

H. Hassanzadeh, M. Pooladi-Darvish, and D. W. Keith, “Modelling of convective mixing in CO₂ storage,” Journal of Canadian Petroleum Technology, vol. 44, no. 10, pp. 43–50, 2005.

H. Hassanzadeh, Mathematical modeling of convective mixing in porous media for geological carbon dioxide storage [Ph. D. dissertation], University of Calgary, Alberta, Canada, 2006.

X. Xu, S. Chen, and D. Zhang, “Convective stability analysis of the long-term storage of carbon dioxide in deep saline aquifers,” Advances in Water Resources, vol. 29, no. 3, pp. 397–407, 2006.

J. Ennis-King and L. Paterson, “Rate of dissolution due to convective mixing in the underground storage of carbon dioxide,” in Greenhouse Gas Control Technologies, J. Gale and Y. Kaya, Eds., vol. 1, pp. 507–510, Elsevier, Amsterdam, The Netherlands, 2003.

J. Ennis-King and L. Paterson, “Role of convective mixing in the long-term storage of carbon dioxide in deep saline formations,” in Proceedings of the Society of Petroleum Engineers Annual Fall Technical Conference and Exhibition, SPE-84344, Denver, Colo, USA, October 2003.

M. A. Hesse, H. A. Tchelepi, and F. M. Orr Jr, “Natural convection during aquifer CO₂ storage,” in Proceedings of the 8th International Conference on Greenhouse Gas Control Technologies, Trondheim, Norway, June 2006.

A. Riaz, M. Hesse, H. A. Tchelepi, and F. M. Orr Jr., “Onset of convection in a gravitationally unstable diffusive boundary layer in porous media,” Journal of Fluid Mechanics, vol. 548, pp. 87–111, 2006.

H. Hassanzadeh, M. Pooladi-Darvish, and D. W. Keith, “Scaling behavior of convective mixing, with application to geological storage of CO₂,” AIChE Journal, vol. 53, no. 5, pp. 1121–1131, 2007.

H. Hassanzadeh, M. Pooladi-Darvish, and D. W. Keith, “The effect of natural flow of aquifers and associated dispersion on the onset of buoyancy-driven convection in a saturated porous medium,” AIChE Journal, vol. 55, no. 2, pp. 475–485, 2009.

S. Rapaka, S. Chen, R. J. Pawar, P. H. Stauffer, and D. Zhang, “Non-modal growth of perturbations in density-driven convection in porous media,” Journal of Fluid Mechanics, vol. 609, pp. 285–303, 2008.

S. Rapaka, R. J. Pawar, P. H. Stauffer, D. Zhang, and S. Chen, “Onset of convection over a transient base-state in anisotropic and layered porous media,” Journal of Fluid Mechanics, vol. 641, pp. 227–244, 2009.

K. Ghesmat, H. Hassanzadeh, and J. Abedi, “The impact of geochemistry on convective mixing in a gravitationally unstable diffusive boundary layer in porous media: CO₂ storage in saline aquifers,” Journal of Fluid Mechanics, vol. 673, pp. 480–512, 2011.

K. Ghesmat, H. Hassanzadeh, and J. Abedi, “The effect of anisotropic dispersion on the convective mixing in long-term CO₂ storage in saline aquifers,” AIChE Journal, vol. 57, no. 3, pp. 561–570, 2011.

K. Zhang, C. Doughty, Y.-S. Wu, and K. Pruess, “Efficient parallel simulation of CO₂ geologic sequestration in saline aquifers,” in Reservoir Simulation Symposium, Proceedings of SPE, pp. 26–28, Houston, Tex, USA, February 2007.

K. Pruess and K. Zhang, “Numerical modeling studies of the dissolution-diffusion-convection process during CO₂ storage in saline aquifers,” Lawrence Berkeley Laboratory Report LBNL-1243E, University of California, Berkeley, Calif, USA, 2008.

C. Lu and P. C. Lichtner, “High resolution numerical investigation on the effect of convective instability on long term CO₂ storage in saline aquifers,” Journal of Physics, vol. 78, no. 1, Article ID 012042, 2007.

R. Farajzadeh, H. Salimi, P. L. J. Zitha, and H. Bruining, “Numerical simulation of density-driven natural convection in porous media with application for CO₂ injection projects,” International Journal of Heat and Mass Transfer, vol. 50, no. 25-26, pp. 5054–5064, 2007.

P.-H. Kao and R.-J. Yang, “Simulating oscillatory flows in Rayleigh-Bénard convection using the lattice Boltzmann method,” International Journal of Heat and Mass Transfer, vol. 50, no. 17-18, pp. 3315–3328, 2007.

Y. Lizhong, Y. Junqi, and W. Yafei, “Simulation of gravity currents using the thermal lattice Boltzmann method,” International Journal for Numerical Methods in Fluids, vol. 66, no. 9, pp. 1174–1182, 2011.

N. A. C. Sidik and A. R. M. Rosdzimin, “Cubic interpolated pseudo particle (CIP)—thermal BGK lattice Boltzmann numerical scheme for solving incompressible thermal fluid flow problem,” Malaysian Journal of Mathematical Sciences, vol. 3, no. 2, pp. 183–202, 2009.

S. Chen and G. D. Doolen, “Lattice boltzmann method for fluid flows,” Annual Review of Fluid Mechanics, vol. 30, pp. 329–364, 1998.

C. K. Aidun and J. R. Clausen, “Lattice-boltzmann method for complex flows,” Annual Review of Fluid Mechanics, vol. 42, pp. 439–472, 2010.

B. Dong, Y. Y. Yan, W. Li, and Y. Song, “Lattice Boltzmann simulation of viscous fingering phenomenon of immiscible fluids displacement in a channel,” Computers and Fluids, vol. 39, no. 5, pp. 768–779, 2010.

D. T. Thorne and M. C. Sukop, “Lattice Boltzmann model for the Elder problem, developments in water science,” in Proceedings of the 5th International Conference on Computational Methods in Water Resources, vol. 55, pp. 1549–1557, 2004.

J. W. Elder, “Transient convection in a porous medium,” Journal of Fluid Mechanics, vol. 27, pp. 609–623, 1967.

J. W. Elder, “Steady free convection in a porous medium heated from below,” Journal of Fluid Mechanics, vol. 27, pp. 29–48, 1967.

K. J. Bardsley, S. Anwar, and M. C. Sukop, “Simultaneous heat and solute transport modeling of ground-water with lattice Boltzmann methods,” in Proceedings of the 16th International Conference on Computational Methods in Water Resources (CMWR '06), Copenhagen, Denmark, June 2006.

S. Anwar and M. C. Sukop, “Lattice boltzmann models for flow and transport in saturated karst,” Ground Water, vol. 47, no. 3, pp. 401–413, 2009.

E. S. Boek, “Pore scale simulation of flow in porous media using lattice-Boltzmann computer simulations,” in Annual Technical Conference and Exhibition, Proceedings of SPE, no. 135506, Florence, Italy, 2010.

C. Chen and D. Zhang, “Pore-scale simulation of density-driven convection in fractured porous media during geological CO₂ sequestration,” Water Resources Research, vol. 46, no. 11, Article ID W11527, 2010.

P. Giuseppe, I. Hallidayb, S. Melchionnac, T. J. Spencerb, and S. Succi, “Lattice Boltzmann method as a computational framework for multiscale haemodynamics,” Mathematical and Computer Modelling of Dynamical Systems, pp. 1–21, 2013.

R. Farajzadeh, H. Salimi, P. L. J. Zitha, and H. Bruining, “Numerical simulation of density-driven natural convection in porous media with application for CO₂ injection projects,” International Journal of Heat and Mass Transfer, vol. 50, no. 25-26, pp. 5054–5064, 2007.

T. J. Kneafsey and K. Pruess, “Laboratory flow experiments for visualizing carbon dioxide-induced, density-driven brine convection,” Transport in Porous Media, vol. 82, no. 1, pp. 123–139, 2010.

J. A. Neufeld, M. A. Hesse, A. Riaz, M. A. Hallworth, H. A. Tchelepi, and H. E. Huppert, “Convective dissolution of carbon dioxide in saline aquifers,” Geophysical Research Letters, vol. 37, no. 22, Article ID L22404, 2010.

M. C. Sukop and D. T. Thorne Jr, Lattice Boltzmann Modeling: an Introduction for Geoscientists and Engineers, Springer, Berlin, Germany, 2006.

O. Malaspinas, “How to impose a Neumann boundary condition with the lattice Boltzmann method,” 2008, http://www.lbmethod.org

O. Dardis and J. McCloskey, “Lattice Boltzmann scheme with real numbered solid density for the simulation of flow in porous media,” Physical Review E, vol. 57, no. 4, pp. 4834–4837, 1998.

S. D. C. Walsh, H. Burwinkle, and M. O. Saar, “A new partial-bounceback lattice-Boltzmann method for fluid flow through heterogeneous media,” Computers and Geosciences, vol. 35, no. 6, pp. 1186–1193, 2009.