The CFD method for multiphase flow describes the mechanism of multiphase fluid flow through the Navier-Stokes equations, while simultaneously coupling convection-diffusion equations regarding phase variables to achieve interface tracking, such as Level Set method, Phase Field method, and Volume of Fluid method, etc. Many studies have shown that the Phase Field method is more suitable for simulating non-miscible oil-water two-phase flow in porous media. Essentially, it remains a continuous method. This method describes the two-phase mixing interface as a diffusion interface layer with a certain thickness. On this interface layer, the phase variable gradually transitions from the oil phase to the water phase. The evolution equation of the phase variable is the Cahn-Hilliard equation derived from the free energy, describing the phase separation process under the influence of convection and diffusion: δ ϕ δ t + u ⋅ ∇ ϕ = ∇ ⋅ γ λ ε 2 ∇ ψ ψ = − ∇ ⋅ ε 2 ∇ ϕ + ϕ 2 − 1 ϕ

In the equation, ϕ represents the dimensionless phase field variable; γ denotes the mobility, reflecting the influence of diffusion on the interface, in m²·s/kg; λ is the energy density of mixing, in N; ε is a coefficient related to the interface, in m; ψ represents the phase field auxiliary variable.

Phase field method is an incompressible two-phase flow interface diffusion model, suitable for mixed-phase, non-miscible, and partially miscible flows. The phase field model identifies different phase regions and transition regions between phases by introducing a phase field variable ϕ . In homogeneous phase regions, the phase field variable ϕ is constant (typically set to −1 and 1), while in the transition regions at interfaces, the phase field variable ϕ varies continuously (between −1 and 1).

The fluid motion equations consist of the mass and momentum conservation equations, which are the continuity equation and the Navier-Stokes equations, respectively. ∇ × u = 0 ρ δ u δ t + ρ ∇ × u u = ∇ − p I + K + F s t

In the equation, ρ represents the fluid density, in kg/m³; u is the velocity vector, in m/s; p denotes pressure, in Pa; I symbolizes the unit vector; K is the viscous stress tensor, in Pa; F s t is the interfacial tension acting on the fluid interface, in N·m^-3.

The viscous stress tensor is defined by the following expression: K = μ ∇ u + ∇ u T

The model presented in this article is based on the built-in phase field two-phase flow basic model in Comsol. Its solution method relies on finite element method. Additionally, in order to achieve rapid solution of high-quality grids under complex pore structures, this article introduces the density and viscosity equivalent models in the mixed-phase region, which can be calculated by the following equations: ρ = 1 + ϕ ρ 1 2 + 1 − ϕ ρ 2 2 μ = 1 + ϕ μ 1 2 + 1 − ϕ μ 2 2

Where, Subscripts 1 and 2 represent fluid one and fluid 2. In the phase field method, F s t can be represented as: F s t = G ∇ ϕ

Where, G represents chemical energy, J·m^-3, which can be expressed as G = λ − ∇ 2 ϕ + ϕ ϕ 2 − 1 ε 2 = λ ε 2 ψ

It can be seen that utilizing the dimensionless phase field variable ϕ and the phase field auxiliary variable ψ , surface tension can be computed as a distributed force on the interface. Therefore, coupling the Cahn-Hilliard equation and the Navier-Stokes equation avoids the complex process of calculating surface tension using surface normals and surface curvature.

The solid boundary condition for the motion equations is a no-slip boundary, and the phase variable does not cross the solid boundary, which is described as: u n = 0 n ⋅ γ λ ε 2 ∇ ψ = 0 n ⋅ γ λ ε 2 ∇ ϕ = ε 2 ⁡ cos ⁡ θ ∇ ϕ

In the equation, u n represents the fluid velocity at the solid wall, in m/s; n denotes the unit normal vector to the solid wall; θ is the contact angle, in rad.

The visualization experiment used to compare the numerical simulations is sourced from Karadimitriou et al. (2012, 2013). The experimental parameters are fully characterized and the model is easy to accurately reconstruct. Therefore, this reference selects its experimental results for comparison. At the same time, this experiment is used for verification in many two-phase flow numerical simulations, such as SPH et al. (Kunz et al., 2016).

