of immiscible two-phase displacement in two-dimensional

: Understanding the dynamic displacement of immiscible ﬂuids in porous media is important for carbon dioxide injection and storage, enhanced oil recovery, and non-aqueous phase liquid contamination of groundwater. However, the process is not well understood at the pore scale. This work therefore focuses on the effects of interfacial tension, wettability, and the viscosity ratio on displacement of one ﬂuid by another immiscible ﬂuid in a two-dimensional (2D) Berea sandstone using the colour gradient lattice Boltzmann model with a modiﬁed implementation of the wetting boundary condition. Through invasion of the wetting phase into the porous matrix, it is observed that the viscosity ratio plays an important role in the non-wetting phase recovery. At the viscosity ratio ( λ ) of unity, the saturation of the wetting ﬂuid is highest, and it linearly increases with time. The displacing ﬂuid saturation reduces drastically when λ increases to 20; however, when λ is beyond 20, the reduction becomes less signiﬁcant for both imbibition and drainage. The front of the bottom ﬁngers is ﬁnally halted at a position near the inlet as the viscosity ratio increases to 10. Increasing the interfacial tension generally results in higher saturation of the wetting ﬂuid. Finally, the contact angle is found to have a limited effect on the efﬁciency of displacement in the 2D Berea sandstone. porous media, the pore structures are extremely complex and irregular, which allow us to further assess the capability of numerical model. We apply the colour gradient LBM for dynamic displacement of immiscible ﬂuids in a complex porous ﬂow structure extracted from a realistic Berea sandstone sample, and study the effects of the viscosity ratio, interfacial tension, and contact angle on displacement process.


