#### Abstract

This paper presents an investigation of the density-driven problem that rises during the CO_{2} 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_{2}-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_{2}. 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_{2} 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_{2} into the sub-surface saline water. Injected CO_{2} 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_{2}. 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_{2} 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_{2} solute convection process and found similarity with previous simulated results.

Lindeberg and Bergmo [1] presented numerical simulations of the long-term fate of CO_{2} 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_{2} 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_{2} 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_{2}. 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_{2} in saline aquifers based on direct numerical simulations. They subjected the top layer to a sudden rise in CO_{2} 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_{2} 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_{2} into brine by fully nonlinear numerical simulations. Their investigation revealed that this dispersion increases the dissolution of CO_{2} 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_{2} geologic sequestration in saline aquifers. They investigated the convective mixing induced by a small increase in brine density due to dissolution of CO_{2}. 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_{2} 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_{2} 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_{2}-rich gaseous phase. They solved this problem numerically using the mass and momentum conservation laws and diffusion of CO_{2} 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_{2} 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_{2} 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_{2} into brine. Then, a convergence is ascertained in order to get the suitable lattice to be used for this problem. Finally, CO_{2} 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_{2} into saline aquifer (brine), Figure 1.

After injection, gravity forces will tend to accumulate the injected CO_{2} at the top of the aquifer. This is modeled by assuming that the CO_{2} 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_{2} 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_{2}, 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_{2} 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)**

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.

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)**

##### 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_{2} storage applications, the geometrical, kinematical, and material properties should be well selected for a healthy injection process.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

##### 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)**

#### 5. Conclusions

In this paper, an investigation into the density-driven problem that rises during the CO_{2} 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.