Pore-Scale Simulations to Enhance Development Strategies in Offshore Weak Water-Drive Reservoirs

1 Introduction

Offshore weak water-drive reservoirs are rich in hydrocarbon resources; however, its overall productivity remains low, and its economic efficiency is poor, significantly limiting large-scale development [1,2]. Weakly water-driven reservoirs lack sufficient aquifer support to maintain pressure, so sustained production relies on artificial energy supplementation to improve recovery. Such reservoirs commonly exhibit poor reservoir quality, strong heterogeneity, and complex pore–throat structures that compartmentalizes flow. Analysis has identified the following key production challenges: (1) The long-term reliance on natural energy development has resulted in severely insufficient formation energy. It is urgent to supplement energy through water or gas injection to restore reservoir pressure and enhance oil recovery. However, the current injection-production schemes lack systematic optimization of key parameters, which negatively impacts the energy supplement effect and displacement efficiency. (2) The viscosity of crude oil varies significantly both horizontally and vertically, with a range of 1 mPa·s–8 mPa·s. This variability leads to difficulties in controlling the mobility ratio during displacement, presenting significant challenges for optimizing injection and production parameters in waterflooding development. (3) The complex fluid properties coupled with weak water energy, result in limited effectiveness of conventional waterflooding, necessitating the exploration of more efficient development methods to enhance both recovery factor and economic benefits. In response to the above-mentioned problem of insufficient understanding of the seepage mechanism and development methods, this study conducts a microscopic flow simulation analysis of weak water-drive reservoirs. The aim is to propose effective enhanced oil recovery (EOR) strategies, thereby providing a scientific basis for formulating development technical policies and ultimately achieving efficient reserve utilization and improved recovery [3,4,5].

Currently, the optimization of waterflooding development and the selection of development strategies lack scientific and rational guidance from reservoir engineering methodologies. The optimization and adjustment of water and gas injection parameters are primarily qualitative, with parameter design lacking a quantitative theoretical foundation [6]. Although numerical simulation techniques have become increasingly advanced, field application design still relies heavily on analogy and extrapolation from specific field experiences. This approach is not based on quantitative optimization algorithms derived from the petrophysical characteristics and fluid interactions of the reservoir itself. This leads to significant uncertainties in enhancing oil recovery (EOR) efforts [7,8,9]. In particular, the dynamic behavior of water/gas injection fronts remains an under-explored area, and systematic understanding of residual oil distribution in low-permeability reservoirs is still insufficient [10,11]. Zhou et al. [12] employed a systematic approach combining nuclear magnetic resonance (NMR) tests with fluorescence analysis of frozen slides to evaluate the effects of alkaline, surfactant, polymer, and ASP flooding systems on oil recovery efficiency and the pore-scale distribution of residual oil. They focused on analyzing the variation characteristics of residual oil saturation across different pore size ranges. Li et al. [13] constructed a Multicomponent and multiphasemulti relaxation time-Lattice Boltzmann model (MCMP-MRT-LBM) to investigate the occurrence states of crude oil in laminated shale under varying water saturation conditions. The results demonstrated that with the increase in water content, the remaining oil gradually changes from continuous state to discontinuous state. The remaining oil was gradually distributed in quartz in columnar and solitary droplets and distributed in kerogen in a membranous form. Water is gradually distributed in clusters throughout the reservoir.

The pore-scale multiphase flow simulation techniques primarily include pore network modeling (PNM) [14,15,16], the lattice Boltzmann method (LBM) [17,18], and VOF [19,20,21] based on direct numerical simulation (DNS), the level set method (LS) [22,23,24], and the phase-field method (PF) [25,26]. In addition to physics-based numerical modeling, data-driven approaches such as Response Surface Methodology (RSM), machine learning (ML), and artificial neural networks (ANN) have also been employed in recent years to predict permeability evolution and optimize flooding performance at larger scales [27]. Among these methods, PNM simplifies the modeling of complex pore structures but has limitations in accuracy. The LBM is prone to numerical instability when simulating multiphase systems with high density ratios and significant viscosity differences. In comparison, direct numerical simulation methods exhibit superior numerical stability. Specifically, the VOF is characterized by accurate interface capturing and computational efficiency, making it highly suitable for simulating complex flow processes in porous media. Huang et al. [28] utilized the interFoam solver in OpenFOAM to simulate the water flooding process, focusing on the effects of varying flow velocity. The study investigated how changes in flow velocity influence the oil-water distribution, residual oil location, and morphology under constant water injection volumes. However, the study was limited to a two-dimensional model and did not account for three-dimensional variations. SONG et al. [29] incorporated the interfacial tension term and developed a pore-scale mathematical model for CO2 flooding using an open-source solver. The results indicated that the recovery rate of miscible CO2-SAG flooding was 16.5% higher than that of continuous CO2 injection. However, the study did not consider the impact of water-gas displacement on the swept volume.

