Pore-scale study based on lattice Boltzmann method of density driven natural con

Abdelmalek Atia *,Kamal Mohammedi

1 Univ.Boumerdes,LEMI Lab,Algeria

2 Univ.El-Oued,LEVRES Lab,Algeria

Keywords:Lattice Boltzmann method Density driven Pore-scale CO2 Mass transfer

ABSTRACT Saline aquifers are chosen for geological storage of greenhouse gas CO2 because of their storage potential.In almostallcases of practicalinterest,CO2 is present on top ofthe liquid and CO2 dissolution leads to a smallincrease in the density of the aqueous phase.This situation results in the creation of negative buoyancy force for downward density-driven natural convection and consequently enhances CO2 sequestration.In order to study CO2 injection atpore-level,an isothermal Lattice Boltzmann Model(LBM)with two distribution functions is adopted to simulate density-driven natural convection in porous media with irregular geometry obtained by image treatment.The present analysis showed that after the onset of natural convection instability,the brine with a high CO2 concentration infringed into the underlying unaffected brine,in favor of the migration of CO2 into the pore structure.With low Rayleigh numbers,the instantaneous mass flux and total dissolved CO2 mass are very close to that derived from penetration theory(diffusion only),but the fluxes are signi ficantly enhanced with high Ra number.The simulated results show that as the time increases,some chaotic and recirculation zones in the flow appear obviously,which promotes the renewal of interfacial liquid,and hence enhances dissolution of CO2 into brine.This study is focused on the scale of a few pores,but shows implications in enhanced oil/gas recovery with CO2 sequestration in aquifers.

1.Introduction

The growing concerns over the global warming due to the increase in the global concentration of greenhouse gases in the atmosphere have increased the interest in examining various techniques to reduce the emission of these gases.A main component of greenhouse gases is carbon dioxide(CO2).A promising long term solution for mitigating global heating is to inject CO2into geological formations;either for CO2sequestration or enhanced oil/gas recovery.A suitable choice of geologicalformations for CO2sequestration includes coalbeds,depleted petroleum and gas reservoirs,deep-sea sediments,and deep saline aquifers.In many practical cases,geological storage of CO2is accomplished by injecting it in supercritical phase into a porous rock formation below the earth surface which is already saturated with a liquid(water or oil)[1,2].The great concern on geological CO2trapping is the long-term fate and security of the stored CO2.More speci fically,researchers are concerned with preventing potential leakage of CO2from underground reservoirs,which mightbe caused by fractures in cap rock orabandoned oil/gaswells[3].Therefore,the analysis ofCO2dissolution phenomena in underground geological formations is of great importance in CO2injection projects.

In most CO2injection studies over the world,the carbon dioxide is injected in the supercritical phase.The supercritical density of CO2is about 70%of water density,and the kinematic viscosity is 0.1-0.25 of water[4].Hence,itis expected thatwhen CO2is injected into the porous formation,initially it accumulates under low-permeability cap rock and subsequently dissolves into the formation liquid by molecular diffusion[5].As a result,the saline density aquifer increases and eventually the CO2-brine interface becomes unstable.For favorable conditions,density driven natural downward convection occurs and CO2-saturated brine moves downward and is replaced by underlying unaffected brine,which enhances the mass transfer of CO2.Numerous experimental and theoretical studies on the density-driven natural convection process in porous media[2,6-8]are available.Rapaka et al.studied the effect of anisotropy in the permeability fields on density-driven natural convection[9].Li and Dong reported experimental studies of carbon dioxide diffusion in porous media under reservoir conditions[10].They presented a new method for measuring the effective CO2diffusion coefficientin oil-saturated porous media under reservoir conditions.The investigation of the effect of heterogeneity on density-driven natural convection of CO2overlying a brine layer was studied by Farajzadeh et al.[11].The spectral method has been used to generate permeability heterogeneity fields.Chen etal.investigated experimentally the velocity distribution in Rayleigh convection in the gas-liquid masstransferofacetone volatilization in acetone-ethyl acetate binary system via particle image velocimetry(PIV)[12].

