Modeling Immiscible Fluid Displacement in a Porous Medium Using Lattice Boltzmann Method

: The numerical investigation of the interpenetrating ﬂow dynamics of a gas injected into a homogeneous porous media saturated with liquid is presented. The analysis is undertaken as a function of the inlet velocity, liquid–gas viscosity ratio (D) and physical properties of the porous medium, such as porous geometry and surface wettability. The study aims to improve understanding of the interaction between the physical parameters involved in complex multiphase ﬂow in porous media (e.g., CO 2 sequestration in aquifers). The numerical simulation of a gaseous phase being introduced through a 2D porous medium constructed using seven staggered columns of either circular- or square-shaped micro-obstacles mimicking the solid walls of the pores is performed using the multiphase Lattice Boltzmann Method (LBM). The gas–liquid ﬁngering phenomenon is triggered by a small geometrical asymmetry deliberately introduced in the ﬁrst column of obstacles. Our study shows that the amount of gas penetration into the porous medium depends on surface wettability and on a set of parameters such as capillary number ( C a ), liquid–gas viscosity ratio (D), pore geometry and surface wettability. The results demonstrate that increasing the capillary number and the surface wettability leads to an increase in the effective gas penetration rate, disregarding porous medium conﬁguration, while increasing the viscosity ratio decreases the penetration rate, again disregarding porous medium conﬁguration.


Introduction
The increase in greenhouse gas (GHG) concentration in the atmosphere over the past decades has generated strong interest in developing technologies that help in reducing CO 2 emissions, especially those generated by fossil fuel combustion in heat and power generation [1][2][3]. Carbon dioxide sequestration, a fundamental part of a larger Carbon Capture and Storage (CCS) initiative [4], is an emerging field of research that is viewed as one of the potential solutions to decrease CO 2 concentration in the atmosphere. CCS is based on two main stages [5,6]. Firstly, it requires the capture of carbon dioxide from industrial and energy sources, and, secondly, it needs to transport and store it for longterm sequestration in underground aquifers, which are essentially porous media saturated with saline water isolated from the atmosphere. The penetration of a displacing fluid into a porous medium depends on the capillary and viscous forces at work, as well as the physical properties of the porous medium itself [7]. The displacement of the nonwetting fluid by the wetting phase is essential and affects CO 2 sequestration in terms of storage efficiency, capacity, and security. However, it is a rather complex phenomenon, difficult to study in laboratory or through purely analytical means. Numerical modeling is a complementary and powerful approach that can greatly help in understanding this phenomenon, provided that the discretized model is capable of capturing the rich physics involved in such a process.
Despite the continuous increase in computational power, modelling multiphase flow in a porous medium is challenging for conventional Computational Fluid Dynamics (CFD) tools due to the complex geometries involved, including small length-scale pores, and, hence, the required use of a computationally prohibitive adapted mesh [8]. For such meso-scale phenomena, the multiphase Lattice Boltzmann Method (LBM) has proved to be a promising alternative to standard CFD techniques [8,9]. Among these techniques, the volume-of-fluid (VOF) [10], the level set (LS) [11] and the phase field (PF) [12] methods are most common. However, according to Yang and Boek [7], these techniques have numerical instabilities at small-scale interface regions.
The LBM reproduces the macroscopic continuous dynamics of fluids, starting from the statistical approach embedded in the Kinetic Theory, in which the microscopic, discrete nature of a fluid can be described by particle distribution functions while avoiding the necessity to track individual molecules [9][10][11]. Additionally, in comparison to traditional CFD tools, the LBM algorithm is intrinsically parallel and effective in dealing with complex boundary geometries [13,14]. Since the very early stages of LBM, several multiphase models have been proposed, namely the color-fluid model, the pseudo-potential model, the free energy model and the mean-field theory model [15][16][17]. Despite the many advantages of LBM, all these multicomponent multiphase models are confined to use for small liquid-gas density ratios due to the presence of spurious currents at the interface, which ultimately promote numerical instabilities [18][19][20][21]. According to Fakhari and Rahiman [22], LBM multicomponent multiphase models have a limited range of gas-liquid density ratios (up to 15) due to these numerical instabilities. Generally, if the magnitude of the unphysical currents corresponds to the order of the local flow velocity, the model generates meaningless results [14].
There have been several attempts to improve the performance of multiphase LBM at high density ratios. For example, Inamuro et al. [23] proposed using a lattice kinetic scheme, which applies a single equilibrium distribution function calibrating the speed of sound, while Lee et al. [24] introduced a new method, which is based on a pressure equation and can be applied to a wide range of density ratios. Additionally, in the approach by Lee et al., spurious velocities at the interface are drastically reduced by increasing the surface thickness. Moreover, Yuan and Schaefer [25] studied the stability of LBM at high density ratios under various equation of states (EOS). It was found that Peng-Robinson (P-R) and Carnahan-Starling (C-S) EOS can model high-density ratios while minimizing the effects of spurious velocities because these EOS allow LBM to correctly capture the fluid flow phenomena in terms of temperature, in addition to other relevant thermodynamic parameters. Furthermore, Fakhari et al. [26] have found that the flow through the pore domain is extremely sensitive to the choice and accuracy of the boundary conditions, especially the entrance conditions.
Well-known numerical studies on porous media apply multiphase color-fluid LBM [6,17,27,28] with a density ratio of 1. However, the pseudo-potential model has shown better performance in terms of stability for fluid flows involving higher density contrasts, as, for example, in the case of CO 2 flow in porous media [25]. According to Li et al. [29], Kupershtokh et al. [30], Chen et al. [31] and Yuan-Schaefer [25], the performance of the model depends on EOS and the forcing scheme. There are different forcing schemes, such as velocity-shift, exact difference method (EDM) and Guo's forcing scheme, depending on how interaction force is introduced into LBM. The accuracy and stability of these forcing schemes have been analyzed in several studies [16,[29][30][31]. According to Li et al. [29], EDM is more stable when relaxation time is less than 1, i.e., τ < 1, while a velocity-shift forcing scheme is more stable when τ ≈ 1. Li et al. [29] also analyzed the effects of different values of τ on the coexistence curves. However, in a previous study, Kupershtokh et al. [30] showed that at τ ≈ 1, EDM and velocity-shift forcing schemes are identical. Due to the fact that LBM is still a newly emerging method, there are contradictions between some studies. Therefore, in this study we compare these two popular forcing schemes in terms of numerical stability and accuracy and make recommendations regarding multiphase flows in complex geometries, where a model's high stability is crucial.
We use P-R EOS in this study, as it is one of the accurate EOS, which helps specify the fluid thermodynamic properties through their acentric factor. We study the displacement process and its efficiency by applying pseudo-potential LBM instead of the widely used color fluid model. Throughout our work, the gas is injected with a uniform inlet profile into a saturated aquifer at different viscosity ratios, capillary numbers and pore geometry. By changing these parameters, their effects on displacement efficiency are investigated. The liquid-gas density ratio is fixed at 5 to ensure negligible spurious velocities at the interface. This apparently low-density ratio is also consistent with carbon dioxide injection condition in deep saline aquifers, which typically occurs at critical or near-critical regimes.
The paper is arranged as follows: The multiphase LBM with the P-R EOS model is proposed in Section 2. Section 3 is dedicated to validating the model, and Section 4 addresses extensive analysis using different fluid properties and pore configuration. The conclusions are presented in Section 5.