Karadimitriou et al. (2013) studied the displacement flow behavior of internal fluorinated liquids and water using visualized microscopic models. This paper compares the established phase-field method model with their experimental results to validate the numerical simulation method. The microscopic model used for experiments is a rectangle with dimensions of 2.5×1 mm. It consists of several circular pillars of varying sizes representing the skeletal particles of the porous medium, with an average pore diameter of 160 µm. The fabrication involves photolithography on a silicon wafer, followed by covering the microscopic model with a PDMS sheet and bonding. The experimental setup consists of three main parts: an optical acquisition system, a pressure sensor, and a pumping system. More details about the experimental setup can be found in the work of Karadimitriou et al. (2013). During the experiment, the initial state fills the geometric model with water as the displaced phase, and the walls of the microscopic model are hydrophilic, then fluorinated liquid is injected for displacement. The microscopic model has a contact angle of 45°, a density of 1800 kg/m³, a viscosity of 4.7 mPa·s for the fluorinated liquid, a density of 1,000 kg/m³, a viscosity of 1 mPa·s for water, and a surface tension of 0.055 N·m⁻¹, with a displacement pressure difference of 1860 Pa. The basic parameters of the verification model are consistent with the experiment, and the results show that most of the residual displaced water phases obtained from the experiment and simulation are similar, indicating that the phase-field method can largely simulate the distribution state of two-phase fluids (Figure 1).

Figure 4 presents the distribution of remaining oil at different displacement moments, where white represents rock particles, red represents the oil phase, and blue represents the water phase. Initially, when water is injected, the interface advances relatively uniformly. Subsequently, the water entering from each inlet rapidly flows outward in search of exit channels. Once an exit channel is found in any direction, a continuous flow forms in that direction, overcoming the frictional resistance of the oil phase. Water in other directions might undergo backflow and then enter the established flow channels, leading to convergence. Additionally, when water passes through the throats, capillary valve effects and capillary fingering effects occur. At the entrance of a throat, capillary force acts as a resistance, making it difficult for water to pass and even resulting in Jamin effect; when water passes through the throat, capillary force acts as a driving force. The displacement pressure difference and strong capillary dynamics facilitate flow, overcoming the inhibitory effect of viscous resistance, allowing water to pass quickly and be easily segmented into discrete droplets. In most simulations in this study, the rock matrix is assumed to be oil-wet; therefore, the wetting phase easily forms film flow on the wall, and the oil film on the throat wall interrupts the water phase.

Ultimately, six types of remaining oil distributions are observed. The distribution patterns of remaining oil can be categorized into dead-end remaining oil, pressure balance remaining oil, wall-bound remaining oil, Jamin effect remaining oil, water-encapsulating remaining oil, and remaining oil in the direction perpendicular to displacement. Their proportions and scales are controlled by the micro-pore structure, wettability, and dominant channels.

The formation of dead-end remaining oil is mainly due to the lack of expulsion channels for the oil phase, and its distribution depends on the pore structure, being relatively dispersed (Figure 4A). In oil-wet conditions, dead-end remaining oil primarily results from the inability of the water phase to enter; in water-wet conditions, under the action of capillary forces, the water phase invades along the walls into the dead ends, separating the oil phase in the dead ends. Thus, in oil-wet conditions, there are basically no effective methods to enhance recovery of dead-end remaining oil; in water-wet conditions, reducing the displacement speed allows the water phase to infiltrate into the dead ends under the action of capillary forces, gradually expelling the oil phase within the dead ends.

Pressure balance remaining oil refers to the scenario during the displacement process where the pressure at both ends of the pore throats is equal, failing to establish a positive pressure differential, resulting in the oil phase within the pore throats being unable to be effectively expelled due to a lack of driving force (Figure 4B). Similarly, there is residual oil in the direction perpendicular to displacement, where the pressure gradient across the pore throats is also minimal, and due to the velocity direction of the water phase being perpendicular to the direction of the pore throats during displacement, it is difficult for the water phase to enter these pore throats and displace the oil phase.

Wall-bound residual oil refers to the condition in an oil-wetting state where the oil phase on the surfaces of pore throats experiences significant capillary resistance (Figure 4C). In situations of low displacement velocity, the driving force is insufficient to overcome the wall capillary resistance, resulting in a thin layer of oil phase remaining on the surfaces.

The Jamin effect residual oil refers to a condition in an oil-wetting state where the volume of oil droplets is large, and when passing through finer throats, they encounter significant capillary resistance (Figure 4D). Despite repeated deformation, the oil droplets are unable to pass through and return to the previous pore space, forming isolated oil droplets.

Water-encapsulated residual oil refers to the phenomenon during the displacement process in complex porous media where, due to the uneven velocity within the water phase, after some of the water phase passes through the throats, the remaining oil phase in the pores is squeezed out and rapidly occupies the position previously held by the water phase, resulting in the formation of isolated water droplets and the occurrence of the water encapsulation phenomenon (Figure 4E).