During the past two decades,Lattice Boltzmann Method(LBM)has been demonstrated to have the advantages in clear mathematical and physical meanings,easy implementation of the boundary conditions in complex geometry like real porous structure,high computational efif ciency,numericalaccuracy and easy parallelization.LBMhas been successfully applied to the simulation of fluid flow including single-and multiphase flows[13,14],natural convection caused by temperature gradient[15,16],double diffusive natural convection in a cavity with a hot square obstacle inside[17],Rayleigh convection generated by a local high concentration gradient in the gas-liquid mass transfer process of CO2absorption into liquid ethanol[18,19],and density driven natural convection in anisotropic and heterogeneous saline aquifers[4].Cheng etal.[3]have simulated the density driven naturalconvection in simple fractured porous media at pore scale by using LBM.However,in real porous media the geometry is rather irregular and complex.

In this study,a two-dimensional LBM with a double distribution model is implemented for simulating density driven natural convection of CO2in porous medium of irregular geometry.The density driven natural convection was manipulated by a high concentration gradient area at the top boundary.The complicated coupling between pore geometry,dissolution of CO2,and density-driven natural convection was investigated.The simulated results are used to investigate the velocity,concentration fields,the instantaneous mass flux and total dissolved CO2mass.

2.Problem De finition and Mathematical Model

2.1.Pore-level view of CO2 injection project

Fig.1 shows a schematic of pore-level image during CO2injection in saline aquifers.Porous rocks in the geological formation have a wide range of pore scales.Injected CO2accumulates under the low permeability cap rock and slowly diffuses into brine.The density of CO2-saturated brine is slightly increased,as a result the density driven natural convection occurs and decreases.However,there is little pore-scale investigation of the process coupling density driven natural convection and CO2dissolution in real porous structure.In the following sections,the LB method will be described for the porelevel study.

Fig.1.Schematic illustration of pore-level CO2 injection,dissolution and natural convection.

2.2.Problem statement

The two-dimensional computational domain of the problem is shown in Fig.2.The geometry is considered as a porous structure packed with irregular grains retrieved by image processing and filled with quiescent brine.The domain has a total dimension of 200×200 nondimensional lattice unit with a resolution of 0.2 mm per unit.Thus,the domain dimension(h×h)in the real system is 0.04×0.04 m2.The thermal effect of absorption of CO2is assumed to be isothermal.Also,the chemical interactions and diffusion between rocks/ fluids are neglected,so the dissolved CO2acts as a conservative solute tracer.The top boundary of this domain is CO2-brine free interface.The bottom boundary is considered as solid wall.Because the natures of velocity and concentration distributions are complex and random in the density driven natural convection process,we simulate a simple case where the top boundary is a free surface with a constantconcentration corresponding to the solubility of carbon dioxide.

Fig.2.Geometry of the computational domain of porous structure.

The physical properties of the brine and the saturated brine-CO2are given in Table 1.In the simulation,an interfacial concentration value of 0.4 kg·m-3at the CO2-brine interface was taken.

Table 1 Physical properties of pure brine and saturated brine-CO2 at T=296.15 K,P=101.0 kPa

2.3.LBM model with double distribution functions

Unlike macroscale differential approaches based on traditional Navier-Stokes equations,the LBM discretizes the fluid into particles which are individually sited on every node of a lattice system.At the end of each time step,the particles move and collide with others,and the distribution probability of particle density at a node after the collision obeys the Boltzmann equation.The underlying concept of the LBM is to incorporate the essential physics of the problem into the simpli fied kinetic equation such that the correct macroscopic behavior of the fluid is recovered[18].A D2Q9(9 velocity vectors in 2-D space)lattice model was adopted in this study(Fig.3).The LB equation for fluid flow is written as[21]

Fig.3.Two-dimensional nine-velocity(D2Q9)model.

where fi(x,t)is the volume-averaged distribution function for the particle at lattice location x and time t traveling along the i-th direction;and eiis the lattice velocity vector corresponding to direction i,which is de fined as

In the above equations,c= δx/δt,where δx and δt are the lattice spacing and the time step,respectively,is both equal to one in LB method.τfis the dimensionless relaxation time related to the kinematic viscosity ν as ν =(τf-0.5)δx2/3δt[22].

The pressure can be calculated by p=ρ.In order to simulate the transport of dissolved CO2,another distribution function gi(x,t)is used:

whereτsis the dimensionless relaxation time related to the solute diffusivity by D=(τs-0.5)δx2/3δt[22].(C,u)is the equilibrium distribution function atlocation x and time t along the i-th direction,which is chosen to recover the macroscopic advection-diffusion equation:

where C is the macroscopic concentration ofdissolved CO2calculated by

The main governing dimensionless parameters controlling density driven natural convection are the Rayleigh(Ra)and Schmidt(Sc)numbers,de fined as

One of the important characterizing parameters for density driven natural convection is the critical time,de fined as the minimum time for the onset of instability at the interface and marked as tcr.We can define the dimensionless critical time aswhere tcis the characteristic time and calculated by

One of the key components for the application of the lattice Boltzmann method to physical problems is the correct conversion from physical world units to lattice units to and vice versa.Because it is not possible to directly input the variables therefore,scaling has to be applied.The fluid viscosity has to be scaled first,we can choose any value less than 0.1[23].Having in hand the real values of fluid physical proprieties and the model parameters,we can de fine the domain high h[m],the fluid kinematic viscosity ν[m2·s-1],the strength of external force(gravity)g[m·s-2],the fluid density ρ[kg·m-3],and the size of an LBM lattice Δx[m](resolution dependent).These parameters are used to set a reasonable time step size Δt[s]and limit the mass difference in the simulation Δm[kg].The values of the LBM simulation are dimensionless and they should be converted to the dimensionless values.The following steps show how the physical units(with a subscript'p')are transformed to lattice units(with a subscript'l'):lattice time:tl=tp/Δt,lattice gravity:gl=gp(Δt2/Δx),lattice density difference:Δρl= Δρp(Δx3/Δm),lattice diffusion number:Dl=Dp(Δt/Δx2).

2.4.Implementation of external forces

Initially,dissolved CO2migrates downward mainly by pure diffusion,which implies that the dissolved CO2slightly increases the brine density.The well-known assumption of Boussinesq approximation is adopted to assume thatmaterialproperties are independentofCO2concentration except the density.Thus,only the body force,buoyancy or gravity,caused by density difference is considered as the external force in the system.The density difference,which was assumed to vary linearly with the concentration difference caused by the mass transfer,is the only source[18].In this problem,the exerted body force can be written as

whereρis the saline aquifer density,ρ0is the density ofpure water,and β is the expansion coef ficient,which is de fined as β =(Δρ/ΔC)/ρ0.It is noticed that only the second term on the right-hand side of Eq.(12)contributes to density variation.The first term on the right-hand side of Eq.(12)is the same for all nodes,thus we can put it into the pressure term as- ρ0g ey.After inserting β into Eq.(12),the body force term can be reduced to

where Δρ is the incrementof saline aquifer density when the difference between the concentrations of pure brine and the concentration of saturated brine-CO2reaches the maximum value ΔC.The dimensionless form of the body force is de fined as follows

This implies that the body force is linearly proportional to the local dimensionless concentration,θ=C/ΔC,with 0≤θ≤1.Incorporating the body force into the model,the velocity in the equilibrium distribution function u∗and flow velocity ũ is modi fied as follows[24]:

2.5.Boundary conditions

In the LB method,distribution functions out of the domain are known from the streaming process.The unknown distribution functions are those towards the domain[25].Fig.4 shows the unknown distribution functions,which are needed to be determined,as dotted lines.

The bounce-back boundary is employed on the south boundary and brine-rock interface for the flow field and the concentration field,which means that incoming boundary populations are equal

to out-going populations after the collision.For instance for south boundary,the following conditions are imposed:

Aperiodic boundary condition was used atthe westand eastboundariesforthe flow field and the concentration field,which meansthatentering distribution functions at west boundary,are the same as the distribution functions leaving from east boundary and vice versa.For instance for west boundary,the following conditions are imposed:

We use the same boundary condition for concentration field.We use the mirror symmetric boundary[26]to treat the CO2-brine free interface.This type of boundary condition can ensure that the gradient of horizontal velocity along the vertical direction is equal to zero,i.e.

where m is the number of nodes in the y direction.

The constant concentration boundary[27]is used on the north boundary to handle the dissolution of CO2into brine,which is described below:

Fig.4.Boundary conditions with domain classi fication:(1) fluid lattice(light gray nodes);(2)inner solid lattice(gray nodes)and(3)boundary solid- fluid interface(black nodes).

3.Model Validation