Lattice Boltzmann Method
In this study, the standard two-dimensional Lattice Boltzmann Method (LBM) with a Bhatnagar-Gross-Krook (BGK) collision operator is used. The Lattice Boltzmann Equation (LBE), which describes the rate of change of the particle distribution function in a generic point x on a Cartesian grid at a generic time t on a single direction of a predefined lattice, is given by Equation (1): where f i is the particle velocity distribution function associated with i-th direction of the lattice, ∆t is the time increment and e i is the lattice speed along the i-th direction. ∆ f i is the body force term. τ expresses the dimensionless relaxation time, which is used to describe the kinematic viscosity v = c 2 s (τ − 0.5)∆t. M is the number of discrete velocities and c s is the speed of sound. The number and directions of the set of M speeds define the lattice, which is typically referred to as DnQM, n being space dimensionality. Since the goal is to obtain the Navier-Stokes (NS) equations, a given lattice is acceptable if and only if it possesses enough symmetries to recover 4th-order isotropy of stress-strain tensors [8]. f eq i is the corresponding equilibrium distribution function described as follows: where u is the macroscopic velocity and w i is the weighting factor [25]. The form of Equation (2) is designed to recover the Euler equation. In this two-dimensional study, we adopt the widely used D2Q9 lattice, which is characterized by the following weighting factors (see Table 1): These weighting factors are necessary to recover isotropy for lattice models with different velocity vector modules. Macroscopic quantities like density ρ and momentum density can be recovered as discrete moments of the distribution functions, which are given by Equations (3) and (4), respectively. Additionally, in LBM, the incompressible NS equations with pressure are recovered via simple EOS rather than by solving a Poisson equation, as in standard CFD [32]. The ideal EOS with standard LBE is described in Equation (5).
The LBE described so far accounts for single-phase problems; to deal with multiphase flows, non-ideal effects must be introduced. In this study, we chose the pseudo-potential, or Shan-Chen [32], model, which consists of introducing a body force depending on the discrete gradient of a function of density, usually referred to as effective mass (Equation (6)).
where G(x,x) is the Green's function, which describes the interaction strength between neighbour particles. Depending on its magnitude, the interaction force can be described as repulsive or attractive. ψ(x) is the effective mass, which controls the EOS and, thus, the nature of the non-ideal fluid to be modeled. The SC model leads to the general EOS given by: The adhesion force between two fluids is expressed by Equation (8): where G ads is the interaction force between the fluid and solid walls. The interaction can be considered as non-wetting or wetting depending on its magnitude. Specifically, it is considered positive for a non-wetting fluid surface and negative for a wetting fluid surface [33]. To calculate the contact angle, the modified Young's equation is used as proposed by Huang et al. [34]: where G c is cohesion parameter, G ads is calculated from Equation (8) and is the density factor.
The Shan-Chen EOS is widely used for multiphase LBM simulations due to its simplicity (Equation (7)). However, in this study, the Peng-Robinson EOS [35] (Equations (10) and (11)) is preferred because it has proven to be capable of dealing with higher density ratios with a reduction in the instabilities associated with spurious currents [25].
where R is the universal gas constant, a (a = 0.45724R 2 T 2 c /p c ) and b (b = 0.0778RT c /p c ) are species-dependent coefficients, α denotes a function of the reduced temperature (T r = T T c , T c is the critical temperature of the fluid in consideration) and ω is the so-called acentric factor. In the present work, the values of a, b and R are set as 2/49, 2/21 and 1, respectively; thus, the critical properties of the LBM fluid are evaluated correspondingly.