Introduction
The simultaneous flow of two immiscible fluids in porous media is often found in nature and industrial processes such as non-aqueous phase liquid (NAPL) contamination of groundwater, carbon dioxide injection and storage, and enhanced oil recovery (EOR) operations. For the spontaneous imbibition of oil, carbon dioxide or water, which is immiscible to the oil phase, is usually imbibed into permeable rock formations which are characterized by low porosity and low in situ permeability. The displacement process is affected by matrix properties such as the heterogeneity, wettability, and fluid properties of each phase [1].
Experimental efforts have been made to understand interactions between fluids and rock during fluid displacement [2,3]. Theoretical works usually use simplified pore geometries which allow theoretical solutions to the microscopic flow patterns [4]. As a complementary method to the experimental and theoretical approaches, accurate and reliable numerical tools are needed to understand the underlying mechanisms of multiphase displacement. Conventional macro-scale simulations have achieved great success in solving the continuity, momentum, and species balance equations which rely on the empirical estimation of relative permeabilities. This approach, however, may not be able to capture effects related to micro-scale structure in a multiphase system. On the contrary, pore-scale simulation can provide local detailed information about the flow field. The previous works on pore-scale simulation of immiscible two-phase flows in porous media [5][6][7][8][9][10] have been shown to be capable of capturing the effects of various factors, e.g., capillary number, viscosity ratio, surface wettability, and media heterogeneity. However, due to the inherent complexity related to interface evolution, heterogeneous geometry, dominant capillarity, and moving contact line, the computational simulation of multiphase displacement in realistic porous media still remains a research challenge.
As a mesoscopic method, the lattice Boltzmann method (LBM) has been widely accepted as a useful alternative for simulating multiphase flows, in particular, with the advantage of dealing with complex geometries such as porous media [11]. A number of multiphase LBM models have been proposed, e.g., the colour gradient model [12,13], the pseudo-potential model [14], the free energy model [15] and the mean-field model [16]. The colour gradient model is known to be the first multiphase model. Later, a new multiphase LBM model was proposed by Shan and Chen [14], where the interparticle forces (whether attractive or repulsive) mimicking the molecular potentials are introduced by a modified equilibrium velocity in the equilibrium distribution functions. The free energy model is thermodynamically consistent, and it considers a generalised equilibrium distribution function that includes a nonideal pressure tensor term. The original free energy model suffers from a lack of Galilean invariance, but later works have mitigated this problem [17,18]. The mean-field model describes interfacial dynamics based on the mean-field theory, and uses a pressure distribution function instead of its density counterpart to reduce discretization errors in the calculation of density gradient and to improve the numerical stability for variable density. Extensive reviews of these LBM multiphase models can be found in previous research [8,19].
The viscosity contrast is large while the density difference is small for the water or carbon dioxide involved oil recovery process. However, most LBM models can only access a limited range of viscosity ratios for dynamic fluid displacement. Yang et al. [20] performed a systematic comparison of three popular multiphase models, i.e., the pseudo-potential model, the free energy model, and the colour gradient model, for flow simulations in porous media. They concluded that the pseudo-potential model is a promising tool for liquid-gas systems, but may not be the optimal solution for the simulation of immiscible flows. Also, it produces high spurious velocities and thick interfaces for immiscible flows [20]. Jonathan et al. [21] used the pseudo-potential model to simulate two-phase flows in a two-dimensional channel with different viscosity ratios. A viscous-fingering-like phenomenon was observed for viscosity ratios of up to 10, and the simulations become unstable for larger viscosity ratios. Hatiboglu et al. [22] studied the spontaneous imbibition of water into a water-wet micro-model which was initially saturated with oil. It was verified that the maximum allowable viscosity ratio is four. The free energy model and colour gradient model show similar capabilities in modelling the flow of binary fluids with high viscosity contrast, but the latter produces a relatively thinner interface than the former one. In addition, the colour gradient model is able to control the interfacial diffusion and adjust the interfacial tension and viscosity independently to facilitate the numerical investigation.
The colour gradient LBM was first proposed by Rothman and Keller [13], and has been widely used for the simulation of immiscible fluids [9,12,[23][24][25][26]. To obtain an accurate and detailed understanding of the pore-scale mechanisms within porous media, Xu et al. [25] proposed a new algorithm for imposing the contact angle on the solid surface in the colour gradient LBM. The simulation results were validated against the experimental observations of Wu et al. [27]. However, these simulations only dealt with the simplified porous media, where the solid grains are regular circles and rectangles. In realistic porous media, the pore structures are extremely complex and irregular, which allow us to further assess the capability of numerical model. We apply the colour gradient LBM for dynamic displacement of immiscible fluids in a complex porous flow structure extracted from a realistic Berea sandstone sample, and study the effects of the viscosity ratio, interfacial tension, and contact angle on displacement process.

Colour Gradient Lattice Boltzmann Model for Two-Phase Flow
The colour gradient LBM was used in this work. In this model, the "Red" and "Blue" distribution functions ( f R i and f B i ) are used to represent two different fluids. The total distribution function is defined as Each colored distribution function undergoes the collision and streaming steps expressed by the following equation: where x and t are the position and time, k = red or blue, i represents the discrete velocity directions for the two-dimensional 9-velocity (D2Q9) model, ∆t is the time step, and c i denotes the lattice velocity vectors in the i-th discrete velocity direction. For the two-dimensional 9-velocity (D2Q9) model used in this work, the lattice velocity vector is given by The collision step includes self-and cross-interactions with the other type of particles which can be written as where (Ω k i ) (1) is the single-phase collision operator, (Ω k i ) (2) is the perturbation operator, and (Ω k i ) (3) is the recolouring operator.

Single-Phase Collision Operator
To achieve high accuracy and good stability, we implemented the collision step in accordance with the two-relaxation-time (TRT) model [28][29][30]-a special multiple relaxation time (MRT) model with only two relaxation rates [31]. With the TRT model, the single-phase collision operator is expressed as where M is the transformation matrix and is given by The diagonal relaxation matrix S is defined as For the conserved moments of density and momentum, the relaxation rates were set to 1. s e and s ν are related to the bulk and shear viscosities, and s ε and s q are free parameters. Following the recommendations of Ginzburg et al. [32], these relaxation rates were taken as where τ is the dimensionless relaxation time. The mass conservation for each fluid distribution function and the total momentum conservation require where ρ is the total density, and ρ k is the density of fluid k, and u is the local fluid velocity. The equilibrium distribution function ( f k,eq i (x, t)) is a truncated Taylor series expansion up to the second order in Mach number of the Maxwell-Boltzmann equilibrium distribution function, which was chosen to satisfy Equation (8): in which, c s = 1/ √ 3 is the lattice sound speed, and the weight coefficients are W 0 = 4/9, W 1−4 = 1/9 and W 5−8 = 1/36 for the D2Q9 lattice model. As we focused on the oil-water two-phase systems, the density ratio of both fluids was set to be unity.
To ensure a constant viscosity across the interface when both fluids were of different viscosities, the following harmonic mean [33] was employed to determine the viscosity of the fluid mixture: where η k (k = R or B) is the dynamic viscosity of fluid k, and η is the dynamic viscosity of the fluid mixture which is related to the dimensionless relaxation time τ by η = c 2 s ρ(τ − 0.5)δt; ρ N is the indicator function used to distinguish different fluids, and it is defined as

Perturbation Operator
(Ω k i ) (2) takes effect in the mixed interfacial region to generate an interfacial tension. Using the continuum surface force (CSF) [34] model, a body force term was added to the macroscopic momentum equation, and it reads as where σ is the interfacial tension coefficient, and K is the local interface curvature calculated by [34] where ∇ s = (I − nn) · ∇ is the surface gradient operator, I is the second-order identity tensor, and n = ∇ρ N / ∇ρ N is a unit vector normal to the interface. Such a body force was then incorporated into LBM through the forcing scheme proposed by Guo et al. [35] which can significantly improve computation accuracy and reduce spurious velocities effectively [36]. The perturbation operator (Ω k i ) (2) is expressed as with where α k is the fraction of interfacial tension contributed by fluid k and satisfies ∑ k α k = 1.
The fluid velocity was corrected to recover the Navier-Stokes (N-S) equations in the interfacial region where two fluids coexist

Recolouring Operator
The recolouring step was then applied to promote phase segregation and maintain a sharp interface [37] which is given as ( where f * i = ∑ k f k * i denotes the post-perturbation, pre-segregation value of the total distribution function along the i-th discrete velocity direction, and f eq i = ∑ k f k,eq i is the total equilibrium distribution function. β is the segregation parameter related to the interface thickness, and its value must be between 0 and 1 to ensure positive particle distributions-it was chosen to be 0.7 here. ϕ i is the angle between the indicator function gradient ∇ρ N and the lattice vector c i , which is defined by To minimize the discretization error, a fourth-order isotropic finite difference was used to calculate the partial derivatives. For example, for a variable ψ, the partial derivatives can be calculated by

Wetting Boundary Condition
For the colour gradient LBM, the most widely-used wetting boundary condition is to set a virtual density at solid surfaces [37]. However, it has been found in recent works [38,39] that numerical errors resulting from nonphysical mass transfer and spurious velocities may be accumulated during simulations which potentially render meaningless results. To overcome this, the algorithm proposed by Xu et al. [25] was implemented to impose the contact angle on solid surfaces. This algorithm is able to precisely control the contact angle for both static and dynamic problems which is essential for accurate simulations of two-phase displacement in porous media.

Relative Permeability Validation: Concurrent Flow in a 2D Channel
We first used the above model to simulate a pressure driven flow through a porous media. In the single-phase case, the intrinsic permeability K can be determined according to the Darcy's law [4] where u is the Darcy velocity in the porous medium. Note that the intrinsic permeability does not depend on the nature of the flow field but only on the geometry of the porous medium. In case of two-phase flows, the relative permeability K r,i can be given as [4] Here, i = w and nw refer to the wetting and non-wetting fluids respectively, µ i is the dynamic viscosity of fluid i, and ∇P i is the pressure gradient of fluid i.
First, we tested whether the LBM model can accurately predict relative permeability for a simple immiscible two-phase cocurrent flow. Specifically, we calculated the relative permeabilities of the layered two-phase flow in a channel under different saturation values. In the simulation, a constant pressure difference was applied on the inlet and outlet boundaries following Zou and He [40], and a non-slip (bounce-back) boundary condition was applied on the upper and lower walls using the halfway bounce-back scheme. The wetting phase flowed along the walls (a < |y| ≤ b) and non-wetting phase flowed in the central region (|y| < a). One can obtain the analytical solutions for the relative permeability, a function of wetting phase saturation S w and viscosity ratio λ [41], where S w is the wetting phase saturation, which is defined as S w = 1 − a/b, and λ is the viscosity ratio of the non-wetting to wetting fluids, i.e., λ = µ nw /µ w . For two typical viscosity ratios of λ = 1 and λ = 0.1, Figure 1 shows excellent agreement between the computed relative permeabilities and the analytical solutions.

Problem Statement
The 2D micromodel is shown in Figure 2, and its porosity is 0.33. The entire micromodel is divided into L × H = 2644 × 2117 lu 2 with a resolution of 0.67 µm per lu. The narrowest throat has a width of 12 lu which is fine enough to produce grid-independent results in two-phase displacement simulations [7].Twenty-two extra layers of lattices are added to both the inlet and outlet to facilitate the implementation of boundary conditions, so the actual domain size is 2688 × 2117 lu 2 . Here, the flow velocity for each fluid is so small that the inertial effect becomes insignificant.
The constant pressure boundary condition developed by Zou and He [40] was imposed on the inlet and outlet, which are located at the left and right sides of the domain, respectively. The halfway bounce-back scheme and wetting boundary condition were applied at the walls that contain the surfaces of solid grains as well as the upper and lower sides of the domain. The densities of both wetting and non-wetting fluids were set to be unity as the density difference between oil and water is insignificant. The dynamic viscosity of the displaced fluid was kept as 0.167 mu·lu −1 ·ts −1 , and the dynamic viscosity of the displacing fluid was varied to achieve different viscosity ratio values (λ = µ displaced /µ invading ). Initially, for x ≥ 22 lu, the domain was filled with the displaced fluid (represented in red) at rest, and a constant pressure p out = 1/3 mu·lu −1 ·ts −2 was imposed at the outlet, while for x < 22 lu, the region was occupied by the displacing fluid (represented in blue) with a higher pressure at the inlet. The simulations stopped when the displacing fluid broke through at the outlet to avoid the capillary end effect [43]. By studying the dependency of porosity on the size of the rock sample, we calculated the porosity using different sizes of Berea sandstone to find its representative elementary volume (REV). The results are shown in Figure 3, which demonstrates that the chosen size of 2688 × 2117 lu 2 could be a good REV of the Berea sandstone. This guarantees that the simulated geometry included a sufficient number of pores for the meaningful statistical average, which is generally required by the continuum concept [4].

The Viscosity Ratio Effect
The viscosity ratio is one of the most important dimensionless numbers in two-phase displacement through a porous media, and its effect was studied for a constant contact angle of 120 • , corresponding to a weak wetting property of the displacing fluid. The inlet pressure is also the same for different viscosity ratios, i.e., p in = (p out + 0.0125) mu·lu −1 ·ts −2 . Figure 4 shows the fluid distributions inside the pore network when breakthrough of the displacing fluid occurs for viscosity ratios of 1, 2, 5, and 10. It can be seen that for all the viscosity ratios considered, two or three preferential paths are generally formed. As the viscosity ratio increases, it becomes harder for the finger to flow into the large pore as highlighted by the black elliptical circles. The front of the bottom finger is finally halted at the position near the inlet for the viscosity ratio of 10. Also, the breakthrough finger switches to a lower position at the outlet when the viscosity ratio increases to 10.  Figure 5 shows the time evolution of the displacing (wetting) fluid saturation for a viscosity ratio changing from 1 to 30. It can be clearly seen that a higher viscosity ratio results in a lower fluid saturation which is less favourable for the recovery of displaced fluid. This is consistent with previous experimental [44] and numerical results [45]. At the viscosity ratio of unity, the highest displacing fluid saturation (around 0.37) is observed, and the displacing fluid saturation keeps increasing linearly with time. Although the displacing fluid saturation reduces with λ, the reduction becomes less significant at larger λ values. We also tested whether the above observations are valid for varying wettability and interfacial tension. Four pairs of contact angle and interfacial tension were considered (see Table 1), and the simulation results are shown in Figure 6. It is clear that the displacing fluid saturation drops rapidly when the viscosity ratio increases from 1 to 20 and then stays almost unchanged for all the cases considered. For each viscosity ratio, increasing contact angle or interfacial tension always leads to a higher displacing fluid saturation or displacement efficiency. In addition, the limitation of the current model was tested by decreasing the viscosity of the displacing fluid, and it was found that the maximum viscosity ratio that it can successfully access is as high as 50.

The Interfacial Tension Effect
The interfacial tension between the displacing and displaced fluid is known to play an important role at the pore scale for two-phase displacement through a porous media. The same viscosity ratio (λ = 1) was used and the contact angle was fixed at 120 • . Figure 7 shows the fluid distributions at breakthrough for the interfacial tensions of 0.03 mu·ts −2 , 0.015 mu·ts −2 , and 0.0005 mu·ts −2 . From this figure, it is seen that when the interfacial tension is as low as 0.0005 mu·ts −2 , the viscous force prevails, the main fingers get thinner and even break up, and many fingers only occupy partial pores or throats that they pass by. Also a number of small blobs of displaced fluid are left behind at low levels of interfacial tension. When the interfacial tension increases to 0.015 mu·ts −2 , larger blobs of the displaced fluid are trapped by the displacing fluid, and the displacing fluid flows into the large pore more easily due to higher capillary pressure, as highlighted by the black elliptical circles in Figure 7. As the interfacial tension is further increased to 0.03 mu·ts −2 , we surprisingly found that the front of the displacing fluid can be flat inside some pores, e.g., the one surrounded by the black rectangle box in Figure 7, though most interfaces remain curved due to large capillary pressure. Meanwhile, much less of the displaced phase is left behind as the invading phase continues to move forward.  Figure 8 shows the evolution of the displacing fluid saturation for an interfacial tension of 0.0005 mu·ts −2 , 0.015 mu·ts −2 and 0.03 mu·ts −2 . It is seen that generally higher interfacial tension results in a higher saturation. The linear relationship between saturation and the evolving time is again observed.

The Contact Angle Effect
The contact angle was measured from the displaced fluid. All the contact angles were larger than 90 • in this part to simulate the imbibition process. Figure 9 shows the fluid distributions when the displacing fluid breaks through at the outlet for the contact angles of 162 • , 150 • , 135 • , and 120 • . It is observed that the invading paths of the wetting fluid do not change significantly except for a few branches. An obvious difference can be seen at the contact angle of 135 • where the bottom finger advances furthermost. Figure 10 shows the evolution of the displacing fluid for different contact angles. The case of θ = 120 • has the lowest wetting fluid saturation but its difference from the others is very small. It can also be seen from Figure 6 that the contact angle does not have a big impact on the saturation. Therefore, the contact angle has a limited effect on the efficiency of displacement in the micro-model of Berea sanstone.

Conclusions
The colour gradient two-phase LBM was used to study dynamic displacement of immiscible fluids in a 2D micromodel of Berea sandstone. The effects of the viscosity ratio, interfacial tension, and contact angle on the fluid distributions at breakthrough and the evolution of displacing fluid saturation were systematically investigated. When the viscosity ratio is no more than 20, a higher viscosity ratio results in a lower displacing fluid saturation which is less favourable for the recovery of the displaced fluid. However, when the viscosity ratio is larger than 20, the saturation is almost constant for both imbibition and drainage, regardless of the interfacial tension. At the viscosity ratio of unity, the displacing fluid saturation has the highest value, and it linearly increases with time.
The interfacial tension affects the flow pattern. When the interfacial tension is as low as 0.0005 mu·ts −2 , thin viscous fingers appear. A number of small drops of displaced fluid are left behind. As the interfacial tension increases from 0.0005 mu·ts −2 to 0.015 mu·ts −2 , the size of the trapped blobs increases, and the number decreases. The displacing fluid flows into large pores more easily due to the increased capillary pressure. When the interfacial tension continues to increase, smaller amount of the displaced phase is left behind as the invading phase continues to move forward. In addition, increasing interfacial tension results in a higher saturation of wetting fluid.
The contact angle generally does not change preferential flow paths except for a few branches in the imbibition, suggesting that the contact angle has limited effect on the efficiency of displacement within 2D Berea sandstone.