The LB model described above was validated with experimental study performed by Fu et al.[18]by simulating under the same operating conditions the 2D Rayleigh convection generated by a local high concentration gradient in the gas-liquid mass transfer process during CO2absorption into liquid ethanol.The height and width of the computational domain are equal to 100(in the lattice unit);with a resolution of 0.4 mm.Thus,the domain dimension in the real system is equal to 40×40 mm2.For boundary conditions,the bounce-back boundary is employed on the south boundary,and periodic boundary on the east and west boundaries for the flow field and the concentration field.On the north boundary,we used the mirror symmetric boundary for fluid flow and the constant concentration boundary for solute transport.The predicted distribution of vertical velocity along vertical line during the absorption process at x=20 mm for different time instants is presented in Fig.5.It can be seen from this figure that simulated results coincide quite well with the experimentalresults.The minor discrepancies between the present simulation and experimental measurements can be attributed to the measurement uncertainties and simulation accuracy.Fig.6 shows the simulated and experimental transient streamline results after 120 s,140 s and 160 s,respectively.These figures show the agreement between simulation and experimental results of circulation flows obtained by particle image velocimetry technique.

4.Simulation Results and Discussion

As described in Section 2.4,the brine density was assumed to be linearly proportional to the concentration of CO2caused by mass transfer process at the interface.Hence,saline aquifer density near the top was higher than that near the bottom.Thus,the density driven natural convection will be obvious and develop if the Rayleigh number is higher than the critical number.Under the real physical proprieties listed in Table 1,the Schmidt number will be around 500,however,this value presents some numerical instability.Hence we set Sc=100 to avoid this problem.In addition,the computational domain has a total dimension of 200 × 200 and we set ν =0.1 for LB simulation.The simulation results of the temporal-spatial of solute concentration contours at different times are illustrated in Fig.7,representing the density driven natural convection of CO2in porous structure.It can be seen that in the early stage,the dissolved CO2migrated into the porous medium mainly by diffusion and the liquid density near the interface is increased.In this scenario,the density of the brine at the interface increases and becomes unstable and tended to move downward.

The density driven naturalconvection could be observed at t=50 s.After that,CO2continues to diffuse and intrude into the underlying unaffected brine.This process increased the interfacial area between dissolved CO2plume and unaffected brine,which is favorable to enhance the migration of dissolved CO2into the porous structure.The density driven natural convection reached the bottom of the domain at about t=200 s,after that it began to spread into the porous structure,as shown in Fig.7g,h and i.

The transient velocity vectors in the porous structure are shown in Fig.8.Itcan be seen in the early stage,thatthere are severalpairs ofvortex cellswith diverse rotationaldirectionsand differentsizesare formed in the brine.Also,with time going on,the flow fields became more disordered.As a result,the vortex centers continue moving downwards and impel the continuum mixing and exchange between the saturated brine and unaffected liquid by interfacial vicinity and promote the mass transfer process.Hence,we can conclude that in reality,the density driven natural convection during a CO2injection project is able to produce obvious irregular motion of this geometry type.

Fig.5.Comparison between simulated(continuous line)and experimental(black points)results for distribution of vertical velocity along vertical line at x=20mm for different time instants.

The criticaltime for the occurrence ofdensity driven naturalconvection as function of Rayleigh number is presented in Fig.9.tcrdecreases when Ra is increased.Speci fically,for the case of Ra=5.7×1012,tcr≈30 s.Several other Ra numbers were used in the simulations.The shape of the tcr-Ra curve was similar to the result obtained by Chen et al.[3],who investigated the relationship between the critical time and Ra number in fractured porous media.

Fig.6.Transient simulated streamlines of present study(left)and transient experimental streamlines(right)for different time instants.

The simulated dimensionless critical time is presented as function of block size,as shown in Fig.10.It can be observed that there is a nearly linear relationship between-1/2 and block size for our problem.

Fig.11 shows the transient behavior of velocity vectors for different values of block size at t=200 s on a grid of 200×200.As can be seen,a set of circulation flow is formed in the liquid bulk.Also,it can be noted that the geometrical morphology of porous media is invariant with variation of block size and the flow patterns are geometrically similar to Fig.8.

The effect of Rayleigh number on interfacial mass transfer can be presented quantitatively by the instantaneous mass flux and dissolved CO2mass into the porous structure.From the simulated concentration field,the instantaneous mass flux Nins,tacross the interface at time t under the condition of CI=0.4 kg·m-3can be estimated by