Forcing Scheme Analysis
According to Kuperstokh et al. [30], the exact difference method (EDM) considers additional body forces to avoid discrepancies of distribution function at the higher-order terms of conservation laws. EDM can be derived from LBM by discretization of LBE in the velocity space. Finally, the body force term of the distribution function can be written as: where, ∆u = F∆t/ρ is change of mass velocity at a node. ∆ f i is equal to the difference in equilibrium distribution functions during time-step ∆t at constant density ρ.
The performance of the pseudo-potential model with both a velocity-shift forcing scheme and an EDM scheme was investigated at different relaxation times. Figure 1 shows the lowest achievable reduced temperature (T/T c ) for both schemes at τ < 1. Results from these tests show that the EDM forcing scheme is more stable. The stability of the velocityshift scheme rises with increasing τ. When τ ≈ 1, the difference between the stability of the two forcing schemes becomes not so significant, in correspondence with the results of Kupershtokh et al. [30].  Table 2 also compares the maximum values of spurious currents at various temperatures predicted by LBM with both forcing schemes. Overall, some little differences are noticed in the magnitudes of the spurious currents. These results suggest that numerical stability improves, allowing to reach lower temperatures (T/T c ≈ 0.79), and that the magnitude of spurious currents decays as relaxation time τ approaches 1.

Validation of the Model
The validation of our model was carried out in consideration of a capillary pressure test with a known analytical solution. The test considered the injection of a gaseous phase from a plenum into two parallel passages (see Figure 2). The inlet and outlet boundary conditions were of uniform velocity and constant pressure, respectively, while the top and bottom walls were prescribed as on-grid bounce-back. The passages had different widths (r 1 Fluids 2021, 6, 89 6 of 14 and r 2 ) and, correspondingly, different capillary pressures (P c1 and P c2 ). Capillary pressure can be calculated by the Young-Laplace equation where σ is the surface tension. Then, the obtained values are compared to the pressure difference (∆p) at the inlet (p in ) and the outlet (p out ). Depending on the conditions present, the fluid can either enter or not enter the passages (note that passage 1 is narrower than passage 2). When ∆p is less than P c2 , the gas cannot enter any of the passages, while if ∆p is in the range from P c1 to P c2 , the gaseous phase can only penetrate through passage 2.
The computational domain was discretized using 280 lattice units (lu) in the x-direction and 140 lu in the y-direction. The passages 1 and 2 are l00 lu in length by 15 lu-width and 30 lu-width, respectively. Figure 3 shows the excellent agreement between the numerical results and the expected flow behavior.

