Water Flooding and Viscous Fingering in Fracture and Porous Media by Lattice Boltzmann Method

Y. Shiri,a H. Hassani,b,* M. Nazari,c and M. Sharifia aDepartment of Petroleum Engineering, Amirkabir University of Technology (Polytechnic of Tehran), Tehran, P.O. Box: 15875-4413, Iran bDepartment of Mining and Metallurgy Engineering, Amirkabir University of Technology (Polytechnic of Tehran), Tehran, P.O. Box: 15916-34311, Iran cFaculty of Mechanical Engineering, Shahrood University of Technology, Shahrood, P.O. Box: 36199-95161, Iran


Introduction
Fluid flow in porous media has a very important role in various industries, such as resin transfer modeling, filter analysis, groundwater transfer, biomechanics, and oil recovery.Replacing a fluid with another fluid for the preferential flow path of injecting fluid, which creates its own unique pattern of fluid flow according to the environmental and fluid characteristics, is one of the challenging issues in industries, especially in oil resource recovery.For instance, tongues of water are of limiting parameters in EOR.Fluid front and viscous fingering in immiscible oil displacement are very important.When a less viscous fluid pushes a fluid with a higher viscosity, the interface between the two fluids becomes unstable and the thickness of the fluid front oscillates.As a results, it is not always consistent with its estimated value [1][2][3] .Disturbances occur at the interface of the two fluids in both the laboratory and the numerical study [2][3][4] .In addition, the stability of the fluid front, because of the continuous development of the flow across the unstable state that causes droplets, is very important.The viscosity of two fluids, surface tension, wettability, the rel-ative velocity of fluids, porous media structure and geometry of grain are the determining factors in viscous fingering and fluid front shape 5,6 .
Cinar et al. have investigated the two-phase flow regime and the effect of gravity, viscosity, capillary force and surface tension of two fluids in two parallel porous layers by finite difference and streamline methods, and they also experimented it on glass beads 7 .Maaref et al. have investigated the effect of wettability, heterogeneity, and viscosity differences in the process of water and oil displacement in glass micro-model and its numerical solution with finite element method.The results showed a strong correlation between the numerical method and laboratory observation, and the increase in EOR caused by higher viscosity and wettability of the injecting fluid as well as the decrease in heterogeneity in porous media 8 .
In the absence of any analytical solution for the fluid flow, the numerical simulation will work.The most noticeable method is discretizing of equation by using finite difference and finite element methods on the macroscopic scale as well as the useful Lattice Boltzmann method (LBM) in the mesoscopic scale 9 .Different Computational Fluid Dynamic (CFD) techniques including Volume of the Void (VOF) and the Level-set method are used to simu-late two-component or two-phase fluids, but their interface can only be assessed in large-scale, and the information is usually ignored in small scales.LBM has recently become one of the most useful methods for calculating fluid dynamics.The convenience of using a variety of boundary conditions along with intrinsic network has valued LBM for simulating flow and heat transfer in complex systems with special geometrical dimensions [10][11][12][13] .Some advantages of LBM are its ability to evaluate the factors involved in viscous fingering in the mesoscopic scale, reduce computational time, perform parallel computations, be implemented easily, and continuously track the fluid front within the time rather than other macroscopic and microscopic methods 14 .The mesoscopic nature of LBM coupled with the easy use of the boundary conditions on the solid wall has made it a powerful tool for simulating fluid flow in porous media in small scales.Also, this method has the ability to simulate fluids both in microscopic and macroscopic scales 15 .Different models of LBM, such as Color Gradient (Rothman-Keller) Method, Pseudo-Potential (Shan-Chen) Model, Free Energy Model, and Mean-Field Model have been developed for two-phase fluid simulation [16][17][18][19] .In the Color Gradient Model, a local gradient in the phase field is used to separate the phases, and the advantages of it are the independent use of viscosity ratio, contact angle, surface tension, and density ratio 14,20 .Pseudo-Potential (Shan-Chen) Model uses pseudopotential coefficient for particle interaction, and its advantages are simplicity and computational efficiency.In this method, surface tension, density, and viscosity ratios are not independent 21 .In the Free Energy Model, the phases' interaction directly defined by the general equilibrium distribution function, which includes the pressure tensor in the collision process.In the Mean--Field Model, the appropriate force is extracted from fluid phase behavior.LBM is one of the hottest topics in recent years and Shan-Chen Model is the main body of simulation in these studies 22 .
At the time of injecting a low viscous fluid into a higher viscous fluid, due to the applied force with the instability in the interface of two fluids, finger-like patterns form, and its shape depends on several factors.Tracking and investigating effective factors of fluid front in viscous fingering in a microfracture is the most common model due to the existence of their analytical solutions.Various studies in the Hele-Shaw microchannel, which represent the smallest component in complex geometries, have been performed experimentally 6,[23][24][25] and numerically using LBM [26][27][28][29][30] .Studies have shown that the fingertip is narrowed by increasing the capillary number and the viscosity ratio, as wells by decreasing the wettability of the injected fluid.It is promis-ing that the Shan-Chen Model of LBM has the ability to study the fluid front 31 .In addition, the study of the fluid front in porous media is also important.In recent years, remarkable studies have been carried out in this regard and some parameters including the viscosity ratio, velocity, surface tension, and adhesive force have been investigated [32][33][34][35][36][37][38] .The purpose of the present study is to use Shan-Chen Model of LBM and investigate its ability to accurately examine the main factors affecting the fluid front in a microfracture and a simple porous media.