The instantaneous mass flux through the top boundary per unit cross-sectional area along the vertical direction can be obtained by penetration theory,and can be written as[28]

The totaldissolved CO2mass accumulated after time t per unitcrosssectional area can be obtained by integrating Eq.(22)to yield

Fig.7.Simulated of temporal-spatial cloud map of solute concentration at different times.

The variation of instantaneous mass flux across the interface with time for different Rayleigh number by penetration theory and simulated concentration is shown in Fig.12.In all cases the instantaneous mass flux along the vertical direction across the interface decreased clearly before 40 s where the interfacial mass transfer only proceeds by molecular diffusion.It can be found that the simulation results agree very well with the values predicted by the penetration theory within that time period,and then increased slowly with time atthe onsetofdensity driven natural convection.Following that,the instantaneous mass flux across the interface after 40 s is very close to the one derived from penetration theory with decreasing Rayleigh number,but it is above that derived from penetration theory with increasing Rayleigh number.As a result,the interfacial mass transfer proceeded in the combined manners of both density driven natural convection and molecular diffusion at high Rayleigh number.Therefore,the interfacial mass transfer flux is enhanced with increasing of Rayleigh number.

Fig.13 shows the effect of Rayleigh number on dissolved CO2mass per unit interface area as a function of time.In the case of Ra=4.56×109,the dissolved CO2mass is close to thatobtained by analytical solution,or the diffusion dominated transport of dissolved CO2.When Ra was increased to 7.12×1010,the dissolved CO2mass in the system was enhanced,because of the intensive downward convection due to the density difference.However,in this case,the dissolved CO2mass in the porous structure does not exceed 0.001 kg·m-2over that by combined diffusion-convection and diffusion-only.As Ra was increased to 4.15×1011,the difference in the total dissolved CO2mass became signi ficant.This implies that the highly disordered fluid motion greatly enhanced the dissolution of CO2at the different regions of unaffected brine.

5.Conclusions

In this study,the LB model with a double distribution functions has been used to simulate density driven naturalconvection in a porous medium with a complex geometry.The LB model was veri fied by experimental study performed by Fu et al.[18],by simulating the Rayleigh convection generated by a local high concentration gradient in the gas-liquid mass transfer process during CO2absorption into quiescent liquid ethanol.The key point for this study lies in the integration of irregular geometry obtained by image treatment in order to investigate real systems.The simulated results with the adopted LBM model show that the density driven natural convection will not occur until a critical time,and this time decreases with increasing Rayleigh number.When critical time is reached,the onset of convective instability will make the brine with a high CO2concentration intrude into the underlying brine by circulation flows with high velocity and vortexes,which increases the interfacial mass transfer and mixing between rich brine by CO2and poor brine,consequently in favor of the migration of CO2into the fracture of porous structure.It was also demonstrated that our adopted LBM simulation approach can effectively capture the instantaneous mass flux across the interface and total dissolved CO2mass,and the mass transfer rate and dissolved mass are effectively enhanced by increasing Rayleigh number.Pore-scale simulation by lattice Boltzmann method is considered as a promising tool for studying the coupling between CO2dissolution,mass transport,and heterogeneity ofpore geometry.Speci fically,the effect of density driven natural convection during geological storage of CO2on underground petro-physical and chemical properties should be investigated.

Fig.8.Simulated transient velocity vectors at different times with enlarged map of rectangular zone.

Fig.9.Critical time for the onset of instability as a function of Rayleigh number for μ=0.93 × 10-3Pa · s,D=1.85 × 10-9m2·s-1 and Δρ =10 kg·m-3.

Fig.10.Dimensional critical time as function of block size(grid 200 × 200)for μ =0.93 × 10-3Pa · s,D=1.85 × 10-9m2·s-1 and Δρ =10 kg·m-3.

Nomenclature

Fig.11.Transient velocity vectors for different values of block size at 200 s(grid 200×200).

Fig.12.Variation of instantaneous mass flux N ins,t with time for LBM simulation and penetration theory at C I=0.4 kg·m3 for different Rayleigh numbers.

Fig.13.Variation of dissolved CO2 mass with time at C I=0.4 kg·m3 for different Rayleigh numbers.

Subscripts