The residual vuggy oil refers to a scenario where the vugs are relatively large, facilitating rapid advancement of the water phase through the central region of the vug, creating preferential pathways. Consequently, due to the influence of surface tension, a significant amount of residual oil gets trapped on the vug walls, making it challenging to displace (Figure 4F).

In this paper, a constant flow rate displacement is set, thus the water flooding efficiency is divided into three stages. The first stage is the linear growth phase before the water phase reaches the outlet; the second stage occurs when the water phase has reached the outlet, but the oil phase in the model has not yet equilibrated; the third stage is the equilibrium phase, where the displacement efficiency no longer changes, and the crude oil remaining in the pores can no longer be displaced. In the simulation process, the length and slope of the first stage are related to the flow rate, while the occurrence of the second and third stages is related to the number of dominant percolation channels. With only one dominant percolation channel, there appears to be one plateau; when there are multiple dominant percolation channels, there are multiple plateaus. The reason for the plateau in stage3 is due to the water phase reaching the outlet and the formation of dominant migration channels, at which point the water phase on other paths recirculates to the dominant migration channels. In this large pore and coarse throat model, when the wettability is oil-wet and the CA is 30°, the final displacement efficiency is close to 37%, with most of the crude oil far from the dominant percolation channels remaining untouched (Figure 5).

(1) Influence of Throat Size

Figure 6 shows the oil-water distribution and microscopic oil removal efficiency changes under different throat sizes, where red, blue, and white represent the oil phase, water phase, and rock matrix, respectively. The transition zone between red and blue is the oil-water interface. It can be observed that when the pore structure consists of large pores and coarse throats, at a timestep of 0.02 s, water enters from each inlet, with three injection points at the upper inlet converging into one flow channel, and the two injection points at the lower part gradually converging and branching. Then, the water phase begins to seek the path of least resistance, where the front end of the water phase shows finger-like features, even developing branches perpendicular to the direction of injection velocity. At a timestep of 0.04 s, the first batch of water has reached the outlet, and a dominant flow channel has been formed. The pressure gradient in other flow channels decreases, and the early branching water seeking paths gradually retract back to the flow channel, with some of the water that did not retract in time forming isolated droplets. Overall, in a water-wet state, besides the dominant flow channel, it is difficult to mobilize the oil phase in other locations, with the main types of remaining oil being pressure equilibrium residual oil, followed by dead-end residual oil and water-encapsulating residual oil.

When the pore structure consists of large pores and fine throats, the water phase migration behaves similarly to that in structures with large pores and coarse throats, exhibiting characteristics such as bifurcation, convergence, and contraction of the water phase front. However, overall, the water phase migration path becomes more tortuous and complex, leading to potentially more throats being affected. Nevertheless, due to the smaller size of the throats, capillary resistance increases, necessitating higher breakthrough pressures. In a time step of 0.5 s, the majority of the crude oil remains unused, with the predominant type of remaining oil being pressure equilibrium residual oil, followed by dead-end residual oil, pore residual oil, and water-encapsulated residual oil (Figure 7).

Comparing the microscopic displacement efficiency of structures with large pores and coarse throats, it is observed that the overall shape of the curves is similar and can be divided into three stages. The quantities in stage 2 and stage 3 are the same because ultimately only one dominant migration pathway exists. Despite the finer throats in structures with large pores and fine throats and the more tortuous water phase displacement routes, they reach the outlet earlier than structures with large pores and coarse throats. From the perspective of final oil displacement efficiency, the oil displacement efficiency for structures with large pores and coarse throats is 37%, while for structures with large pores and fine throats, it is 40%. This implies that smaller pore/throat sizes lead to higher oil displacement efficiency, contrary to traditional understanding. This discrepancy may be related to the results of the constant injection velocity displacement discussed in this study. The constant injection velocity leads to the generation of pressure building, resulting in higher pressure differentials within Type B. This will be further discussed in subsequent analyses.

(2) Influence of Wettability

As indicated in section 4.3.1, the efficiency of oil displacement depends on whether the oil-water interface advances uniformly. The movement of the oil-water interface is influenced by various forces, including inertia, displacement pressure, capillary force, wet phase viscosity, non-wet phase viscosity, gravity, and inlet effects. For oil-water flow in porous media, inertia and gravity effects can be neglected at low Reynolds numbers. Assuming the case of a single capillary as a straight cylindrical shape, the capillary force below the single capillary can be derived from the Laplace-Young equation: p c = 2 σ ⁡ cos ⁡ θ R

