A Pseudopotential Lattice Boltzmann Method for Simulation of Two-Phase Flow Transport in Porous Medium at High-Density and High–Viscosity Ratios

In this work, the capability of a multiphase lattice Boltzmann method (LBM) based on the pseudopotential Shan-Chen (S-C) model is investigated for simulation of two-phase flows through porous media at high-density and high–viscosity ratios. The accuracy and robustness of the S-C LBM are examined by the implementation of the single relaxation time (SRT) and multiple relaxation time (MRT) collision operators with integrating the forcing schemes of the shifted velocity method (SVM) and the exact difference method (EDM). Herein, two equations of state (EoS), namely, the standard Shan-Chen (SC) EoS and Carnahan-Starling (CS) EoS, are implemented to assay the performance of the numerical technique employed for simulation of two-phase flows at highdensity ratios. An appropriate modification in the forcing schemes is also used to remove the thermodynamic inconsistency in the simulation of two-phase flow problems studied at low reduced temperatures. The comparative study of these improvements of the S-C LBM is performed by considering an equilibrium state of a droplet suspended in the vapor phase. The solver is validated against the analytical coexistence curve for the liquid-vapor system and the surface tension estimation from the Laplace Law. Then, according to the results obtained, a conclusion has been made to choose an efficient numerical algorithm, including an appropriate collision operator, a realistic EoS, and an accurate forcing scheme, for simulation of multiphase flow transport through a porous medium. The patterns of two-phase flow transport through the porous medium are predicted using the present numerical scheme in different flow conditions defined by the capillary number and the dynamic viscosity ratio. The results obtained for the nonwetting phase saturation, penetration structure of the invading fluid, and the displacement patterns of two-phase flow in the porous medium are comparable with those reported in the literature. The present study demonstrates that the S-C LBM with employing the MRT-EDM scheme, CS EoS, and the modified forcing scheme is efficient and accurate for estimation of the two-phase flow characteristics through the porous medium.