This study takes a typical offshore weak water-drive reservoir as the research object. Based on CT scanning and high-precision digital core reconstruction, a geometric model was established, and a pore-scale oil–water (gas) flow mathematical model was developed. The VOF method is then used to simulate waterflooding optimization and development method selection. By quantitatively analyzing the effects of injection velocity, oil-water viscosity, miscibility degree, and injection patterns on oil recovery efficiency and front stability, the micromechanical mechanism of development parameter optimization is clarified. This study aims to provide a theoretical basis and quantitative design for transitioning development methods and injection-production adjustment in offshore weak-water drive reservoirs.

2 Mathematical Methods and Geometric Modeling

2.1 Mathematical Formulation

The Navier-Stokes equations are a set of partial differential equations used to describe the motion of viscous fluids and are fundamental to fluid mechanics. They account for the effects of fluid viscosity, pressure gradients, inertial forces, and external forces, enabling accurate simulations of fluid behavior under different conditions [30,31]. By coupling the continuity equation with the N-S equations, the pressure and velocity of incompressible two-phase flow in Newtonian fluids can be solved. The governing equations are as follows: ∇⋅u=0(1) ∂ρu∂t+∇⋅ρuu=−∇p+∇⋅μ∇u+∇Tu+Fsa(2) where u is the velocity vector, in m/s; p is the pressure, in Pa; ρ is the density, in kg/m3; μ is the viscosity, in mPa·s; and Fsa is the interfacial force acting on the phase interface, in N/m.

The core idea of the VOF method is to define the volume fraction of the reference phase within a control volume. If the control volume is entirely filled with the water phase, the volume fraction is 1; if it is entirely filled with the oil phase, the volume fraction is 0. At the interface, the volume fraction within the control volume lies between 0 and 1, as shown in Fig. 1 [32].

Figure 1: Volume fraction schematic.

To maintain the mass and momentum balance of each phase, the density ρ and viscosity μ in Eqs. (3) and (4) can be expressed as functions of the α: ρ=αρ1+1−αρ2(3) μ=αμ1+1−αμ2(4) where subscripts 1 and 2 represent the water phase and the oil phase, respectively.

By substituting the density and viscosity calculated in the two-phase interface region into the continuity Eq. (1), the phase equilibrium equation for the two-phase fluid flow can be derived:

The definition of the interfacial force Fsa in Eqs. (2) is given by [33,34]: Fsa=∫τσkδx−xsn^dΓxs(6) where Γ represents the two-phase interface, δ x − x s is the Dirac delta function acting on the two-phase interface, k is the interface curvature, n ^ is the unit vector, and σ is the interfacial tension between the two phases, in N/m. The expression for the interface curvature k is as follows:

The contact angle θ satisfies: n^⋅n^s=cosθ(8) where n ^ is the unit vector perpendicular to the interface between the two phases, and n ^ s is the unit vector perpendicular to the wall.

Microscopic flow simulation analysis was carried out by pore-scale CFD on the digital core. The governing equations were implemented using interFoam, an open-source solver based on OpenFOAM. InterFoam is specifically designed for multiphase flow problems and employs the VOF method to track the interface between different phases. In this study, pressure–velocity coupling was achieved through the PIMPLE algorithm. Temporal discretization used an implicit Euler scheme, and the time step was dynamically adjusted to maintain the Courant number below 0.5. For spatial discretization, Gauss linear (second-order) schemes were applied to gradient and interpolation terms. The convection term of velocity was discretized using linearUpwind, while the VOF equation employed the vanLeer scheme with interface compression, ensuring bounded and sharp interface capturing.

A representative core sample from the target reservoir was selected, with a diameter of 6 mm and a length of 4 mm. After undergoing oil-washing, salt-removal, and drying treatments, the core was scanned using micro-CT to construct a digital rock geometry model of a weak waterflooding reservoir. Through conventional core physical property testing, the rock sample was measured to have a porosity of 19.9% and a gas permeability of 30.3 × 10−3 μm2.