Therefore, in the hydrophilic state (CA 0°–90°), the capillary force is positive, manifesting as driving force; whereas, when wetting transitions to hydrophobic (CA 90°–180°), the capillary force becomes negative, representing resistance. Capillary force acts as the driving force in water-wet media, while in neutrally wet or oil-wet media, it presents resistance due to dynamic contact angle. Therefore, water-wet media exhibit higher final oil recovery efficiency, and the microscopic oil displacement efficiency in neutrally wet or oil-wet media is more sensitive to water phase breakthrough.

Figure 8 shows the distribution of remaining oil as the wettability gradually changes from strongly oil-wet to strongly water-wet in the large-pore fine-throat model under an injection velocity of 0.01 m/s. Comparing Figure 6 and Figure 8, it can be observed that when the system is strongly oil-wet, the main driving force for displacement is the displacement pressure differential. The water phase tends to spread along the dominant migration pathway, and the remaining oil mainly consists of pressure equilibrium residual oil, with the majority of the remaining oil in the pores unused. When the system exhibits weak oil-wettability, the displacement pressure remains the primary driving force, but the capillary resistance decreases significantly. Consequently, the area covered by the water phase increases substantially, and the remaining oil mainly consists of pressure equilibrium and dead-end types. When the system reaches a neutral wettability state, the capillary resistance further diminishes, leading to a slight increase in the area covered by the water phase, with the remaining oil predominantly comprising pressure equilibrium and dead-end types. In the water-wet state, the main driving forces are displacement pressure and capillary force, resulting in a larger area covered by the water phase, with dead-end residual oil being the primary type of remaining oil.

From the perspective of oil recovery efficiency, due to capillary force being the driving force in water-wet media, the ultimate oil recovery efficiency is higher in water-wet media. In an oil-wet reservoir with CA 150°, the final oil recovery efficiency is 40%, while in a water-wet reservoir with CA 30°, the final oil recovery efficiency is 75%. Regarding the change in curves, the shape of the curve in Stage 1 remains consistent, and the time taken for the water phase to reach the outlet is essentially the same. However, there is a significant difference in the shape of Stage 2. In strongly oil-wet reservoirs, Stage 2 appears as a single step, and due to water phase refluxing to the dominant migration pathway, there is a noticeable drop between Stage 2 and Stage 3. As the reservoir transitions to water-wet conditions, Stage 2 exhibits distinct multiple slope characteristics, indicating the opening of new migration pathways connecting to the outlet during the displacement process (Figure 9).

(3) Influence of Injection Velocity

Under the strong displacement force conditions of water flooding, the displacement force is more significant than capillary force, and the viscosity in small pores is stronger. Therefore, the leading edge propagation velocity in large pores is greater than in small pores, thereby trapping oil in the small pores. Conversely, under hydrophilic conditions, in the weak displacement force conditions of water flooding, capillary force plays a more significant role, with clear capillary advancement, often showing faster advancement in small pores and slower advancement in large pores. However, in reservoirs with small throats, the associated viscous resistance is also stronger, even exceeding capillary force, so capillary advancement does not always occur. To achieve a more uniformly advancing water flooding interface, it is necessary to comprehensively analyze the evolution of water flooding oil based on the combined characteristics of capillary diameter, displacement rate, and wettability, thereby improving reservoir development efficiency.

However, under oil-wet conditions, capillary forces act as resistance. Therefore, when the displacement velocity increases, the displacement pressure can overcome capillary resistance and viscous resistance, potentially displacing the oil phase into more fine throats. Figure 10 illustrates the oil-water distribution under different displacement velocities. It can be observed that when the displacement velocity increases from 0.01 m/s to 0.02 m/s, the isolated water phase originally produced by backflow in the lower-left corner is broken through, adding a new migration pathway, but the overall difference is not significant. As the displacement velocity continues to increase to 0.05 m/s and 0.1 m/s, another migration pathway is added, and the water phase coverage area increases slightly, but the increase is not significant. Although the oil recovery efficiency increases, the predominant type of remaining oil remains pressure balance.

Figure 11 compares the microscopic oil displacement efficiency under different displacement velocities. It can be observed that under oil-wet conditions, as the displacement velocity increases, the water phase reaches the outlet earliest and the equilibrium section appears earliest. When the displacement velocity is 0.1 m/s, the characteristics of Stage 2 are not prominent. From the perspective of final oil recovery efficiency, as the displacement velocity increases, the final oil recovery efficiency increases. When the displacement velocity increases from 0.01 m/s to 0.1 m/s, the oil recovery efficiency increases from 40% to 55%. In other words, in oil-wet carbonate reservoirs, higher speed water flooding is more suitable.