Results and Discussion
The schematic setup of the porous domain is shown in Figure 4a,b. Inlet and outlet boundary conditions are uniform velocity and constant pressure, respectively, while the top and bottom walls are prescribed as periodic boundary conditions. The computational domain has a size of 401 × 401 (lu 2 ) and is filled with staggered columns of circular pores with a radius of 10 lu or square pores at 20 × 20 (lu 2 ), separated by 10 lu. Each column of pores is spatially arranged uniformly with just a small perturbation (+1 and −1 lu) among contiguous pores to facilitate the break-up of flow symmetry during the simulation and resemble a more realistic flow behavior. Additionally, one pore is removed from the 1st column to observe its effect on the fingering phenomenon. The uniform size and position of the pores in our study allows us a more precise control of the problem, which is different to what happens in experiments where pores often have different positions, pore sizes or variations in wettability. Initially, the porous medium is filled with the liquid phase, while the gaseous phase is introduced from the left to the right at a constant and uniform velocity. All simulations are conducted at a liquid-gas density ratio of 5, in order to have negligible spurious currents at gas-liquid interfaces [25].
The penetration of the gaseous phase is expected to depend on multiple dimensionless parameters, namely the capillary number (C a ), the viscosity ratio (D) and surface wettability. The capillary number C a describes the relationship between viscous and interfacial forces: where u G and n G are the average velocity and the dynamic viscosity of the gaseous phase, respectively. The viscosity ratio is defined as the ratio between the viscosity of the liquid and gaseous phases: where n L is the dynamic viscosity of the liquid phase. Furthermore, a common dimensionless time t * defined in Equation (16) is used in the comparison of gas penetration for cases with different injection velocity.
where t is the total time-steps and H is the width of the domain. As a first stage in the modeling process, we performed a forcing scheme analysis to ensure its proper selection. For this purpose, the model with surface contact angle θ = 70 • , C a = 0.038 and circular pore geometry was selected. In this test, two different forcing schemes (namely velocity-shift and EDM) were analyzed under the same conditions at τ = 1. As can be seen from Figure 5, all simulations have approximately the same results, as previously proven in Section 2.2. Thus, further calculations at τ = 1 will be presented using the original velocity-shift forcing scheme. Nevertheless, the authors would like to suggest that for low-temperature and low-viscosity ratio fluid flow conditions not within the scope of this paper, the advantages of EDM over velocity-shift should be scrutinized, since the former might possess better numerical stability. Multiple simulations were performed to identify the effect of C a on the penetration of the gaseous phase into the porous media. Initially, the viscosity ratio was equal to 1 and the surface contact angle was set to θ = 70 • , indicating a wettable condition. Three different capillary numbers were examined in this study, namely C a1 = 0.038, C a2 = 0.076 and C a3 = 0.115. Different C a can be obtained by varying velocity at the inlet boundary. Figures 6 and 7 illustrate the different penetration lengths at the common nondimensional time t * = 0.175 at which the same amount of gas has been injected despite the injection velocity. Figures 6 and 7 show that a larger C a favors faster penetration of the gaseous phase into the porous media, as it was expected. Additionally, the fingering can be observed for each selected array of pore geometry near the entrance of the domain. Furthermore, at a lower C a (Figures 6a and 7a), it can be observed that the fingering favored the span-wise direction (with a larger pore throat) due to lower entry pressure, while it moved towards the stream-wise (front) faster with increasing C a . Nevertheless, backward fingering, also known as capillary fingering, is hardly appreciated in our simulations, despite previous works observing that phenomenon being driven by the capillary force [28]. Lenormand et al. [36] also pointed out that capillary fingering may occur through a large pore throat in backward direction due to low entry pressure, but this backward-fingering is more noticeable for non-homogeneous porous media where the difference in pore throat between the cylinders is higher than that of the spacing used in the present investigation. pore throat in backward direction due to low entry pressure, but this backward-fingering is more noticeable for non-homogeneous porous media where the difference in pore throat between the cylinders is higher than that of the spacing used in the present investigation.  It is evident that, depending on pore geometry, the effective gas penetration, measured by length T (Figure 8), changes slightly, as shown in Table 3. Figure 8 illustrates the domain length (L), the slip distance (S) and the effective penetration length (T). In our case the slip distance is the same for all simulations. It is evident that, depending on pore geometry, the effective gas penetration, measured by length T (Figure 8), changes slightly, as shown in Table 3. Figure 8 illustrates the domain length (L), the slip distance (S) and the effective penetration length (T). In our case the slip distance is the same for all simulations.   Table 3 shows that, for low C a , gas penetration is almost the same for both geometries, but as C a increases, the penetration through the square-shaped pores is augmented significantly more than for the circular pores. We noted that extremely low C a promotes the growth of spurious currents associated with numerical instabilities, consistent with findings by Raeini et al. [37].
Additionally, our results overcame the limitations found by Lenormand et al. [36], where it was stated that in a heterogeneous pore domain, a high C a will force the fluid to occupy a limited fraction of the domain. Instead, our model has periodic boundary conditions, uniform size and uniform positioning of pores, removing that pitfall.
Furthermore, the effect of the viscosity ratio D is examined by varying gaseous phase viscosity while keeping the viscosity of the liquid phase constant. Additionally, C a and θ are set as constant and equal to 0.038 and 70 • , respectively. Figures 9 and 10 show that increasing the viscosity ratio makes the gas injection more difficult and increases the time needed to displace the liquid phase out of the domain. The results corroborate previous studies, indicating that, at low D, the fingering narrows and spreads sideways faster, while it locally thickens when D increases [17]. Moreover, our results show the expected [27] easier front penetration at lower viscosity ratios, depicting how the intruding phase progresses faster in a finger-like manner as shown in Figures 9a and 10a. Our results also corroborated that increasing the viscosity ratio favors a uniform frontal displacement (Figures 9c and 10c). It should be noted that pore geometry affects the effective gas penetration length. According to Table 4, in the circular-shape pore domain, the gas penetration length is shorter than in the square-shape pore domain.   The effect of contact angle on gas penetration is investigated by changing the adhesion parameter between wall and fluid. As was previously discussed, depending on capillary pressure (Equation (13)), the gross evolution of the gaseous phase through the pores might be predicted. Consequently, the surface contact angle effectively affects the displacement of the gaseous phase into the liquid-filled porous domain. Previous studies [6,17] on porous media were done on a constant contact angle. Therefore, multiple simulations with various surface contact angles, namely θ = 83 • , θ = 70 • , θ = 50 • and θ = 33 • at fixed C a , were carried out. Figure 9 demonstrates the evolution of the gaseous phase into the porous domain at t * = 0.225. All simulations were performed for wetting conditions, i.e., θ < 90 • . Figure 9 shows how gas penetration length increases while the surface contact angle decreases, as was expected. However, the effective gas penetration in the circular-pore domain ( Figure 11) is higher than that in the square-pore domain ( Figure 12) due to the narrower inter-pore pas-sages (see Figure 13).