2.2 Model Parameter Settings

2.2.1 Geometrical Model Parameters

In this study, Avizo software, in combination with the blockMesh and snappyHexMesh utilities in OpenFOAM, was employed to generate the pore-scale mesh. The workflow in Avizo involved a series of image-processing steps, including artifact removal, filtering and noise reduction, threshold segmentation, local contrast-based segmentation, and binarization, which collectively enabled the construction of a three-dimensional digital rock model. The resulting digital rock model was exported in STL format. The surfaceFeatures utility in OpenFOAM was then used to extract the characteristic edges and surface features of the solid skeleton. Subsequently, the computational mesh of the pore structure was generated using the snappyHexMesh method, with surface refinement (level 2) and snapping applied to accurately capture the intricate pore boundaries while preserving the original pore morphology [19,32,35].

A region of 150 × 150 × 350 voxels was extracted as the representative elementary volume (REV), with a single-voxel resolution of 3.0 μm. The geometrical model, shown in Fig. 2, has a porosity of 19.5% and a permeability of 27.6 × 10−3 μm2. These parameters are in good agreement with the results of conventional core property measurements, thus providing a reliable basis for subsequent pore-scale numerical simulations. Additionally, Fig. 3 shows the pore-throat radius frequency distribution for the connected pore network of the 3D digital rock [36,37]. The radii are mainly less than 30 μm but broadly dispersed, indicating pronounced heterogeneity in the pore structure.

Figure 2: Digital rock geometrical model: (a) before and (b) after mesh generation.

Figure 3: Pore-throat radius distribution of the connected pore network.

2.2.2 Mesh Sensitivity Analysis

Previous studies utilizing the VOF method with interface compression schemes have demonstrated the necessity of assessing grid sensitivity [19,31,38,39,40]. Accordingly, to ensure numerical stability and accuracy, a grid independence test was conducted. Three mesh systems with different resolutions were generated based on the same digital rock geometry: Coarse (985,486 cells), Medium (1,224,178 cells), and Fine (1,583,120 cells). The water flooding process was simulated using these three mesh systems under identical boundary conditions.

Fig. 4 compares the recovery factor curves obtained from the three mesh systems. As observed, the results from the coarse mesh show a noticeable deviation, predicting a higher recovery factor (50.8%), likely due to numerical diffusion effects associated with insufficient grid resolution. However, the curves for the medium and fine meshes are in excellent agreement, particularly in the later stages of injection. Table 1 quantifies the results at 1.0 PV. The final recovery factors for the medium and fine meshes are 47.7% and 46.8%, respectively. The relative error between them is less than 2.0%, indicating that the solution has achieved grid independence.

Considering the balance between computational accuracy and efficiency, the medium mesh system (1,224,178 cells) was selected as the baseline for all subsequent simulation schemes in this study.

Figure 4: Comparison of recovery factors simulated using different mesh systems.

2.2.3 Wettability Settings

Based on core samples collected from the target offshore reservoir, a laboratory rock wettability test was conducted to determine the wetting characteristics. Wettability is generally characterized by the contact angle θ, which is defined as the angle formed at the three-phase contact line, measured through the liquid phase between the solid–liquid and liquid–gas interfaces [41]. As shown in Fig. 5, the reservoir exhibited a contact angle of 85.9°.

Figure 5: Experimentally measured contact angle.

3 Model Validation and Scheme Design

3.1 Model Validation

To validate the accuracy of the VOF method in simulating interfacial flow, a laminar two-phase flow model between parallel plates was solved and compared with the analytical solution [42,43]. Fig. 6 illustrates the schematic of the stratified two-phase flow between parallel plates. Two immiscible fluids are distributed between two infinitely long, parallel plates, with the channel extending to a height of 2b in the y-direction, and the channel centerline located at y = 0. Fluid 1 occupies the region 0 < |y| < a, while Fluid 2 occupies the region a < |y| < b. Two constant pressure gradients, G1 and G2, are applied to Fluid 1 and Fluid 2, respectively. Both fluids have the same density ρ, while their viscosities are μ1 and μ2, respectively.