(4) Characteristics of Water Flooding in Fractured and Fractured-Vuggy Reservoirs

Under strongly oil-wet conditions, although the fracture width distribution in fractured reservoirs is relatively uniform compared to pore-type reservoirs, the water phase still moves rapidly along the main fracture towards the outlet. The water in the upper fractures advances approximately 0.2 mm forward before the water phase reaches the outlet, but as the dominant migration pathway through the lower main fracture is established, the water phase in the upper fractures communicates through the vertical fractures and no longer advances forward. In fractured-vuggy reservoirs, the overall migration pathway is similar to that in fractured reservoirs, but most of the oil in the pores on the dominant migration pathway is displaced, and the oil away from the dominant migration pathway cannot be mobilized (Figure 12).

From a microscopic perspective, the oil displacement efficiency of fractured-pore reservoirs (33%) is slightly higher than that of fractured reservoirs (30%). Comparing with pore-type reservoirs, it is found that the presence of Stage 2 is almost indistinguishable in fractured and fractured-pore reservoirs. This is because the capillaries in fractured reservoirs are relatively wide, resulting in lower capillary resistance and the absence of backflow situations (Figure 13). The early high oil displacement efficiency in fractured reservoirs may not result from involving a larger pore space but could be related to the pore space area. Fracture-pore type reservoirs have larger pore area, potentially leading to a relatively smaller oil displacement efficiency ratio, with a larger area occupied by the water phase.

(5) Influence of Pore structure

The pore structures within the study area reservoirs are significantly more complex than those described in the previous model. These reservoirs typically include cast mould pores (Type E), visceral foramen (Type F), primary intergranular pores (Type G), and intergranular dissolution pores (Type H), each exhibiting substantial differences in pore structure and consequent oil displacement efficiency.

Firstly, under constant velocity boundary conditions (0.05 m/s) and CA 90°, the water displacement efficiency of different pore structures is compared. In mould pores, most of the oil phase remains unrecovered, especially those distant from dominant migration pathways, with remaining oil primarily of the pressure equilibrium type. In visceral foramen, most remaining oil is concentrated within biological pore, with significant residual oil in almost every large pore, resulting in relatively higher overall displacement efficiency, primarily pressure equilibrium residual oil. In primary intergranular pores, there is minimal remaining oil, as most crude oil is effectively displaced, which is attributed to their homogeneity, with remaining oil mainly comprising pressure equilibrium and dead-end residual oil. In intergranular dissolution pores, remaining oil is also minimal, with most crude oil effectively displaced and primarily concentrated within dissolution pores, mainly comprising pressure equilibrium residual oil (Figure 14).

Contrary to traditional beliefs, smaller throat/pore sizes exhibit higher oil displacement efficiency, possibly influenced by the results of constant flow displacement, which generates backpressure. Therefore, this paper discusses the pressure drop conditions under constant velocity displacement in different types of pore structures (Figure 15; Figure 16). It is evident that while each pore structure type exhibits driving capabilities, mould pores demonstrate a backpressure reaching 12 MPa, increasing to 70 MPa under oleophilic conditions (CA120°), rendering displacement impractical for actual reservoir water flooding processes, especially when CA reaches 150°, making displacement impossible. In visceral foramen, as fluid exits, pressure differentials reach 1.6 MPa, increasing nearly tenfold when wetting transitions from neutral to strongly oleophilic. However, in primary intergranular pores and intergranular dissolution pores, pressure drops remain below 0.1 MPa under CA90°, increasing to only 0.5 MPa under strongly oleophilic conditions (CA30°).

Next, the injection condition is altered to constant pressure displacement (0.2 MPa). It can be observed that under neutral conditions, crude oil in mould pores remains completely immobile, with capillary resistance far exceeding the displacement pressure. In visceral foramen, most crude oil within biological pores remains stagnant, with the aqueous phase only infiltrating along dominant pathways. In primary intergranular pores, there is minimal remaining oil, as most crude oil is effectively displaced, primarily comprising pressure equilibrium residual oil. In intergranular dissolution pores, the water displacement situation is inferior to that of intergranular pores, with significant residual oil remaining within numerous dissolution pores, primarily comprising pressure equilibrium residual oil (Figure 17).

Finally, a comprehensive comparison was made between the oil displacement efficiency of different types of pore structures under constant velocity displacement and constant pressure displacement. From the perspective of actual reservoir displacement efficiency, intergranular pores > intergranular dissolved pores > visceral foramen > mould pore, with this trend being more pronounced under oleophilic wetting conditions (Figure 18).