Lattice Boltzmann method
Lattice Boltzmann equation LBM can be simplified by the Lattice Gas Automata (LGA) or by the finite difference solution of the Boltzmann equation.Both LBM and LGA can be obtained from the Navier-Stokes Equations on the macroscopic scale.For the first time in 1964, LGA was used to simulate fluid flow 39 .However, this method did not yield good results due to its statistical noise and its dependence on velocity.Later in 1988, LBM was introduced for simulating the fluid flow, which apart from the advantages of LGA, was numerically more efficient and more appropriate to simulate flow with high Reynolds number and more applicable in glacier flow 40,41 .Then, LBM was obtained by discretizing and solving the Boltzmann equation 42 .

Shan-Chen model of LBM
Shan-Chen model of LBM is used to simulate fluid flow with several phases and components, and is known as the Pseudo-Potential model that is based on the concept of the intermolecular interaction of fluid particles.Shan-Chen model of LBM for each phase is in equation (1), and is updated after initialization in a cycle consisting of collision, streaming, applying boundary conditions, and determining fluid variable until the stop criterion is achieved. (1) The density, velocity, and distribution function of each phase can be determined by: In this method, it is assumed that the total force to one of the phases which consists of the effective nearest neighbor, gravity, and adhesion forces ( ) are as follows 43 . (4) In Equation ( 4 (7) For densities much smaller than 1, this function is equal to the density ( 1( ) r ψ r r << ⇒ = ).The total force can be used as follows: In this equation, u′ is the velocity of all phases and is defined as follows: (9) By considering intermolecular force, fluid behavior will be non-ideal, and the equation of state is as follows: In this equation, the right term is comprised of the ideal and non-ideal fluid behavior, respectively.

Laplace test and surface tension estimation
Among the fluid parameters, the surface tension is directly related to G and the density profile at the interface of the two fluids, and they are expressed by the equation of state 44 , but it cannot be determined directly and requires simulations.Taking into account simulation with circle droplets in different radii, the surface tension is equal to the slope of the inside and outside pressure differences of the droplets and the inverse of the radii 45 .The periodic boundary conditions are on the four sides of the model, and a simulation is started by considering a square droplet of the first fluid inside the second fluid.After reaching the equilibrium condition, the radius of the circle and inside and outside pressure are calculated.By repeating this simulation and calculating the pressure difference in different radii, the surface tension is obtained from the slope of the pressure difference versus inversion of the radius.To calculate the pressure, Equation ( 10) is used.To measure the radius graphically, the density graph of the two fluids is drawn from the center of the network toward the outside direction by radius of and the intersection point is equal to the radius of the droplet.The following formulae is also used to calculate the radius of the droplet in which r is the average densi- ty in the entire network for the first fluid, or the fluid that forms the droplet 46 .
The variations in the surface tension are dependent on the G and the density of the fluids.Smaller G causes and bigger G causes instability.Therefore, the size of G must be large enough not only to avoid the thickening of the interface, but also to stabilize it enough.The results of this simulation for two fluids with a density and relaxation time of 1 and G = 3 are shown in Fig. 1, and the estimated surface tension is 0.066.

Laplace test, Contact angle determination
By substitution of surface tension and molecular adhesion of two fluids in the Yang equation (Equation ( 12)) with the parameters of LBM, the Equation ( 13) is obtained.In equation ( 13), r 1 and r 2 are the densities of the first fluid and the second dissolved fluid, respectively.For the values of G ads1 and G ads2 , different intervals exist, but the least error in using this equation and LBM simulation occurs when the adhesion parameters of the two fluids are equal with opposite signs.The positive and negative signs indicate repulsion and attraction, respectively.The results of the contact angle determination with various G are shown in Fig. 3.In this simulation, a rectangular droplet with the size of 4 % of cross-sectional area is placed in 100 × 200 lattice network.The density of each fluid is 2 and the density of 0.06 is considered for dissolved fluid.The viscosity of both fluids is 0.1667.The results of this simulation are in complete correlation with the results of Huang et al. 47

Grid independence test
The effect of grid size on the simulation has been considered with four sizes of porous media.The results of this simulation with Reynolds num-ber of 168 for fluids with surface tension of 0.666, which resulted from the G of 3, relaxation time of 1, density of 1 for both fluids (they are equivalent to viscosity ratio of 1), and adhesion coefficient of 0 or 90 degrees contact angle (they are equivalent to natural wettability) are shown in Fig. 4. As it can be seen, in four grid sizes of 255×255, 561×561, 937×937, and 1403×1403, the fluid front only had a slightly different shape in grid size of 255×255.This deviation can also be seen in Fig. 4(e), in which the water saturation is plotted versus dimensionless time (t D = t/t breakthrough ), and there are no significant changes in other sizes.Therefore, these figures suggest that the grid sizes above 561×561 have negligible impact on water saturation and fluid front shape.

Flow regimes
The fluid front in porous media depends on a variety of factors.In this study, the viscous fingering is considered according to three parameters of the displacing fluid, including the discharge or capillary number ( / Ca V rυ s = ), dynamic viscosity ratio (µ rυ = ), and surface tension (s).Fig. 5 shows the viscous fingering in a microfracture, in which L and H are the fracture's length and width, respectively.W in is the internal width, W out is the outer width, T is the length, and S is the length of the contact line of the viscous fingering.
Various parameters affect the viscous fingering.Fig. 6 shows the distribution of the fluid front shape with respect to the changes in the viscosity ratio and the capillary number that is divided into three different regions 48 .These results were also obtained by Color-Gradient model and Free-Energy model of LBM simulation and experimentally later 34,49,50 .
In the stable displacement region, the main force is derived from the viscosity of the displacing fluid, and the capillary effect and the pressure drop  48 .
In the viscous fingering region, the main force is applied by the displaced fluid, and the capillary effect and the pressure drop are negligible.The viscous fingering is similar to tree branches without loops and extends to the entire porous media toward the outlet 48 .
In the capillary fingering region with low capillary numbers, viscosity is not noticeable for both fluids, and the main force is due to the capillary effect with a finger-like pattern, porous media extends in all directions, even to the opposite side and to- wards the entrance, which form loops of various sizes in which the displaced fluid is being trapped.The three boundaries among these three regions are dependent on pore size and porous media, but the overall shape does not change 48 .
Several studies have been conducted in this regard.According to the results, for a viscosity ratio more than one in low capillary numbers, the fluid is displaced by capillary fingering, in which the fluid tends to move towards the wider pores in different directions, and the fluid may move perpendicularly and even backward of the fluid flow direction.By the increase in capillary number, the displacement becomes stable.For a low viscosity ratio less than one in low capillary numbers, the capillary fingering is regenerated, while viscous fingering will prevail in high capillary numbers due to the injection of very low viscous fluid into the highly viscous fluid.In this case, there are several connected and disconnected flow paths toward the outlet in which the capillary force and the pressure drop in the fluid movement are negligible.The main force of the movement is the viscous force.In moderate capillary numbers, both the phenomena are observed and the region is known as Crossover Zone 34,51,52 .The simulation of current study takes place in the viscous fingering and in the transient zone of Fig. 6.The results of microfracture were consistent with Hele-Shaw microchannel results [26][27][28][29][30] and the same conditions were followed in the simulation of porous media.

Viscous fingering and flow rate
In immiscible displacement, the capillary force is dominant in the low discharge, leading to the formation of blobs trapping.By increasing capillary number and flow rate, the viscous effects are dominant and cause viscous fingering 53,54 and it is effective in the displacement efficiency at the breakthrough time.According to the study of Liu et al., the increase in capillary number due to instability in the fingers of the fluid front, the injected fluid flow is not in a quasi-steady state and the slow decrease in displacement efficiency is observed, then, with the increase of the capillary number with quasi-steady state of flow, the displacement efficiency increases linearly 32,49 .
In this simulation, the Fracture width is 40 lattices, the Fracture length is 200 lattices, and the porosity of porous media is assigned to 90 %.A constant velocity boundary scheme proposed by Zou and He was used at both left inlet and right outlet of the model 55 and the upper and lower and grains' solid wall was half-way bounce back boundary condition.The results of this simulation with velocities of 0.001, 0.005, 0.01, 0.05, 0.1 for fluids with surface tension of 0.666, which resulted from the G of 3, relaxation time of 1, density of 1 for both fluids, which are equivalent to viscosity ratio of 1, and adhesion coefficient of 0 or 90 degrees contact angle, which are equivalent to natural wettability, are shown in Fig. 7.By applying very low velocities equivalent to low capillary numbers, the fluid front will have a slight curvature, and by increasing the flow rate, this curvature increases and breakthrough time decreases.As it can be seen, for better visualization of fluid front of the channel in higher flow rate, the contour lines of fluid front are plotted in 100 time steps instead of 500 time steps used in lower flow rate.Therefore, the increases in the velocity will persuade the fluid front from the piston form to the viscous fingering.However, as it can be seen, the amount of trapped droplets of the saturated fluid depends on the velocity, and most of them are observed at the velocity between the piston form front and the viscous fingering.

Viscous fingering and viscosity ratio
In this simulation, the Fracture width was 40 lattices, the Fracture length was 200 lattices, and the porosity of porous media was assigned to 90 %.A constant velocity boundary scheme proposed by Zou and He was used at both left inlet and right outlet of the model 55 and the upper and lower and grains' solid wall was half-way bounce back boundary condition.The density of both fluids was 1, and the relaxation time of the displaced and displacing fluids were 1 and 0.666, respectively.The viscosity ratio changed between 1/5 and 1, the capillary number was fixed at 0.063, the wettability of the fluids was neutral (adhesion coefficient of 0 or 90 degrees contact angle), and G was determined by the stability of the simulation.Due to the changes in G and in the viscosity of the displacing fluid, the suitable velocity by fixing the capillary number was calculated and applied (Table 1).If the dynamic viscosity ratio, which is the ratio of the displacing fluid to the displaced fluid, is greater than 1, the displacement  is called "favorable", and if it is less than 1, the displacement is called "unfavorable" 48 .The results of this simulation by various viscosity ratios in the fracture and porous media are shown in Fig. 8. Viscous fingering length increases in both fracture and porous media.The breakthrough time decreases with the increase in the dynamic viscosity differences.It is important to note that Shan-Chen model of LBM is sensitive to the dynamic viscosity ratio, and it is unstable in viscosity ratios above 5.Fluid-fluid interaction coefficient (G) lower than its critical value will yield a stable solution and higher values rather than this value rising until numerical instability occurs 43 .By higher density ratios, G must be reduced to avoid instability and by Shan-Chen model of LBM (equation ( 7)), this ratio should be less than 5.This limitation has been partially removed using different equations of state 56,57 and using different models of LBM, such as phase-field model 58 .
According to the dynamic viscosity ratio, three zones are defined.In zone I, in which the dynamic viscosity ratio is very low, the viscous force of the displacing fluid in comparison to the viscous force of the displaced fluid is negligible.In zone II or the transitional zone in which the pressure drops, the viscosity of both fluids plays an important role.In zone III, dynamic viscosity ratio is high, the pressure drop due to the viscosity of displaced fluid is negligible and viscous force of the displacing fluid is dominant 48 .In a constant capillary number with low viscosity ratio and without gravity force, displacing a wetting fluid by a non-wetting fluid causes a thin viscous fingering pattern and occupies only small pore spaces.In some pore spaces, the displacing fluid is trapped, and is known as "daughter bubbles".By increasing the viscosity of displacing fluid or viscosity ratio, the thickness of the fingers and the daughter bubbles' size increase.With continuous increase, all pores will fill and the displacing fluid will not be trapped in the pores, and therefore, the capillary fingering or stable conditions will form 49 .According to Fig. 6, the simulation of this study takes place in the viscous fingering and in the transient zone.As it can be seen, with the decrease in the dynamic viscosity ratio, the viscous fingering becomes very thin and takes up only small pore spaces, while the trapped bubbles of displaced fluid decrease.

Viscous fingering and wettability
Considering the adhesion force between the fluid and the solid surface, its effects on the fluid front are important.Fluid injection, by assuming that the fracture and porous media are saturated with the second fluid, with different wettability equivalent to contact angles of 30, 60, 90, 120, 150 degrees, is shown in Fig. 9.In this simulation, with the dynamic viscosity ratio of (1/3), the density of both fluids was 1, and the relaxation time of the displaced and displacing fluids were 1 and 0.666, respectively.The fluid-fluid interaction coefficient (G) was 2, the width of the fracture was 40 lattices, the length of the fracture was 200 lattices, the porosity of porous media was 90 %, and the injection velocity was 0.021.A constant velocity boundary scheme proposed by Zou and He was used at both left inlet and right outlet of the model 55 and the upper and lower and grains' solid wall was half-way bounce back boundary condition.As it can be seen, contact angle or wettability was another factor affecting the shape of the fluid front due to its effects on the capillary force and the non-tendency of nonwet fluids to enter the pore spaces.By injecting fluid with neutral wettability, smaller bubbles of displaced fluid were trapped in the displacing fluid, and therefore, the dendritic structures had decreased.By increasing the wettability and decreasing the breakthrough time, these bubbles were absent and the flow regime became stable in piston form.The results are consistent with the results of Cottin and Liu et al. 49,51 Liu et al. showed that in gas injection by the viscous fingering mechanism at the breakthrough time with different wettability, the sweep efficiency of the gas was a function of the wettability.By increasing the wettability from non-wet to neutral, the sweep efficiency decreased, and from a neutral fluid to wet, the sweep efficiency increased.In addition, the lowest sweep efficiency occurred with neutral wettability.The sweep efficiency obtained from this simulation (Fig. 10) confirmed this, and was in accordance with the results of Dong et al. 32,49 Water flooding in a fracture and porous media saturated with Arvandan crude oil In this simulation, the width of fracture was 40 lattices, the length of the fracture was 200 lattices, porosity of porous media was 90 %, and the distilled water as the displacing fluid was 1, and the density and the dynamic viscosity of Arvandan crude oil as the displaced fluid were 0.875 and 6, respectively.The contact angle of water droplet on glass sheet in the oil-saturated media as the second fluid was 55°.A constant velocity boundary scheme proposed by Zou and He was used at both left inlet and right outlet of the model 55 and the upper and lower and grains' solid wall was half-way bounce back boundary condition.After converting these parameters from SI units to the lattice units, the simulation was performed with a capillary number of 0.063.The results (Fig. 11) showed that the fluid front pattern is fingering and some bubbles of saturated fluid remain in the pores.

Conclusions
The fluid front in EOR process of oil resources is important.Along with the experimental results, the use of accurate simulations, which are consistent with experimental results, will help petroleum engineers to assess petroleum reservoirs.Shan--Chen model of LBM is one of the mesoscopic-scale simulation tools by which we can study the fluid front in the pore scale.
In this study, the parameters of velocity, dynamic viscosity, with which the capillary number is defined, and wettability were investigated in formation of viscous fingering in a microfracture and a F i g .9 -Wettability effects on fluid front profile with injection velocity of 0.021, viscosity ratio of 1/3, G of 2 and the contact angles of CA = 30 (1), CA = 60 (2), CA = 90 (3), CA = 120 (4), CA = 150 (5) degrees in a microfracture (a) in the intervals of 500 time steps and in porous media (b) in breakthrough time simple porous media.Independently, the increase in velocity and decrease in dynamic viscosity ratio resulted in the formation of a fingering pattern, while wettability reduced it.In addition, the least amount of sweep efficiency occured when the displacing fluid had neutral wettability.
The results of this simulation show the strength and accuracy of the Shan-Chen Model of LBM in simulating the multi-component immiscible displacement, EOR, and fluid front tracking in the pores of porous media.
) s and s are different fluid phases.The amount of G ss represents the fluid-fluid particle interaction or the amount of attraction or repulsive force between fluid particles.The zero value indicates the miscibility of the phases, and its increase causes significant reduction in their miscibility.s G s represents fluid-solid interaction and determines the wettability.The effective density, which is introduced by Shan and Chen is as follows 17 .

F i g . 2 -
Contact angle F i g .3 -Simulation results of contact angle determination with different adhesion values are negligible.The shape of the fluid front is flat with few irregularities and fluid droplets are trapped in pore space

F i g . 5 -
Fluid flow front and viscous fingering in a microfracture F i g .6 -Fluid regimes and fluid front48