Assuming that the fluid flow between the parallel plates reaches a steady state and satisfies Poiseuille flow, the analytical solution for the velocity profile along the y-direction cross-section, within the positive half-axis between two infinitely long parallel plates, is given by: ux(y)=A1y2+C1,0⩽y⩽aA2y2+B2y+C2,a⩽y⩽b(9) where the coefficients A1, A2, B2, C1, and C2 are, respectively:

Figure 6: Stratified two-phase flow model between parallel plates.

The computational domain measures 10 μm × 200 μm. The initial velocity field was set to zero. Periodic boundary conditions were applied along the x-direction, while no-slip boundary conditions were imposed at the walls. The densities of the two fluids are ρ1 = ρ2 = 1000 kg/m3, with the viscosities of μ1 = 4 mPa·s and μ2 = 50 mPa·s, respectively. The interfacial tension between the two phases is 0.07 N/m. The distance between the parallel plates is 2a + 2b = 200 μm. A constant pressure gradient of G1 = 1 × 106 Pa/m was applied to Fluid 1, and G2 = 0 was applied to Fluid 2. Fig. 7 presents a comparison between the numerical velocity profile across the parallel plates cross section and the corresponding analytical solution. The numerical results agree well with the analytical solution derived from Poiseuille’s law, thereby validating the accuracy of the proposed model and numerical methodology.

Figure 7: Comparison of numerical and analytical solutions for the velocity profile between parallel plates.

3.2 Simulation Scheme Design

Based on the target reservoir core, a three-dimensional digital rock model was constructed (geometrical representation shown in Fig. 1) to conduct pore-scale numerical simulations. The effects of injection velocity, oil–water viscosity ratio, and miscibility degree on microscale flow behavior were systematically analyzed. Specifically, Schemes 1–3 compare cases with different injection velocities, Schemes 2, 4, and 5 examine variations in the viscosity ratio, and Schemes 6–8 focus on different miscibility degrees and injection patterns. The multiphase flow simulation parameters and corresponding scheme distinctions are summarized in Table 2.

4 Results and Discussion

To further indicate the controlling mechanisms of oil–gas–water flow in the target reservoir, multiphase VOF simulations were conducted based on the digital rock model shown in Fig. 1. The study systematically analyzed the effects of injection velocity, oil–water viscosity ratio, miscibility degree, and injection patterns on displacement efficiency. By examining the distribution of residual oil, the influence mechanisms of these factors on reservoir recovery were clarified, thereby identifying the most effective development strategy for the target reservoir. The simulation results not only provide a theoretical basis for reservoir development but also offer technical guidance for subsequent engineering applications.

4.1 Effect of Injection Velocity on Displacement Efficiency

To investigate the effect of injection velocity on the recovery factor and residual oil distribution of the target reservoir, pore-scale two-phase flow simulations were employed based on the model established in the previous section. The injection velocity was set as a variable parameter, with values to 0.0002 m/s, 0.001 m/s, and 0.002 m/s. Water was continuously injected along the positive X-axis, while zero-gradient boundary conditions were applied at both the inlet and outlet. The oil–water distribution results under different injection velocities are shown in Fig. 8, and the corresponding variations in residual oil saturation across slices are shown in Fig. 9a. Variations in injection velocity exerted a significant influence on fingering and channeling phenomena in the water flooding process. According to the calculations and comparison with the log M–log Ca phase diagram summarized by Liu et al. [44], the region corresponding to V = 0.002 m/s falls within the crossover zone, where both viscous fingering and capillary fingering coexist. As V decreases, the displacement behavior becomes increasingly dominated by capillary fingering. As the injection velocity increases, the width of the fingering front decreases, resulting in a non-uniform oil–water distribution and leaving portions of the reservoir unswept and ineffectively displaced. In contrast, lower injection velocities allow water to permeate the reservoir more uniformly. Although residual oil saturation varies with different injection velocity, the overall differences are relatively small. A pronounced reduction in oil saturation was observed near the core inlet, whereas substantial residual oil remained in the central and outlet regions.

Figure 8: Oil–water phase distributions under different injection velocities, with water shown in blue and oil in red.

Figure 9: (a) Residual oil saturation profiles across slices; (b) Variation of recovery factor at different injection velocities.