Conclusions
A multiphase LBM model was employed to study the influence of capillary number, viscosity ratio and surface contact angle on the hydrodynamics of gas penetration into a homogeneous 2D porous medium saturated with liquid. Previous studies [6,17,27,28] were done based on the phase-field (PF) method in a framework of a color-fluid multiphase LBM model with a density ratio of 1. In the present work, the Peng-Robinson EOS was successfully integrated into LBM in order to increase the stability of the model. As a result, P-R EOS enables the modeling of the gas penetration phenomenon with a liquid-gas density ratio of 5, resulting in negligible instabilities (i.e., spurious currents) at the interface. Additionally, the validity of the prescribed forcing scheme was assessed as one the main strategies to increase the stability of the model. It was found that at τ = 1 the velocityshift and EDM forcing schemes performed similarly in terms of stability, and, therefore, the former was chosen to develop the whole set of simulations in porous media.
The numerical results demonstrated that increasing the capillary number (C a ) and surface wettability enhanced the effective gas penetration at constant viscosity ratio D, while an opposite effect is noticed when increasing D whilst maintaining constant C a and wettability. Furthermore, the shape of the pores influences fingering length and effective gas penetration. In the present study, the effective gas penetration was higher in the squarepore domain than that in the circle-pore domain by approximately 10% and 4% for explored C a,max and D min , respectively. Moreover, for the largest wettability explored in this study, the effective gas penetration in the square-pore domain was less than 3% than that in the circle-pore domain. Future work is planned to investigate gas-liquid penetration behavior in heterogeneous pore domains while increasing the viscosity ratio. What is more, in future study, we intend to use the EDM forcing scheme for better stability at low relaxation times.

Data Availability Statement:
The data presented in this study are available in the displayed figures and tables within the article, but can also be obtained on request from the corresponding author within a reasonable time-frame.