Introduction
The lattice Boltzmann method (LBM) is known as a powerful mesoscopic numerical scheme for simulation of the multiphase flows transport through the porous medium [1][2][3][4]. The kinetic nature of the LBM allows to model the phase separation, interparticle interactions, and interfacial phenomena in multiphase flows [5,6]. Due to these favorable features, numerical solution methods based on the multiphase LBM do not need an additional algorithm for tracking or capturing the interfacial dynamics, in contrast with the Navier-Stokesbased flow solvers. The efficiency and remarkable versatility of multiphase LBM make this method widely used for simulation of fluid flows in the applications involving interfacial dynamics in porous geometry.
Several multiphase models developed in the LBM framework [7][8][9][10][11]. Among the existing models, the pseudopotential Shan-Chen (S-C) model [8,9] is still attractive due to its capability for simulation of multicomponent multiphase flows and simplicity in implementation with appropriate accuracy and performance [12]. The S-C model has the ability to compute the phase transition [8] and predict properties of multicomponent multiphase systems [9,13,14]. Despite the advantages of the S-C model, this scheme suffers from spurious currents in the interfacial region of two-phase flows at high-density ratios [15]. The unphysical dependency of the magnitude of the density ratio to the value of the viscosity ratio is another restriction of the standard S-C model for the simulation of multiphase flows [16]. To apply the S-C LBM for the simulation of multiphase flow transport through porous media, it needs to review the improvements made in this scheme to eliminate the mentioned restrictions.
Several techniques have been proposed to improve the stability and accuracy of the LBM incorporated with the pseudopotential S-C model for studying multiphase flows at high-density and -viscosity ratios. Implementation of different realistic equations of state (EoS) [16][17][18] and using higher-order spatial and temporal discretization [19,20] are the attempts made for the numerical solution of multiphase flow problems at high-density ratios. Another category includes the extension of multiple relaxation time (MRT) technique for S-C LBM [21,22] to eliminate the numerical instability of this approach for simulation of flow problems with high-viscosity ratios. The unphysical dependency of the viscosity-ratio to the value of the density-ratio in the S-C LBM is also resolved by employing appropriate forcing schemes to compute the interparticle interaction forces [23,24]. More details of these improvements are provided below.
A review of the literature relevant to the S-C LBM shows that the implementation of EoSs, e.g., Peng-Robinson (P-R) and Carnahan-Starling (CS) EoSs instead of the standard Shan-Chen (SC) EoS, is an attractive and efficient technique to achieve higher density ratios by the elimination of spurious currents. Although, Ezzatneshan and Vaseghnia [18] have recently represented that considering wall boundary conditions with wettability effects impacts the performance of the S-C LBM incorporated with the EoSs. A thermodynamic inconsistency can also appear in the solutions based on the standard S-C model at high-density ratios which leads to a considerable deviation of numerical results from the theoretical data in the phase change diagram [25,26]. To solve this problem, Gong and Cheng [25] proposed a free parameter that is added into the interaction force computed based on the S-C model to control the strength of interparticle interactions. This technique is employed by Stiles et al. [27] to compute the coexistence curve for a realistic high-density ratio two-phase system, including the air and water. Recently, the method proposed by Gong and Cheng is also implemented for simulation of complex bubble dynamics in a cavitating flow problem [28].
In the numerical solution of single-phase fluid flows at high Reynolds numbers using the LBM, the MRT collision operator is introduced as an efficient alternative approach to the standard single relaxation time (SRT) technique [29,30]. In the same way, the MRT scheme is also implemented in the S-C LBM to eliminate the numerical instability of this model for simulation of two-phase flows with high-viscosity ratios. Kuzmin et al. [21] have investigated the advantages of dealing with the MRT collision operator in the study of capillary effects and the formation of the droplets. Yu and Fan [31] have shown that using the MRT-LBM leads to a considerable reduction in the spurious currents in the interfacial region of two-phase flows. Recently, the MRT S-C LBM is also implemented for studying various flow problems, e.g., the simulation of two-phase flow through porous media [32,33], investigation of cavitation bubble dynamics [18,34], and studying the heat transfer in multiphase and multicomponent fluid flows [34,35]. Evaluation of the performance of the S-C LBM with the implementation of the MRT collision operator for simulation of multiphase flows is an interesting subject which this article deals within different flow conditions.
In addition to the aforementioned improvements made for the development of high-density and -viscosity ratio S-C model, some unphysical properties of this model are also resolved by incorporating appropriate forcing schemes. The calculation method of the interparticle interaction force in the pseudopotential LBM has a significant effect on the accuracy and performance of the numerical solution of multiphase flows. A main restriction of the standard S-C LBM in this category is the unphysical dependency of the magnitudes of the density ratio and viscosity ratio of a multiphase system which can be eliminated by the implementation of forcing schemes. There are three well-known forcing schemes that are widely used by the S-C model. The standard Shan-Chen forcing scheme is known as shifted velocity method (SVM) [8,9], the Guo et al. method [23], and the exact difference method (EDM) [24]. Evaluation of these schemes for the simulation of multiphase flows by using the S-C LBM shows that the EDM is numerically stable in a wide range of density ratios in comparison with the SVM and Guo schemes [16]. It has also been shown that the EDM forcing scheme gives more accurate predictions for high-density ratio multiphase flow characteristics with employing the Carnahan-Starling EoS [36]. Recently, Peng et al. [37] have studied the performance of the SVM, EDM, and Guo schemes for the simulation of the phase separation phenomenon in a two-phase flow system. They found that all these schemes are consistent with the Laplace law; however, the maximum density ratio which can be achieved by the EDM scheme is greater than that of obtained by implementation of other forcing schemes.
This wide range of improvements shows that it is crucial to understand the capability and performance of the modified techniques of the S-C LBM, including the EoSs, collision operators, and forcing schemes, to employ an efficient algorithm for simulation of multiphase flow transport through a porous medium. However, to the best knowledge of the authors of this article, there is not a thorough comparative study in the literature to simultaneously evaluate the SRT and MRT collision operators incorporated with appropriate forcing schemes and realistic EoSs for the numerical solution of two-phase flows using the pseudopotential LBM. In this paper, the accuracy and robustness of the S-C LBM are investigated by the implementation of the SRT and MRT collision operators with integrating the SVM and EDM forcing 2 Geofluids schemes. The capability of the standard SC EoS for simulation of two-phase flows at a high-density ratio is examined in comparison with the CS EoS. In the present study, the method proposed by Gong and Cheng is also used to remove the thermodynamic inconsistency in the simulation of twophase flow problems studied. The comparative study of these improvements of the S-C LBM is performed by considering an equilibrium state of a droplet suspended in the vapor phase. The numerical technique is validated against the analytical coexistence curve for the liquid-vapor system and the surface tension estimation from the Laplace Law. Then, according to the results obtained, a conclusion has been made to choose an efficient numerical algorithm, including an appropriate collision operator, a realistic EoS, and an accurate forcing scheme, for simulation of multiphase flow transport through the porous medium using the S-C LBM employed. Finally, the characteristics of two-phase flow transport through a porous medium are computed based on the present numerical scheme in different flow conditions defined by the capillary number and the dynamic viscosity ratio. The performance of the implemented numerical method based on the S-C LBM is also investigated by predicting the parameters of a flow transport through the porous medium at high-density and -viscosity ratios. The present paper is organized as follows: In Section 2, the LBM with SRT and MRT collision operators is briefly presented. The pseudopotential multiphase LBM coupled with the EoSs and the incorporating of forcing schemes for computing the interparticle interaction forces are introduced in this section. The numerical results obtained for different liquid-vapor flow problems and the relevant discussions are presented in Section 3. Finally, some conclusions are summarized in Section 4.