The variations in recovery factor under different injection velocities are shown in Fig. 9b. The results indicate that the incremental recovery achieved at higher injection velocities is relatively limited. This phenomenon occurs because high injection rates readily induce viscous fingering of the injected water, which promotes subsequent flow preferentially through the previously formed dominant channels. Consequently, oil in adjacent low-permeability zones remains ineffectively displaced, leading to unswept residual oil in certain regions and reduced volumetric sweep efficiency. Although lowering the injection velocity can improve displacement efficiency and enhance oil recovery, field applications often encounter challenges such as reduced liquid production rates and insufficient well productivity, which significantly diminish economic benefits.

4.2 Effect of Oil–Water Viscosity Ratio on Displacement Efficiency

With increasing crude oil viscosity, the viscous resistance and shear stress arising from intermolecular friction become more pronounced, thereby impeding the relative movement between the two fluid phases. To further investigate the influence of oil-water viscosity ratio on recovery efficiency and residual oil distribution in low-permeability reservoirs, pore-scale two-phase flow simulations were conducted based on the model established in the previous section. The crude oil viscosity was set to μo = 3 mPa·s, and different oil-water viscosity ratios (M = 1, 3, 8) were considered, corresponding to injected water viscosities (μ) of 3 mPa·s, 1 mPa·s, and 0.375 mPa·s, respectively, while keeping other conditions unchanged. The oil–water distribution results under different viscosity ratios are shown in Fig. 10, and the corresponding variations in residual oil saturation across slices are shown in Fig. 11a. Based on the calculations and comparison with the log M–log Ca phase diagram, the region corresponding to M = 1 lies within the crossover zone, where both viscous fingering and capillary fingering coexist. As M increases, the displacement process becomes increasingly dominated by capillary fingering. Under the condition of M = 1, the displacement front remains relatively stable, with the water phase advancing uniformly toward the core center, thereby enabling most of the oil to be effectively displaced. Remaining oil decreased most rapidly at the front section of the core and gradually declined toward the middle. As the viscosity ratio increases, channeling phenomena become more pronounced, leading to extensive retention of uncontacted oil at the outlet end.

Figure 10: Oil–water phase distributions under different oil–water viscosity ratio, with water shown in blue and oil in red.

Figure 11: (a) Residual oil saturation profiles across slices; (b) Variation of recovery factor at different oil–water viscosity ratio.

The variations in recovery factor at different oil–water viscosity ratios are shown in Fig. 11b. The results indicate that as the viscosity ratio decreases, the water–oil mobility ratio is reduced, expanding the water swept volume and mitigating viscous fingering. Consequently, the slope of the recovery curve increases progressively, ultimately improving oil recovery. At M = 8, the recovery factor is the lowest, with the growth rate markedly slowing during the mid-stage due to severe instability of the displacement front and early water breakthrough through high-permeability channels. Considering the adverse effects of high oil–water viscosity ratios, mobility control can be achieved by adding polymers or other thickening agents to the injection water. Under optimized low mobility ratio conditions, microscale interfacial dynamics governed by viscous and shear forces markedly enhance the stability and uniformity of the flood front, thereby improving oil mobilization within pore spaces and significantly increasing remaining oil recovery.

4.3 Effect of Miscibility Degree and Injection Patterns on Displacement Efficiency

Based on the pore-scale two-phase flow model established in the previous section, different degrees of miscibility were characterized by varying the interfacial tension and diffusion coefficient. Under the conditions of 120°C and 30 MPa, the density and viscosity of CH4 were calculated as 144.35 kg/m3 and 2.12 × 10−5 Pa·s, respectively. Three CH4 flooding schemes were designed: partially miscible flooding, miscible flooding, and water-alternating-gas miscible flooding (i.e., injecting 0.5 PV of water followed by 0.5 PV of gas). By comparing the effects of different miscibility degree and injection patterns on recovery factor and residual oil distribution in low-permeability reservoirs, the pore-scale flow mechanisms of gas injection processes were revealed.

The oil–water (gas) distribution results under different miscibility degrees and injection patterns are shown in Fig. 12, and the variations in residual oil saturation across slices are shown in Fig. 13a. Under partially miscible flooding, the injected gas (yellow) advances preferentially through high-permeability channels, exhibiting significant fingering, and large portions of low-permeability and dead-end regions remain unswept. Under miscible flooding conditions, where the interfacial tension is reduced to zero, the dissolution and diffusion of CH4 into crude oil are enhanced, resulting in a markedly improved displacement efficiency compared with partially miscible flooding. The lighter yellow coloration indicates a broader gas distribution and a significantly expanded sweep volume. However, due to the strong heterogeneity of the reservoir, residual oil still remained in some low-permeability zones. During the water-alternating-gas flooding process, the water phase initially occupies the inlet and the high-permeability channels. Due to the Jamin effect, gas is subsequently diverted into regions that are otherwise difficult to sweep. This resulted in further mobilization of remaining oil in the front section of the core. The overall residual oil saturation reaches its minimum and exhibits the most uniform distribution, with most slice profiles concentrated in the low-saturation range. This indicates that the synergistic displacement of water and gas significantly enhances both sweep efficiency and microscopic displacement efficiency.

