Equation of State’s Crossover Enhancement of Pseudopotential Lattice Boltzmann Modeling of CO 2 Flow in Homogeneous Porous Media

: The original Shan-Chen’s pseudopotential Lattice Boltzmann Model (LBM) has continu-ously evolved during the past two decades. However, despite its capability to simulate multiphase ﬂows, the model still faces challenges when applied to multicomponent-multiphase ﬂows in complex geometries with a moderately high-density ratio. Furthermore, classical cubic equations of state usually incorporated into the model cannot accurately predict ﬂuid thermodynamics in the near-critical region. This paper addresses these issues by incorporating a crossover Peng–Robinson equation of state into LBM and further improving the model to consider the density and the critical temperature differences between the CO 2 and water during the injection of the CO 2 in a water-saturated 2D homogeneous porous medium. The numerical model is ﬁrst validated by analyzing the supercritical CO 2 penetration into a single narrow channel initially ﬁlled with H 2 O, depicting the fundamental role of the driving pressure gradient to overcome the capillary resistance in near one and higher density ratios. Signiﬁcant differences are observed by extending the model to the injection of CO 2 into a 2D homogeneous porous medium when using a ﬂat versus a curved inlet velocity proﬁle.


Introduction
Carbon Dioxide (CO 2 ) and other greenhouse gases (i.e., methane, nitrous oxide, among others) are key contributors to climate change. According to the Energy Information Administration (EIA) estimates, the global emission of CO 2 alone can rise to 6.41 billion tonnes by 2030 [1]. Therefore, various mitigation measures have been proposed to tackle global climate change. Among those proposed solutions, we use CO 2 in different applications such as Hot-Dry Rock (HDR) systems, heat engines, selective extraction processes, reclamation processes of metals from industrial effluents, and chemical reactions [2,3]. Nevertheless, CO 2 sequestration is recognized as one of the key strategies for reducing CO 2 emission rates. Due to its low critical temperature (31 • C), CO 2 is in supercritical (for example, in CO 2 sequestration) or in near-critical (as in air conditioning and cooling systems [4,5]) state.
The wide usage of CO 2 around the critical state requires a deeper investigation of the fluid thermodynamics under and above critical conditions, affecting fluid transport properties [1]. In addition, it facilitates understanding the industrial and laboratory processes wherever the fluid is involved since the behavior of fluids in supercritical or near-critical conditions dramatically differs from those with constant properties.
The fluid flow under different conditions can be studied experimentally and numerically. However, as the experimental study may not always be feasible and is timeconsuming, the numerical analysis can be used as a starting point. It can even replace the experiments to provide valuable insights into the complex flow patterns.
Computational Fluid Dynamics (CFD) techniques such as the finite volume, volume of fluid, and level-set methods solve the discrete Navier-Stokes equations. However, these methods become computationally expensive when applied to the interface tracking of transient multiphase flows [6,7]. The mesoscale Lattice Boltzmann Model (LBM), specifically the LBM employed in the current work, does not explicitly track the phase interface; instead, it includes inter-particle forces naturally preserving the interfaces [8]. Additionally, it accommodates complex boundaries such as porous media and narrow channels. Furthermore, the pseudopotential LBM can incorporate real Equations of State (EoS) such as Van der Waals (VdW), Carnahan-Starling (CS), and Peng-Robinson (PR) EoS [9], which, in turn, allows the model to match the fluid characteristics more closely.
Despite its advantages, the pseudopotential LBM still faces issues when high density and viscosity ratios are included in multicomponent-multiphase flows in complex geometries. For example, studies on displacement patterns in porous media do not consider density contrast and differences in the fluids' critical properties [10][11][12][13][14]. Hence, most numerical studies do not consider that CO 2 and H 2 O are in different thermodynamic states. For example, CO 2 is usually in a supercritical state during sequestration, while H 2 O is in a liquid state, limiting the accuracy of the numerical results.
In addition, most studies apply classical EoS to analyze the thermodynamic properties of the fluid in different ranges of temperature [9,[15][16][17]. Though classical cubic EoS predict well the fluid properties away from the critical region, they fail to accurately capture the fluid behavior in the near-critical region or at a supercritical state because of long-scale fluctuations in density that cover a wide range of temperatures [1,[18][19][20][21][22]. In this region, analytical solutions provided by classical EoS are no longer accurate. As a remedy for this pitfall, a so-called "crossover" formulation of the EoS was proposed, incorporating the non-analytical scaling laws and predicting fluid properties in the critical region more accurately [1,18,21]. Any classical EoS can be re-formulated in a crossover fashion. In this work, the crossover formulation is incorporated into PR EoS. As a result, PR EoS better characterizes the specific fluid through an acentric factor. This approach's advantage is a hybrid equation that accounts for both non-analytical and analytical fluid behavior, close and far from the critical point, respectively. It can also transform back to ideal gas EoS at a zero-density limit.
Our previous study [23] incorporated the crossover PR EoS into the standard pseudopotential multicomponent LBM. We showed that crossover PR EoS improved the model accuracy for fluids at near/supercritical conditions compared to the classical PR EoS. Another previous work [24] showed that a higher water wettability promotes a higher effective CO 2 penetration through a porous medium. We then applied the improved model formulation to study the transport of a supercritical CO 2 bubble within a simplified hydrophilic porous medium saturated with H 2 O, proving the capacity of the newer model crossover formulation.
This work extends our work by introducing the PR EoS multicomponent LBM to study the penetration mechanisms of a continuous front of supercritical CO 2 through a pore-scale homogeneous domain saturated with H 2 O. Firstly, a 2D narrow channel is considered to model the penetration of a supercritical CO 2 bubble through the channel with a significant density ratio. Secondly, the sequestration of a continuous flow of CO 2 is modeled by injecting it through a water-saturated 2D homogeneous hydrophilic porous medium (contact angle of 70 • ), varying the pressure gradient around the pressure drop needed to overcome the capillary resistance. Finally, we study the CO 2 injection under a fixed flow regime and a mildly hydrophobic 2D porous medium (contact angle 85 • ) to demonstrate the effect of the injection profile in the CO 2 penetration into the watersaturated porous medium. The porous media was constructed by placing circle-obstacles homogeneously. Hence, this study represents the first investigation, to the very best of our knowledge, featuring a crossover PR EoS formulation of LBM to shed some light on the CO 2 -H 2 O-surface interaction and CO 2 sequestration in a porous domain at near supercritical conditions. This paper is organized as follows: Section 2 introduces the multicomponent LBM and crossover PR EoS formulation. Section 3 presents the model's performance for the static H 2 O droplet immersed in CO 2 and dynamic conditions of CO 2 injection at supercritical conditions into H 2 O-saturated capillary domains. The final discussion and conclusions are presented in the last section of the paper.

General Lattice Boltzmann Model
The LBM equation represents the evolution in time of distribution functions because of fluid particles streaming and colliding.
The most general formulation of LBM for a multicomponent system is presented in Equation (1): where f i j (x, t) is the density distribution function at point x and time t with the discrete velocity in i direction for component j. In Equation (1), the left-hand side represents the streaming step, while the right-hand side is the collision step, which can assume different forms: in this work, we opted for the simplest one, the Bhatnagar-Gross-Krook (BGK) collision operator [25], which each fluid species is characterized by a single relaxation time τ j which is related to the whole mixture kinematic viscosity through the relation (Equation (2)) [26]: where c s is the lattice sound speed and χ is the local value of the mass fraction. The collision operator relaxes the local distribution function towards its equilibrium value, f j,eq i (x, t) which depends on macroscopic quantities, as prescribed by Kinetic Theory, and is given by Equation (3): In this work, we used a two-dimensional formulation with nine local speeds at each computational node. This set of speeds, referred to as a lattice, must be carefully chosen to recover the Navier-Stokes equation. If the lattice vectors possess different modules, weighting factors w i need to be defined to achieve this goal. The number of dimensions and speeds defines a generic lattice. We chose the standard D2Q9 one, depicted in Figure 1, characterized by weighting factors 4 9 (i = 0), 1 9 (i = 1, 2, 3, 4), 1 36 (i = 5, 6, 7, 8). Macroscopic quantities, including density ρ j , can also be obtained as discrete moments of the particle distribution f i : where u is the macroscopic velocity vector. The difference between continuous distribution function f and f i is the discrete nature of the last one; subscript i refers to the discrete set of velocities e i , and superscript j refers to a component.

Pseudopotential Lattice Boltzmann Model
In Shan-Chen's or pseudopotential model, the interaction between the particles is modeled using pairwise interaction potentials. These interaction potentials are functions of density for each component. If we define the interaction matrix G jk between two generic species j and k at a generic node x, ψ as a pseudopotential ("effective mass"), then the force acting on jth component can be written as [27]: The interaction force is given by Equation (6), and it is applied between the nearest neighbor particles and leads to the pressure P j in the EoS: where the first term represents the ideal EoS and the second term accounts for repulsive (positive g jk ) and attractive forces (negative g jk ) between the fluid particles. For a two-component system, usually the diagonal terms of interaction matrix, accounting interactions between particles of same species are equal to zero, while interparticle terms are chosen positive to enforce repulsion, as shown in Equation (8).
The pseudopotential ψ in Equation (7) can be tailored to obtain a given equation of state. Pseudopotential can also be obtained by rearranging Equation (7): Yuan and Schaefer incorporated different EoS such as VdW, PR, and CS into the pseudopotential model through Equation (9). In this case, the term g is used to ensure the positive sign under the square root [9].
The fluid spreading on a solid surface depends on the strength of fluid-solid interaction. The surface wettability is modeled using Equation (10) proposed by Martys and Chen [27] in multicomponent form: The g wall represents the fluid-wall interaction strength. For example, when g wall is positive, the wettability decreases, while if the g wall is negative, the wettability increases. s(x + e i ) is the "switch" function equal to 1 in the presence of the solid surface at a lattice node x; otherwise, it is 0. The following method is based on density instead of potential ψ(x) [28].
The total force acting on a fluid particle is described with the equation: The F b is the external body force, which can have a constant value if applied to drive the flow, for example, or zero, if no external body force is present. The present work adopts the velocity-shift force scheme, originally used by Shan and Chen [8]. In this scheme, forces are incorporated by updating the velocity u in the equilibrium distribution function, replacing it with u eq , giving: where the actual physical velocity u is taken from Equation (5). Other available schemes are the Exact Difference Model (EDM) [17,29] and Guo's [30]. These schemes result in a further term at the right-hand side of (1), the source term S i (x, t).

Crossover Peng-Robinson EoS
According to Landau's theory, free energy is an analytic function of the order parameter. However, phase transition disregards the order of the parameters, and the partition function becomes the sum of all possible values. Consequently, an analytical solution in the region near the critical point cannot be formulated [18,31]. On the other hand, the non-dimensional expression of Helmholtz's free energy is expressed in Equation (13) in the crossover formulation of EoS. Moreover, it is made up of two terms [18], where the first term represents the non-analytical part ∆A, and the second term considers the analytical part µ 0 (T) that depends on temperature. where is the non-dimensional pressure at the critical isochore (V = V c ), and both ∆T = T r − 1 and ∆V = V/V c − 1 are dimensionless deviations from classical critical parameters. The critical parameters are found by solving Equation (14): where Z c is the critical compressibility factor. After imposing the condition expressed in Equation (15), the critical part of the Helmholtz free energy takes the form shown in Equation (16) with the introduction of dimensionless coefficients b and critical exponent α: The analytical part is given as: The term A 0 (T) in Equation (17) is expressed as A 0 (T) = A 0 (T)/RT and defined as a non-dimensional function of temperature that plays a role in reproducing ideal gas behavior in the limit of zero density. Helmholtz free energy in crossover form uses the renormalized form of classical reduced parameters ∆T and ∆V as shown in Equations (18) and (19): where τ = T/T crit − 1 is the temperature difference from the real critical temperature (T crit ) and ∆η = V/V crit − 1 is the deviation of the volume from the real critical volume (V crit ). Both Equations (19) and (20) use the real critical temperature and volume, and their second term accounts for the difference between classical and real critical parameters. Moreover, these shift parameters are expressed as ∆τ crit = T crit /T c − 1 and ∆η crit = V crit /V c − 1. As in [22], we assume that T c = T crit and therefore, ∆τ crit = 0. Another important parameter in Equations (17) and (18) is a crossover function, Y, that accounts for volumetric changes in the EoS depending on how far the system is from the critical region. If Y(q) → 0 , when the system is at the critical point, the crossover function should approach zero, Y(q) → 0 , and Y(q) → 1 when the system is away from the critical point. The crossover function proposed by Feyzi et al. [21] satisfies these requirements, and we use the same formulation in Equation (20).
where q is the renormalized parameter that measures the distance between the system volume and the critical one; q is found solving the Equation (21), as in [20].
and Gi is the Ginzburg number, while d 1 , v 1 , and m 0 are systemdependent parameters. Equation (21) has been solved numerically using MATLAB, and as it has several solutions, we take the highest positive value among all possible solutions. After finding q and Y(q), the renormalized parameters τ and ∆η are substituted into Equation (16) to obtain the critical form of the Helmholtz free energy.
Finally, the crossover PR EoS is obtained after substituting all the parameters and taking the differentiation with respect to the volume, Equation (22): As mentioned earlier, the critical shift ∆τ crit is set to 0 and ∆η crit , Gi, d 1 , v 1 , m 0 , a 20 , and a 21 are parameters that depend on the type of fluid, and their values are shown in Table 1 for CO 2 and H 2 O. The crossover PR EoS is incorporated into the LBM through Equation (22). For CO 2 , all parameters are taken from [21], while all other fluid parameters are taken from [20]. Furthermore, the physical critical parameters are derived from the NIST database [32].

Results
Our previous study [23] conducted several 2D-droplet tests at different reduced temperatures, and their respective coexistence curves have been constructed. The detailed comparison of the thermodynamic behavior of CO 2 and H 2 O at various temperatures predicted by PR and crossover PR EoS showed that PR EoS dramatically improves the model accuracy, which is more noticeable in the liquid phase. PR generally overestimates the density of the curve's liquid branch, although it predicts the gas phase accurately. As temperature reduces, the difference in liquid density predicted by the PR EoS increases, following the experimental data's behavior [23]. Similar results are also noted in other studies [33,34]. It is also observed that the crossover PR EoS matches experimental data much better than the classical PR EoS at different temperatures close to the critical point. The results show that crossover PR EoS is more accurate than classical PR EoS because it considers the non-linear behavior of fluid species in the critical regions [23].

Multicomponent System
Modeling the multicomponent system with actual fluid parameters is challenging mainly due to a lack of thermodynamic consistency and potential large spurious currents at interfaces. Moreover, the magnitude of these spurious currents increases with the density ratio. Therefore, studies that involve complex systems such as flow in porous media are mostly limited to fluids with similar densities [10,35,36].
In our previous study [23] and the current work, we implement the improvements to the MCMP model following Stiles and Xu's [15] approach, which improves the model in terms of numerical stability and allows a density ratio higher than unity. The improved model is then applied to simulate the stationary 2D droplet immersed in supercritical CO 2, accounting for density contrast and the fluid's critical temperature difference. We use a so-called scalar multiplier as suggested by Stiles and Xue [15] with the assumption that ∆t = ∆x = 1. The pressure term in the numerator of Equation (9) is replaced with P * i = kP i . We ensure stability of the model by using relaxation time of 1, which is a requirement of the formulation.
In addition, using crossover PR EoS helps better match the binodal curve. It considers fluid properties by applying the acentric factor and system-dependent parameters to improve the critical region's curve adjustment. The acentric factor for H 2 O is 0.344, and the set of parameters used for CO 2 is given in Table 1.
Stiles and Xu investigated the effect of scaling parameters on the interaction force model and the coexistence curve (P versus V). As expected, the curve's amplitude decreases when k < 1 and the volume of the liquid phase remains constant. This parameter only affects the volume of the gas phase. This effect becomes more intense at lower temperatures, thus limiting the range of application to moderate and high temperatures [15,23]; in this study, we assume that k = 0.75.
To account for differences in real fluids, critical temperatures T c, H 2 O T c, CO 2 ≈ 2 , we use the set of parameters presented in Table 2. Previous validation [23] shows supercritical CO 2 co-existing with subcritical H 2 O at 0.071 (lattice units-LU-temperature), corresponding to T r of 2.07 and 0.97, respectively. Hence, we first set the critical lattice temperature of H 2 O to 0.07292 and then use the relation in Equation (20) to find the critical temperature of CO 2 .  Figure 2 presents a 2D static water droplet surrounded by CO 2 in the center of a rectangular domain with the size of 201 × 201 LU 2 , where the system's temperature is 0.071 in LU. Here, g ij = g ji = 0.001, as it gives good stability and is recommended in [15], and the drop radius is 30 LU.  Figure 2 shows a slight trace of H 2 O in the CO 2 and vice versa, but the CO 2 , being the primary phase, mainly dictates the interfacial dynamics. The presence of trace elements in small amounts helps to reduce the initial gradient at the phase interface and improves numerical stability.

Contact Angle Test
This section presents the contact angle tests and shows that the improved model is valid for simulating fluid flows involving fluid-solid interactions. Figure 3 shows that the multicomponent multiphase model can reproduce a wide range of contact angles. The properties of both fluids are the same as described in Section 3.1.
The interaction strength parameter between the solid and fluid g CO 2 ,wall is 0.2, and for the g H 2 O,wall is −0.2, accordingly. In this case, the H 2 O droplet behaves as a wettable fluid, and the CO 2 is non-wettable. The values of the g CO 2 ,wall and g H 2 O,wall used in this study to perform the contact angle tests are presented in Table 3. The results show that the pseudopotential model can model a wide range of contact angles by adjusting the solid-fluid interaction parameters, g CO 2 ,wall and g H 2 O,wall .

Penetration Process in 2D Narrow Channel
In this section, we analyze the performance of the multicomponent LBM model using the crossover PR EoS formulation in the visualization of the penetration of CO 2 into an H 2 O-saturated 2D funneled channel that helps approximate the contractions in a pore network. CO 2 is a non-wetting phase, and H 2 O is set as a wetting phase. The geometry of the domain is depicted in Figure 4a-c filled with immiscible H 2 O and CO 2 fluids subject to a streamwise pressure gradient (body force) and inlet/outlet periodic boundary condition. The domain size is 41 × 151 LU 2, and the channel is narrowed by placing two-block obstacles with a height of 40 LU and a width of 10 LU. Figure 4a also presents the initial condition with two large CO 2 bubbles immersed in H 2 O and placed initially at the entrance and exit of the channel. H 2 O and CO 2 fractions are 0.61 and 0.39, respectively, and the system's temperature is 0.071 in LU.
The density ratio H 2 O/CO 2 is 4, and H 2 O is set as the wetting phase with a contact angle of 70 • . The top and bottom sides of the domain are set as periodic and left, and the right walls are set with a no-slip boundary condition. The streamwise body force is then applied throughout the domain with magnitude G f = 4.8 × 10 −5 just enough to overcome the capillary resistance. Figure 4b-d show the evolution of the flow with the bubble squeezing in a periodic pseudo-steady regime. The results show that the velocity of the CO 2 bubble is almost twice faster in a narrower section of the channel, as observed in Figure 5. From further analysis, it is observed that bubbly flow tends to form when CO 2 saturation is lower or equal to 0.4. At CO 2 saturation between 0.4 and 0.7, the annular flow is formed where wetting phase water stays closer to the channel wall. These findings are similar to what was observed in another similar study where free energy LBM was used instead [37,38].  The model is further challenged by varying the pressure gradient from insufficient to a much larger value than needed to overcome the capillary resistance. Figure 6 presents the snapshots after a long transient simulation (297,500 ts) to reach a pseudo-steady solution. The boundary conditions and wettability are like those in the previous case. Figure 6a shows that under a pressure gradient of 3.5 × 10 −5 , the CO 2 cannot enter into the narrow channel due to insufficient driving force to overcome the capillary resistance. The body force is further increased to 4 × 10 −5 (Figure 6b), slightly above the needed capillary pressure. The CO 2 squeezes through the narrow channel until an annular flow is developed with a stable symmetrical shape. The centerline velocity in the middle (narrow) domain section is 0.016 LU/ts, twice higher than in the broader section of the domain. When the body force increases to 6 × 10 −5 , the multicomponent flow becomes unstable, and the inner CO 2 starts oscillating in the spanwise direction, depicting a wavy shape surrounded by a film of H 2 O (Figure 6c).

Penetration Process in the 2D Pore Network
The modeling CO 2 penetration into a 2D porous media under different body forces was also studied using the crossover PR EoS formulation embedded in the multicomponent LBM. The porous media is constructed by placing several previously saturated cylinders with H 2 O. The schematic diagram of the porous region is shown in Figure 7a, which has a porosity of 0.78 and a pore throat of 10 LU. The domain size is 191 × 81 LU 2, with the horizontal length of the pore network being 90 LU. Top and bottom walls are set as a no-slip boundary condition, while both left and right sides are set as periodic. The computational domain is subject to a streamwise pressure gradient to drive the CO 2 into the porous medium saturated with H 2 O (wettable fluid) with a contact angle of 70 • . The density ratio is set to 1, and the system's temperature is 0.071 in LU. Figure 7b presents the case with body force 2.5 × 10 −5 at 113,400 ts. In this case, the driving force is not enough to overcome the capillary resistance. We further increase the driving force to 4.5 × 10 −5 , and it is observed from Figure 8 that CO 2 can then easily intrude the porous media with a flow pattern that resembles a stable displacement process. Next, we compare the effect of the front profile of the intruding fluid (CO 2 ). According to Fakhari et al. [39], the distance from the inlet to the porous section of the micromodel has a significant effect on numerical simulations. The parabolic velocity profile better matches the initial conditions between the numerical simulations and experiments. To achieve a parabolic profile, we extended the length of the inlet region to 85 LU. To observe the effect of the front shape of intruding fluid on displacement patterns, we also placed the inlet region closer to the pore network at 7 LU from the first row. The initial conditions for both cases are shown in Figures 9a and 10a. In both cases, the overall domain size is 401 × 201 LU 2 , with the pore network length of 229 LU in the streamwise direction. The inlet velocity is 0.02 LU/ts as proposed by Zou and He [40], applied on the left bound of the domain (inlet), and the outflow boundary condition is applied on the right end. Both top and bottom walls are set as a no-slip boundary condition, and the contact angle is set to 85 • . The density ratio is again set to 1, and the system's temperature is 0.071 in LU.    Figure 9d shows that the flow pattern resembles a fingering displacement at 9000 ts. The process occurs at near-constant pressure conditions around the evaporation region. Thus, the model is expected to be capable of capturing the phenomenon despite the inherent isothermal assumption. Figure 10a presents the case with the flat front of the intruding fluid placed 7 LU away from the porous section. Compared with the previous case, the front profile of the velocity remains flat. In contrast, with the parabolical front profile, the velocity is higher in the middle section of the porous media. These results suggest that the distance from the inlet to the porous section can affect the displacement pattern.
A grid-independence analysis is developed by introducing two more refined models of 601 × 301 LU 2 and 801 × 401 LU 2 , with geometry, inlet speed, viscosity, and timestep span scaled accordingly. The average velocity magnitude along a line drawn over the last row of obstacles, as shown in Figure 11, was utilized as a comparison parameter for the three grid sizes. Results of the velocity probe are shown in Figure 12, taken from lattice time of 0 to 7100.
As shown in Figure 12, velocity magnitude is very similar in all grid sizes, deviating only slightly in later time-steps, demonstrating the suitability of any of the three grids in our analysis.
In another comparison between the grids, we used the CO 2 flux through the probe line for an extended period to permit the CO2 to reach the probe line. This parameter was calculated as the fraction of CO 2 multiplied by the time-average velocity magnitude integrated over the whole probe line. Results, shown in Table 4, highlight the small error between contiguous grid resolutions, justifying our initially selected grid for the analysis.   Figure 11) in the CO 2 penetration LBM model. Table 4. Time-space average CO 2 flux over a probe line ( Figure 11) integrated over 7100 timesteps.

Discussion and Concluding Remarks
In this study, the crossover PR EoS was incorporated into the pseudopotential LBM, and an improved model was further applied to study the CO 2 penetration process in capillary water-saturated domains. The improved model's performance is analyzed through several tests. Firstly, a static weightless water droplet is suspended in equilibrium and immersed in CO 2 to prove the stability of the numerical model. Secondly, a sessile water droplet in equilibrium with the surrounding CO 2 is presented to demonstrate the stability and accuracy of the model in determining the static contact angle. Both static simulations showed the model's ability to describe different wettability regimes with sufficient numerical stability.
Then, several cases are presented to analyze the injection of supercritical CO 2 in a water-saturated capillary domain. In the first place, a validation of the model around critical conditions is assessed through a single pore or capillary channel, demonstrating the accurate prediction of the threshold pressure gradient needed to overcome constraints to allow a CO 2 bubble to squeeze through the narrow channel. Hence, the injection process was analyzed at different pressure gradients and contact angles. After the necessary pressure gradient was met, the non-wettable CO 2 traveled through the center of the channel surrounded by the wettable fluid (H 2 O), depicting a centerline velocity at the narrow section almost twice higher than through the main channel. These results agree with previous work [37]. Increasing the pressure gradient much farther beyond the minimum capillary pressure affected the hydrodynamic stability of the CO 2 stream, demonstrating the existence of a secondary stability threshold of pressure gradient above which a transient spanwise oscillation sets in.
Next, the CO 2 sequestration process was further analyzed around supercritical conditions where the non-wetting fluid is injected into a complex, but homogeneous 2D pore network saturated with wetting H 2 O. Thus, we investigated the effect of a supercritical CO 2 front profile on its displacing through a 2D water-saturated pore network. The results show that the intruding fluid's front profile plays an important role and influences the overall CO 2 -H 2 O flow pattern through the pore networks. That observation was also highlighted in [41], where it was observed that the parabolical profile is more likely to happen in actual CO 2 sequestration, and this velocity profile returns results that are closer to experimental data than obtained for the flat inlet profile. Therefore, our results proved that with enough development length allocated in the inlet region of the computational domain, it is reasonable to expect a more realistic parabolical profile at the entrance of porous networks during the CO 2 sequestration numerical modeling.
These preliminary results show that the improved model is numerically stable when applied to reproduce moderately complex geometries and physical flow patterns, as expected in realistic conditions. However, our study is limited to qualitative results and a viscosity ratio of one and isothermal conditions. In future studies, we plan to continue exploring the numerical stability of the new formulation for more complex 3D porous media and even more realistic fluid-fluid conditions.
For a detailed expansion of the crossover PR EoS, please refer to Appendix 4 of [18]. Informed Consent Statement: Not applicable.

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.