Multiphase Lattice Boltzmann Method
In this section, the governing equations of LBM with SRT and MRT collision operators are introduced. Then, the multiphase pseudopotential S-C model [8,9] with the standard formulation and its improved form with the realistic CS EoS is presented. Two forcing schemes, including SVM and EDM methods, are outlined for the incorporation of interparticle interaction force in the numerical solutions. The modification proposed by Gong and Cheng [25] is also reviewed to eliminate the thermodynamic inconsistency in the present computations.
2.1. SRT-and MRT-LBM Equations. The governing equation of the standard single relaxation time (SRT) LBM is given as [38]: where f a ðx, tÞ is the particle distribution function which demonstrates the distribution of particles with the velocity e a in position x at time t. The lattice speed c = Δx/Δt is the ratio of the lattice space Δx to the time step Δt, and τ is the relaxation time which is related to the kinematic viscosity ν = c 2 s ðτ − 0:5ΔtÞ [3]. The subscript a indicates the discrete velocity directions. According to the two-dimensional (2D) lattice structure D2Q9 employed in the present work, the set of discrete velocities e a is defined as The right-hand side of Eq. (1) represents the SRT collision operator which leads the system to the local Maxwellian equilibrium on the time scale τ. The f eq a ðx, tÞ is the equilibrium distribution function which is defined as The speed of sound c s for the LBM in D2Q9 framework is equal to 1/ ffiffi ffi 3 p , and the weighting coefficient w a is given by w α = 4/9 a = 0, 1/9 a = 1 − 4, 1/36 a = 5 − 8: The macroscopic characteristics of the fluid flow then can be computed as For the implementation of MRT collision operator in the LBM presented in Eq. (1), the left-hand side of the equation is rewritten as Herein, m a = f a * M and m eq a = f eq a * M are the distribution function and the equilibrium distribution function in moment space, respectively. The transformation matrix M and the relaxation matrix S can be defined as follows [39]: where the values of s 1 , s 4 , and s 6 should be nonzero in order to maintain the influence of the trapezoidal integration in Eq. (6) [29]. In the present work, these parameters set to be s 1 = s 4 = s 6 = 1, s 2 = s 3 = 0:8, s 5 = s 7 = 1:2, and s 8 = s 8 = 1/τ.

Pseudopotential Shan-Chen Model.
Simulation of the phase transition process and the interfacial dynamics in the multiphase flows need to consider the adhesive and cohesive forces between the fluid-fluid and solid-fluid particles. The S-C model is known as an efficient method for computing these interparticle interaction forces in the multiphase LBM framework. Herein, the standard S-C model proposed by Shan and Chen [8,9] and its improved form by incorporating the CS EoS are presented. The numerical technique proposed by Gong and Cheng [25] is also applied to remove the thermodynamic inconsistency in the present solutions.
In the standard S-C LBM, the interaction and adhesive forces are imposed by shifting the macroscopic velocity in the equilibrium distribution function f eq a ðx, tÞ as where the total force F total is a summation of the interaction force F int , adhesive force F w , and external body force F B : This form of incorporation of the forcing term into the S-C LBM is known as the standard shifted velocity method (SVM) which is discussed in Section 2.3.
Yuan and Schaefer [17] proposed the following relation for the interaction force F int : where c 0 is a constant parameter and equals to 6 for the D2Q9 lattice structure. Also, the interparticle interaction potential g indicates the attractiveness (repulsiveness) of the interparti-cle interaction force for g < 0 (g > 0). In the D2Q9 model, the ∇ψ can be discretized as below: where c 1 and c 2 are equal to 1/3 and 1/12, respectively. The mean-field potential function ψ defines the relevance between the interaction strength and the mass of particles.
In the S-C model, ψ depends on the local density ρ as [8].
The physical velocity of the multiphase flow is given by and the macroscopic pressure is expressed as which is known as the standard SC EoS in this paper. For this EoS, T c = 4:5ρ 0 and ρ 0 = 1. To incorporate a realistic EoS in the S-C model, Eq. (14) can be rewritten as Then, the pressure obtained from the EoS in the form of p = pðρÞ can be used for computing the function ψ. For the simulation of multiphase flows at high-density ratios by using the S-C LBM, a variety of EoSs is proposed in the literature [18]. According to the recent studies in the literature [16,18], it is shown that the CS EoS provides more realistic results, particularly when the two-phase flow problem includes a wall boundary condition with wettability effects. Also, it is reported that the CS EoS provides more numerical stability by reducing spurious currents in the interfacial region that leads to achieving higher density ratios. Herein, the noncubic CS EoS [40] is implemented to describe the relationship between the fluid properties by the following equation: In the CS EoS, a = 0:4963R 2 T c 2 /p c and b = 0:18727RT c /p c are the attraction and repulsion parameters. The critical temperature is equal to T c = 0:3773ða/ðbRÞÞ, where R is the gas constant. In the LBM framework, these parameters are set to be a = 1, b = 4, and R = 1 [17]. Accordingly, the critical density ρ c and critical temperature T c will be 0.13045 and 0.0943, respectively.

Geofluids
The S-C LBM coupled with any EoSs may lead to the thermodynamic inconsistency in the simulation of multiphase flows at high-density ratios [41]. In order to eliminate such a numerical problem, Gong and Cheng [25] have proposed a modification in computing the interparticle interaction force F int based on Eq. (10). This equation can be rewritten as follows: In terms of the numerical implementation, the F int computed by Eq. (10) includes both the local and the nearestneighbor lattice points, while the present form based on Eq. (17) only deals with the local lattice nodes. In the modification technique proposed by Gong and Cheng, F int is calculated by a combination of Eqs. (10) and (17) as where β is a weighting factor that controls the effectiveness of each relation in the calculation of the interaction force. The value of β depends on which EoS is employed in numerical solutions [25]. In this paper, the effect of parameter β is investigated on the accuracy and stability of the present numerical scheme based on the S-C LBM coupled with the CS EoS. Finally, the adhesive force F w between the two-phase fluid particles and solid nodes is expressed as where ψ w is calculated by considering a wall density ρ w for solid surface. The desired wettability characteristic of the surface defines the value of ρ w . A value of ρ w close to the liquid density ρ l corresponds to the hydrophilic wall boundary and consequently, the contact angle will be θ w < 90 ∘ . In contrast, for a ρ w value close to the vapor density ρ v , the solid surface acts as a hydrophobic surface with θ w > 90 ∘ . In Eq. (19), the indicator s is used to turn the adhesion force on for the solid nodes by s = 1 and turn it off on the fluid nodes by s = 0 in the numerical solutions performed.