Figure 12: Oil–water (gas) phase distributions under different miscibility degrees and injection patterns, with water shown in blue and oil in red.

Figure 13: (a) Residual oil saturation profiles across slices; (b) Variation of recovery factor at different miscibility degrees and injection patterns.

The recovery factor under partially miscible flooding is the lowest, reaching only 30.9%. In the mid-to-late stages of development, the recovery curve slope decreases rapidly, and the residual oil distribution exhibits high peaks and large fluctuations, indicating an overall non-uniform displacement process with severe gas channeling, leaving substantial residual oil in the reservoir. Compared to water injection, the miscibility between CH4 and crude oil leads to the disappearance of the oil-gas interface, significantly enhancing the swept volume of CH4 and inhibiting CH4 breakthrough. During the transition from water flooding to gas injection, CH4 gradually dissolves in the crude oil, and the “Jia-Min effect” further contributes to improving the microscopic sweep efficiency in the transition process, achieving the optimal displacement effect with a recovery factor of 49.3%. Therefore, for areas with poor water injection development, transitioning from water flooding to gas injection presents significant application prospects.

As shown in Fig. 14, for immiscible displacement, the change in the volume sweep coefficient corresponds to the change in the recovery factor. However, the mass transfer and diffusion effects of CH4 miscibility enhance the CH4 sweep efficiency, and the highest recovery factor is achieved through the synergistic effect of gas and water. However, the water phase restricts the mass transfer miscibility, swelling, and viscosity reduction between CH4 and crude oil, resulting in a reduced CH4 sweep volume compared to miscible displacement.

Figure 14: Volumetric conformance efficiency under 8 different schemes (corresponding to the 8 schemes described in Table 2).

5 Conclusions

This study investigated pore-scale multiphase flow simulations based on the VOF method within the OpenFOAM Finite Volume framework, using plug samples from the target reservoir. The research focused on optimizing water flooding development and identifying the most suitable development methods. The main conclusions are as follows:

• Injection velocity exerts a significant influence on fingering and channeling phenomena during the waterflooding process. Higher injection velocities promote the formation of preferential flow channels, resulting in non-uniform oil–water distributions. Conversely, lower injection velocities allow for more uniform water penetration and improved recovery. Although variations in residual oil saturation occur under different injection velocities, the overall differences remain relatively minor. A pronounced decrease in oil saturation is observed at the core inlet, while substantial residual oil remains near the outlet.
• With increasing oil–water viscosity ratio, water preferentially advances along high-permeability channels, leaving substantial residual oil in unswept regions and reducing displacement efficiency. At M = 1, the water phase advances uniformly, and residual oil is mainly concentrated in low-permeability dead-end regions. In contrast, at M = 8, the recovery factor is the lowest, owing to severe instability of the displacement front and a marked increase in residual oil, which together lead to poor displacement efficiency.
• Partially miscible flooding exhibits the lowest recovery factor, characterized by pronounced gas channeling and high residual oil peaks. Miscible flooding demonstrates better displacement efficiency; however, residual oil still remains in low-permeability regions. The recovery factor during the water-alternating-gas (WAG) flooding process reached 0.49, representing absolute increases of 18.4% and 8.0% compared with the partially miscible and miscible flooding cases (0.31 and 0.41), respectively. This process improves sweep uniformity by blocking high-permeability channels with the water phase and enabling the gas phase to displace unswept regions, thereby significantly enhancing oil recovery and microscopic displacement efficiency.

Current research primarily focuses on multiphase flow simulations at the pore scale. In the future, it should integrate both micro and macro scales to develop more accurate multiscale models. By combining key parameters at the micro scale with those at the macro scale, upscaled mathematical models can be constructed to more accurately describe the flow characteristics of reservoirs, thereby improving the prediction accuracy of oil recovery efficiency and the effectiveness of different development strategies.