Forcing Schemes.
The total force F total , including the interaction, adhesive, and body forces, is usually implemented via forcing schemes to the pseudopotential S-C LBM for simulation of multiphase flows. Therefore, employing a forcing scheme can affect the performance and numerical stability of the S-C model. In the present study, the standard SVM and the EDM approaches are considered to evaluate the effects of the forcing term incorporation into the S-C LBM.
In the SVM forcing scheme with both of the SRT and MRT collision operators, the F total is imposed on the velocity field through Eq. (8). Then, the equilibrium distribution function f eq a presented in Eq. (3) is computed by the u eq . However, in the EDM approach, the force term is directly added to the right-hand side of Eq. (1) as follows: To impose the forcing scheme based on the EDM with MRT collision operator, the force term needs to define in the moment space. Then, Eq. (20) is rewritten for MRT-LBM with EDM scheme as In the present study, the effect of the implementation of the EDM scheme for the incorporation of the force term in the S-C LBM is investigated for simulation of multiphase flows at high-density and -viscosity ratios. For the sake of clarity, the flowchart of the present solution based on the governing equations is presented in Figure 1.

Results and Discussions
In this section, the numerical algorithms presented based on the standard and modified S-C LBM are evaluated for the simulation of multiphase flows. For this comparison study, a combination of numerical schemes with SRT and MRT collision operators, SC and CS EoSs, SVM and EDM forcing schemes, and also the effect of control parameter β is considered. At first, the accuracy and performance of different algorithms applied based on the S-C LBM are investigated by the simulation of two benchmark two-phase flow problems, including the equilibrium state of a droplet suspended in the vapor phase and equilibrium state of a droplet on a solid wall with wettability effects. The results obtained based on the present solutions are compared with the analytical and available numerical results. According to this comparison study, a numerical scheme is proposed for simulation of two-phase flow problems at high-density and -viscosity ratios. Then, the efficiency of the proposed solution technique is validated by the simulation of two-phase flow transport through a duct. Finally, the numerical solution of the two-phase flow transport through a porous medium is performed by employing the proposed S-C LBM, and the flow characteristics are studied in various flow conditions.
3.1. Study of Accuracy and Performance. For a 2D liquid droplet suspended in the vapor phase, the surface tension σ causes a pressure difference Δp between the outside and inside of the droplet. The Laplace law [42] defines a linear relation between σ and Δp for a certain droplet radius R as In order to verify the present numerical results obtained based on the S-C LBM by the Laplace law, a stationary droplet placed in a fully periodic square domain is considered, and 5 Geofluids simulations are performed for various droplet radiuses. Herein, the droplet is initialized in the center of the flow field ðx 0 , y 0 Þ, and the density distribution is defined by Eq. (23) to preserve the stability of numerical solutions at the beginning iterations [12].
ð23Þ Figure 2 illustrates the results obtained for Δp versus 1/R by using the SRT-SVM, SRT-EDM, and MRT-EDM schemes coupled with the SC EoS. The linear relation between these two parameters in this figure shows that the present numerical results based on all the mentioned schemes satisfy the Laplace law. The slope of the line represents the value of the surface tension of the liquid-vapor system studied which is obtained σ = 0:05923, 0:06442, and 0:06292 by the SRT-SVM, SRT-EDM, and MRT-EDM schemes, respectively. The satisfaction of the Laplace law in this study demonstrates the accuracy of the employed S-C LBM coupled with different forcing schemes for predicting the characteristics of a stationary two-phase flow system. The slight difference between the values obtained for the surface tension depends on the efficiency and accuracy of these schemes which is discussed below with more details.
The accuracy of the present numerical schemes based on the S-C LBM is evaluated by using SRT and MRT collision operators, SC and CS EoSs, and SVM and EDM forcing schemes with the simulation of the stationary droplet at high-density and -viscosity ratios. Figure 3 demonstrates the density of the coexisting phases calculated in the equilibrium state of the droplet at different values of reduced temperature T/T c with relaxation time τ = 1:0. It should be noted that the parameter g defines the temperature effects in the SC EoS by g = −1/T, which gives T/T c = −4:5ρ 0 /g. Herein, the analytical solution based on the theory of Maxwell equal-area construction [43] is used to verify the present coexistence curves. Also, to consider the effect of the β -modification to remove the thermodynamic inconsistency in computing the interaction force by Eq. (18), this parameter is set to be β = 1 for the standard calculations and β ≠ 1 for the modified S-C LBM.
The numerical results presented in Figure 3 indicate that all of the S-C LBM schemes applied based on both the SC and CS EOSs accurately predicted the density of the liquid phase in comparison with the analytical solution. Also, at higher reduced temperatures correspond to the low-density ratios, the numerical results for the vapor phase are in good agreement with the analytical solution. However, the density values computed for the vapor phase at higher density ratios (low reduced temperatures) have a significant deviation from the analytical solution by all the schemes coupled with the standard interparticle interaction force (β = 1). The x − axis of Figure 3 is plotted in a log scale to have a better view of the difference between the results. It is obvious that below the reduced temperature T/T c = 0:85, the deviation between the numerical results and the analytical solution becomes substantial. Due to the higher value of the vapor density predicted by the SC EOS, the difference between the vapor density of this model with the analytical solution does not a significant impact on the value of the density ratio of the whole two-phase flow system. However, for the CS EOS, a slight deviation of the vapor density from the analytical solution leads to a considerable variation of the flow density ratio because of very small values computed for the density of the vapor phase by this EOS. Therefore, it is crucial to improve 6 Geofluids the consistency of these EOSs with the analytical solutions for accurate simulation of two-phase flow systems at highdensity ratios. As mentioned before, the β-modification is considered in the present work to resolve this problem. As shown in Figure 3, the results confirm that the improved interaction forcing scheme by the β parameter has a significant effect on the accuracy of the SC and CS EOSs. In this case, the appropriate value obtained for β by CS EOS is between 1.15 and 1.18 depending on the SRT and MRT collision operators employed. Although, there is no considerable difference between the numerical results obtained based on these two β values. For the SC EoS, the value of β = 0:886 is employed. Using different values for β according to the type of EOSs is also reported by Ping and Cheng [25]. By employing the β-modification, the density values obtained by the present numerical simulations for the vapor phase is in good agreement with those obtained by the analytical solution in a wide range of reduced temperatures. Although, the results obtained by the SC EOS still demonstrate some deviation below the T/T c ≃ 0:55 which can be due to the strong spurious currents in the interfacial region at such density ratios [44]. It is worth to mention that the SRT and MRT collision operators have not significant effect on the accuracy of the results obtained for the density ratio in comparison with the analytical solution. However, the MRT collision operator influences the numerical stability by decreasing the spurious currents at low reduced temperatures.
For more detail discussions regarding to the unphysical currents, the magnitude of maximum spurious velocity jU s j max obtained by the standard and modified interaction forcing schemes based on the SC and CS EOSs is presented in Figure 4 for the equilibrium state of a droplet at various reduced temperatures. This figure indicates that with employing both the SC and CS EoSs, the SRT-SVM and SRT-EDM models lead to almost the same values for jU s j max with the implementation of both the standard and improved interaction forces. However, with applying the MRT collision operator, the computed values for jU s j max are lower than those of obtained by SRT collision operator for all the schemes at various reduced temperatures. Increment of the magnitude of the maximum spurious velocity with decreasing the reduced temperature is obvious in Figure 4. The results presented in this figure also show that by applying CS EoS for simulation of the droplet at high reduced temperatures (low-density ratios), the parameter β has no remarkable effect on the magnitude of spurious currents. In contrast, at high-density ratios, the jU s j max is significantly decreased by using the improved interaction force, particularly with employing the MRT-EDM model. This result demonstrates that the CS EoS is more realistic and efficient than SC EoS for simulation of two-phase flows at highdensity ratios due to producing less spurious currents at the interfacial region.
As described in relation ν = c 2 s ðτ − 0:5Þ, the viscosity value is determined by the value of relaxation time τ in numerical solutions based on the LBM. Then, τ ⟶ 0:5 corresponds to lower viscosity magnitudes which can be used for simulation of two-phase flows at high Reynolds numbers. However, in flow simulations based on the standard LBM, serious numerical instabilities arise at the values of relaxation time close to 0:5 [45]. It is expected that this numerical limitation will be intensified when the two-phase flow simulations are performed at high-density ratios. To study the mentioned restriction and examine the capability of the present numerical method based on the S-C LBM with SRT and MRT collision operators, SC and CS EoSs and SVM and EDM forcing schemes for the simulation of two-phase flows at low viscosity magnitudes, the simulation of the stationary droplet is considered at τ = 0:52. This value for τ is the minimum relaxation time that can be achieved a stable numerical solution of the stationary droplet at reduced temperature T/ T c = 0:5 by the present MRT-EDM model and CS EoS. Figure 5 presents the coexistence curve for the stationary droplet simulations at τ = 0:52 by both the SC and CS EOSs and implementation of the numerical schemes employed. At the low viscosity magnitude considered (τ = 0:52), it is obvious in Figure 5 that the MRT-EDM scheme is more efficient than SRT-SVM and SRT-EDM techniques to achieve lower reduced temperatures that lead to the simulation of the two-phase flow at high-density ratios. The obtained results show that the minimum achievable reduced temperature is T/T c = 0:6 by applying the SC EOS with the MRT-EDM scheme. However, by employing the CS EOS with MRT-EDM and modified interaction force, this magnitude is T/T c = 0:475, that consequents a density ratio of 1597.8. To study the numerical stability of the present methods, Figure 6 illustrates the spurious currents calculated in terms of the reduced temperature. It can be seen that the MRT-EDM with the improved forcing scheme produces minimum spurious currents among the other methods. By considering the results presented in Figures 5 and 6, the advantages of employing the MRT  7 Geofluids scheme with β − modification of the forcing scheme are obvious for the simulation of a two-phase flow system at low relaxation times (τ ⟶ 0:5).
It should be noted that using the EDM scheme in the present study eliminates the unphysical dependency of the magnitudes of the density ratio and viscosity ratio of the two-phase flow system considered. To demonstrate this feature, simulation of a droplet in the equilibrium state is performed by the S-C LBM with SRT-SVM, SRT-EDM, and MRT-EDM schemes at T/T c = 0:8. Figure 7 shows the variation of the density ratio of each model at different relaxation times. As can be seen, by using the SVM scheme, the density ratio is decreased by increasing the magnitude of the relaxation time at a defined reduced temperature, which is unphysical. However, by using the EDM scheme, the density ratio remains constant as the reduced temperature T/ T c changes. The notable difference between the MRT-EDM and SRT-EDM schemes is the minimum achievable value of the relaxation time for the simulation of the twophase flow by using these techniques. As given in Table 1, the present study shows that for simulation of a droplet at the equilibrium state, the present numerical solution is still stable at τ = 0:5001 (very low viscosity value) by using the MRT-EDM scheme, which is more less than other two methods considered.

Model Validation.
Predicting the interfacial dynamics by considering wetting surfaces in a liquid-vapor system is essential for simulation of multiphase flow transport through porous media. In order to demonstrate the capability of the S-C LBM employed for the numerical solution of such a flow problem, the simulation results of the equilibrium state of a droplet on a solid wall and a two-phase flow through a channel are verified with considering wettability effects. According to the accuracy and performance study of the different techniques considered in the previous section, the MRT-EDM scheme is used in this validation study by applying the CS EoS and employing the modified interaction force with the β parameter.

Geofluids
The simulation of the equilibrium state of a droplet on a solid wall is performed in a square domain with the grid size 300 × 300, in which the periodic condition is considered in the left and right sides, and the wall boundary condition is imposed on the top and bottom of the flowfield. The study is carried out for a droplet placed on the bottom wall with radius 45 ðluÞ at different wetting conditions. To verify the numerical scheme employed, the contact angle θ w between the droplet interface and the solid surface is computed for different normalized density D ρ which is defined by In Figure 8, the results obtained based on the present solution algorithm are compared with those reported in Ref. [12]. As shown in this figure, the wall surface becomes hydrophilic by increasing the ρ w , which leads to the lower contact angle θ w , and vice versa. It should be noted that the linear dependency between the value of D ρ and the computed contact angle is also reported in the literature [20,46,47] which the present results show good agreement.
The second test case for model validation is the simulation of the two-phase fluid flow through a duct which can be considered as a simplified pore of a porous medium. The transmitting of a two-phase flow within a channel causes the wetting phase to flow near the hydrophilic wall which forms a film of wetting fluid on the surface. Consequently, 9 Geofluids the nonwetting phase is forced to flow through the central core of the channel. The schematic of the two-phase flow inside a 2D channel is demonstrated in Figure 9. The parameters a and b in this figure indicate the height of the nonwetting phase and the height of the channel from the centerline, respectively. Herein, the periodic boundary condition is considered in the streamwise direction, and the bounce-back boundary condition is imposed on the solid walls.
For this two-phase flow problem, the fully developed flow condition is obtained by applying an external body force F. Therefore, the analytical solution of the fully developed flow velocity for the two-phase flow in the cross-section of the channel is defined as follows [46]: where υ w , ρ w , υ nw , and ρ nw are the kinematic viscosity and density of the wetting and nonwetting phases. According to Eq. (25), the velocity profile of the fully developed twophase flow in the channel is related to the dynamic viscosity ratio M = μ nw /μ w and the saturation of each phase [46], which can be defined as In the present study, the fully developed velocity profile of the two-phase flow through the channel obtained based on the numerical solution is compared with the analytical data in Figure 10     10 Geofluids predicted by the present MRT-EDM scheme has about 9.5% and 6% deviation from the analytical profile at M = 0:1 and M = 0:32, respectively. Huang and Lu [48] have reported about 7% difference between their numerical solution and the analytical data for the velocity profile of the two-phase flow through the channel. This comparison shows that the present numerical technique is accurate enough for the simulation of the two-phase flow pattern in 2D pore scales.
To assess the capability of the implemented S-C LBM for computing the two-phase flow characteristics in the pore scale, the relative permeability of the wetting and nonwetting phases through the channel is computed at different flow conditions. Herein, the definition proposed by Yiotis et al. [46] is used for calculating the relative permeability of each phase as where k r,w and k r,nw denote relative permeability of wetting and nonwetting phases, respectively. These relations define 11 Geofluids that in a 2D channel, k r,w only depends on the saturation of the wetting phase, while k r,nw is related to the saturation of the nonwetting phase and the dynamic viscosity ratio. Figure 11 illustrates that the relative permeability obtained based on the present numerical scheme is in excellent agreement with those computed with the analytical relations at different wetting phase saturations. A similar study is performed at different relative temperatures T/T c to examine the accuracy of the S-C LBM for predicting the relative permeability magnitude in various density ratios. Figure 12 shows the comparison of the relative permeability computed by the present numerical scheme with those obtained by the analytical solution as a function of the reduced temperature.
This study shows a very good agreement between the numerical and analytical solutions with the implementation of the modified interaction forcing scheme with the parameter β = 1:18, while using the standard interaction force (β = 1) causes a serious deviation between the results at the flow conditions with T/T c < 0:7, which corresponds to the highdensity ratios. Therefore, the choice of the interaction forcing scheme becomes important for such flow conditions. This validation study shows that the S-C LBM with employing the MRT-EDM scheme, CS EoS, and the modified forcing scheme with parameter β is efficient and accurate for estimation of the two-phase flow structures, and characteristics  12 Geofluids through the channel and can be used for simulation of such flow problems in porous media.

Two-Phase Flow Transport through a Porous Medium.
The present numerical scheme developed based on the S-C LBM with employing the CS EoS, MRT-EDM, and modified forcing techniques is applied for studying the two-phase flow transport through a porous medium. The capillary number C a and the dynamic viscosity ratio M are considered as two effective parameters that characterize the two-phase flow structure in this study [49]. Herein, the capillary number is defined by where μ nw and u nw are the dynamic viscosity and velocity of the advancing nonwetting fluid, respectively, and θ is the fluid-fluid contact angle. Herein, a homogeneous porous medium is generated by a distribution of obstacles in a flow domain with 1000 × 300 grid points. The pore space area and the height of the throat between the obstacles include 20 × 20 lattice and 4 ∼ 7 grid points, respectively. These considered grid points are required for the spatial discretization of a pore to perform the solution accurately enough for capturing the interfacial dynamics in the pore scale. Hence, the existence of very small pores in the geometry of a porous medium needs a high-resolution grid size and increases the computational cost. The geometry of the porous medium and the boundary conditions of the flow domain are shown in Figure 13. The velocity inlet and pressure outlet boundary conditions are imposed in the left and right boundaries, respectively, and the top and bottom sides are considered periodic. The wetting fluid is defined by the higher density (red color) which invades the nonwetting fluid with lower density (blue color) in the flow domain. Such an invasion produces some fingering regimes that their origin can be the surface or geometrical characteristics of a porous medium in nature, e.g., the nonuniform surface wettability of a sandstone [50]. Herein, the fingering regime during the two-phase flow simulations through the considered porous medium is occurred due to the random generation of the shape of obstacles and the wetting condition of the solid surfaces. Figure 14 demonstrates the two-phase flow patterns obtained by the present numerical scheme based on the S-C LBM at the dynamic viscosity ratio M = 1017 with different capillary numbers. To define the capillary number, the surface tension is calculated by using the Laplace law analysis, which is obtained σ = 0:02316. Also, the fluid-fluid contact angle is measured about 175 ∘ which shows fully hydrophobicity of the invading fluid. This study is performed at the reduced temperature T/T c = 0:482, and the parameter β is set to be 1.241 to preserve the thermodynamic consistency of the solutions. From the results presented in Figure 14, it can be seen that the penetration of the invading fluid at a certain time t = 50000 is decreased by reducing the capillary number. These patterns are discussed with more details along with the explanations of Figures 15 and 16 at the end of this section.
To compare the obtained results for the two-phase flow characteristics through the porous medium in a similar situation, a normalized time T * is defined by division of the time step to the total simulation time of each case. The total 13 Geofluids simulation time is counted until the first finger of the invading fluid touches the outer side of the porous medium for each case, which is called breakthrough [49]. Figure 17 illustrates the variation of the pressure difference Δp between the inlet and outlet of the porous medium versus the normalized time T * at different capillary numbers. As seen in this figure, by increment of the capillary number in a constant dynamic viscosity ratio, the invading fluid has enough pressure difference to overcome the throat threshold pressure p c presented by Lenormand et al. [49], which is known as a resistance factor. The throat threshold pressure is defined by the following relation: where σ is the surface tension, θ is the fluid-fluid contact angle, and r is the throat diameter. This parameter also affects the nonwetting phase saturation. In the present study, by considering the minimum and maximum throat radius, the threshold pressure is altered between -0:006591 and -0:011535 for r = 7 and r = 4, respectively. The negative sign indicates the resistance pressure against the invading fluid. Figure 18 shows the nonwetting phase saturation as a function of the normalized time T * . The results obtained demonstrate that the saturation of the nonwetting fluid is enhanced by increasing the capillary number, which is consistent with the previous works [49][50][51]. The final values of the nonwetting phase saturation and the pressure difference at the breakthrough condition are illustrated in Figure 19 at    different capillary numbers. It is obvious that the variation of the nonwetting phase saturation is almost linear which is in agreement with the results obtained by Zhang et al. [51] and Liu et al. [50]. The agreement of the present results in comparison with those reported in the literature confirms that the implemented S-C LBM with the MRT-EDM scheme besides the modified interaction forcing method and the CS EoS is capable for simulation of two-phase flows at highdensity and -viscosity ratio through porous media at different flow conditions.
The structure of the two-phase flow through the porous medium is also investigated at different flow conditions. The displacement pattern of the invading fluid is categorized as the capillary fingering, viscous fingering, and stable displacement, depending on the capillary number and the viscosity ratio between two phases [52]. A phase diagram is used to describe these flow structures by log ðCaÞ and log ðMÞ for a gas-water flow system (see Figure 15). In the present simulations based on the S-C LBM, the magnitude of the viscosity ratio is set as log ðMÞ > 0. Therefore, the capillary fingering and stable displacement patterns are expected in the porous medium considered. As can be seen in Figure 14, the two-phase flow pattern is altered from the stable displacement to the capillary fingering by a decrement of the capillary number. Figure 16 demonstrates the present results for instantaneous two-phase flow displacement at different times for log ðCaÞ = −1:85, in which the pattern behavior is almost stable displacement. In this flow condition, the interface of the invading fluid penetrates into the porous medium along a straight line, despite some   irregularities due to the pore-scale geometry. This behavior can be the result of the dominant effect of the viscous force and the higher pressure difference that cause to fill the porous medium uniformly and simultaneously. By decreasing the capillary number, it is observed that the invading fluid is extended in the shape of irregular fingers in the flow field due to the dominant effects of the capillary force in comparison with the viscous force. Such a capillary fingering of twophase flow through the porous medium is shown in Figure 20 at log ðCaÞ = −3:5. In this case, due to the uniform pressure distribution throughout the domain, the nonwetting fluid only penetrates the pore-scale area with a lower throat pressure threshold. As a result, the direction of the flow displacement in the porous medium is irregular. In Figure 20, it can be seen that the capillary fingers of the nonwetting fluid form several loops that trap some slugs of the wetting fluid among the porous medium. These two-phase flow structures obtained using the present numerical scheme are in consistent with those reported in the literature [49,50]. The interfacial area between two phases is quantified for the two-phase flow through the porous medium at different capillary numbers. This parameter is an important factor that influences the mass and energy transfer among the phases [50,51]. In the present 2D studies, the interfacial length is measured instead of the interfacial area. Herein, due to the extreme hydrophilicity of the wetting fluid, the length of the wetting films around the obstacles is also added to the length of the measured interface between the phases for computing the interfacial length. To have comparable results for this parameter, a specific interfacial length is calculated as the ratio of the measured interfacial length to the pore area of the porous medium. Figure 21 shows the variation of the specific interfacial length as a function of the nonwetting phase saturation at different capillary numbers. The present results demonstrate that there is a linear relation between two parameters in all the capillary numbers considered which is in good agreement with the previous experimental finding of Zhang et al. [51].

Conclusion
In this work, several different combinations of collision operators, EoSs, and forcing schemes of the pseudopotential multiphase LBM are employed to investigate their performance and capability for computing some two-phase flow characteristics at high-density and -viscosity ratios with the aim of applying for porous media applications. Herein, the SRT and MRT collision models, Shan-Chen (SC) and Carnahan-Starling (CS) EoSs, and SVM and EDM forcing schemes are considered to study their accuracy and numerical stability. An appropriate modification in the forcing schemes is also used to remove the thermodynamic inconsistency in the simulation of two-phase flow problems studied at low reduced temperatures. The simulation results obtained for the equilibrium state of a droplet, the two-phase flow characteristics in a channel, and through a porous medium are presented and discussed. The conclusions and remarks regarding the present study are as follows: (1) The Laplace law is satisfied for simulation of the stationary droplet by using the S-C LBM coupled with the SRT-SVM, SRT-EDM, and MRT-EDM schemes which demonstrates the accuracy of the numerical schemes employed for predicting the characteristics of a stationary two-phase flow system (2) The results obtained for the stationary droplet at lowdensity ratios (high reduced temperatures) by the SRT-SVM, SRT-EDM, and MRT-EDM schemes with employing both the SC and CS EoSs are in good agreement with the analytical solution on the coexisting curve. However, at high-density ratios (low reduced temperatures), the results have a significant deviation from the analytical solutions with the implementation of the forcing schemes in the standard form. By employing the β-modification in the forcing schemes, the density values obtained by the present numerical simulations are in good agreement with those obtained by the analytical solution in a wide range of reduced temperatures (3) The results presented in this study show that with applying the CS EoS for simulation of considered two-phase flow problems at high-density ratios, the spurious velocity is significantly decreased in the flowfield. This result demonstrates that the CS EoS is more realistic and efficient than SC EoS for simulation of two-phase flows at high-density ratios due to producing less spurious currents at the interfacial region (4) It is shown that the MRT-EDM with the improved forcing scheme produces minimum spurious currents among the other methods. By considering the results presented in this study, the advantages of employing the MRT scheme with β − modification of forcing schemes are obvious for simulation of a two-phase flow system at low relaxation times (τ ⟶ 0:5) (5) According to the previous remarks, the present study demonstrates that the S-C LBM with employing the MRT-EDM scheme, CS EoS, and the modified forcing scheme is efficient and accurate for estimation of the two-phase flow characteristics. Therefore, this combination of the implemented numerical techniques is recommended and applied for studying two-phase flows through the channel and porous medium (6) The validation study of two-phase flow through the channel shows that the S-C LBM with employing the MRT-EDM scheme, CS EoS, and the modified forcing scheme with parameter β is efficient and accurate for estimation of the two-phase flow structures and characteristics in pore scales (7) The agreement of the present results obtained for the simulation of two-phase flow through the porous medium at different conditions in comparison with those reported in the literature confirms that the suggested S-C LBM is capable for predicting the 16 Geofluids displacement patterns, the nonwetting phase saturation, and the penetration structure of the invading fluid in porous media applications. The extension of the present method in the three-dimensional framework is straightforward and will be presented in future works

Data Availability
The data of this work are available by the request of researchers via email to the corresponding author.

Conflicts of Interest
The authors declare that they have no conflicts of interest.