GAGD, as an exceptional gas injection and extraction process, effectively integrates horizontal wells with stabilized interfacial flow to achieve enhanced recovery. According to the existing literature, the majority of gas-assisted gravity flooding projects have been recorded in the United States, Canada, and the Middle East. Among these, notable successful cases include the West Hackberry Oilfield and Hawkins Oilfield in the United States. The West Hackberry salt dome reservoir, located in the United States, exhibits high dip angle and high permeability, making it susceptible to fingered flow under traditional water flooding techniques. The GAGD process injects gas to establish an artificial gas cap at the top, enabling a dual-drive mechanism that combines gas drive and gravity drive (Abdelaal et al. 2023; Rao et al. 2004; Kasiri and Bashiri 2009; Al-Mudhafar 2016). The field test shows that the recovery ratio is close to 90%, which is obviously improved compared with the traditional water drive recovery ratio of 50%–60% (Backmeyer et al. 1984). Hawkins oilfield is a fault-block reservoir with a dip angle of 6 and sufficient water energy. After gas-assisted gravity flooding, the oil recovery efficiency is increased by 20% compared with water flooding (King and Lee 1976). The core experiment also proves that gas-assisted gravity flooding can still displace the remaining oil in the water invasion area and greatly reduce the residual oil saturation, which is also the first gas-assisted gravity flooding project applied to fault-block reservoirs (Carlson 1988). Most large carbonate reservoirs in the Middle East employ top gravity- stabilized flooding as the initial stage, achieving a superior recovery ratio of over 50%, which is significantly superior to that of water-flooding development reservoirs (Ma et al. 2015). The field application of gas-assisted gravity flooding has achieved remarkable success (Kulkarni and Rao 2006).

The recovery ratio of GAGD technology depends on several factors, including rock wettability, gas injection velocity, and rock heterogeneity (Al-Mudhafar 2018; Khorshidian et al. 2018; Akhlaghi et al., 2012; Kong et al. 2020; Xu et al. 2024; Wang et al. 2024; Yang et al. 2023). Consequently, when evaluating the performance of GAGD technology in tight glutenite reservoirs, it is essential to consider the impact of structural heterogeneity and the intricate wetting behavior of glutenite. Misagh’s numerical simulations reveal that, in homogeneous reservoirs, the maximum gas-assisted gravity flooding recovery efficiency stands at 40%. However, in heterogeneous reservoirs, this recovery efficiency is contingent upon the fracture density. Al-Mudhafar reinforced this notion through comprehensive core displacement experiments, thus highlighting the potential of gas-assisted gravity flooding in heterogeneous reservoirs (Al-Mudhafar and Rao 2017; Al-Mudhafar 2018). The wettability of glutenite reservoirs is a highly complex issue. The presence of various gravel types and distinct mineral compositions within the rock pore spaces leads to spatial variations in rock wettability, giving rise to heterogeneous wettability. This phenomenon would lead to the complexity of the interplay between capillary force and wettability (Wang et al. 2023). Khorshidian changed the wettability of the model with the help of microscopic etching visualization experimental technology. The results show that the recovery degree of oil-wet reservoir is better than that of water-wet reservoir (Khorshidian et al. 2018). Through laboratory experiments and grid model research, Morrow found that the recovery efficiency of mixed wet reservoir is higher than that of water wet reservoir (Mason and Morrow 1991).

The dimensionless method has become widely used in studying the GAGD technique due to the multitude of influencing factors involved. Kulkarni and Rao examined the impact of important dimensionless groups on final recovery in various mixed-phase and unmixed-phase gas gravity drainage experiments (Kulkarni and Rao 2006). For free-fall gravity drainage, Grattoni et al. proposed a new dimensionless group that combines the effects of gravity, viscosity, and capillary forces (Grattoni et al. 2001). Wood et al. introduced 10 dimensionless groups to characterize CO₂ drive in an inclined waterflood reservoir and utilized them for experimental design to establish screening criteria applicable to Gulf Coast reservoirs (Wood et al. 2006). Rostami et al. investigated the impact of Bond number and capillary number on the efficiency of forced (gas-assisted) gravity drainage (Rostami et al. 2010).

The gas displacement process within glutenite is a complex flow behavior that is fundamentally determined by the microscale characteristics of the glutenite (such as pore structure heterogeneity and wettability of flow channel walls) as well as external factors (injection velocity, fluid properties, stress field distribution). Previous research on the GAGD process has mainly focused on laboratory experiments. However, due to the unique strong heterogeneity and complex wettability of glutenite, it is challenging to obtain direct and dynamic visualization of the flow process through physical experiments, which are also characterized by high randomness and difficulty. Therefore, it is necessary to build upon previous research and explore the microscale flow characteristics within tight glutenite pores using numerical simulation methods. There are two types of numerical methods, namely direct (Xie et al. 2023a, b; Xie et al. 2024) and indirect (Qin et al. 2016) simulation approaches, to solve the Navier–Stokes equations. This paper employs the finite element method as a direct numerical tool, coupling the Navier–Stokes equations with the phase field method, to monitor two-phase flow interfaces and investigate two-phase flow phenomena of the GAGD technique in tight glutenite at the pore scale. A quantitative analysis is conducted to assess the structural and wettability heterogeneities of glutenite, as well as to identify the key factors influencing the GAGD process in heterogeneous glutenite. Additionally, by incorporating dimensionless numbers such as capillary number, Bond number, gravity number, gas–water viscosity efficiency, and heterogeneity coefficient, we can comprehensively evaluate the impact of these dimensionless numbers on the flow characteristics of two-phase flow and the effectiveness of oil displacement.

2.1 Physical problem

The Mahu oilfield in Xinjiang, China, is the largest glutenite oilfield in the world, with estimated reserves exceeding one billion tons. It has attracted significant attention globally. (Xia et al. 2021). This reservoir belongs to the category of tight glutenite reservoirs, posing unique challenges compared to conventional conglomerate reservoirs (Fig. 1). Glutenite oil reservoirs exhibit strong heterogeneity, complex wettability distribution, which significantly complicate exploitation efforts (He et al. 2021; Xie et al. 2023a, b; Pu et al. 2022).

Location of Mahu glutenite reservoir

The Ma Lake tight glutenite reservoir is characterized by low permeability. Compared to conventional sandstone reservoirs, the differences in gravel size result in a more complex micro-pore structure, exhibiting pronounced structural heterogeneity. Additionally, the glutenite is composed of gravel with varying material properties, which further contributes to its inherent non-homogeneous wettability. As revealed in Fig. 2, the XRF scanning results of the Mahu tight glutenite rock sample demonstrate that the gravels are predominantly composed of Al, Si, K, and Fe elements. Different mineral properties and content can induce various heterogeneities (Yang et al. 2024). Clay minerals often exhibiting hydrophilicity. Alinejad found through comparisons that samples with higher quartz content than clay minerals tend to develop oil-wettability (Alinejad and Dehghanpur 2021). Zhu also demonstrated the development of multiple oil-wet states with higher quartz content (Zhu et al. 2020).

Microscopic structure of glutenite. a XRF scanning results; b Schematic diagram of mineral wettability; c Element results

2.2 Theoretical model

The Navier–Stokes equations are used to govern the two-phase flow of oil and gas, including the conservation of mass and momentum. As shown in Eqs. (1) and (2).

where ρ is density, kg/m³; μ is dynamic viscosity, Pa/s; u is viscosity, m/s, p is pressure, Pa; g is gravity acceleration, m/s²; F_st is the surface tension acting on the oil–water interface, Pa/m, calculated by Eq. (3).

The phase field method is used to track the two-phase interface. The phase field method is one of the commonly-used methods in multiphase flow, which is mainly based on the change of phase field variables to respond to the change process of the whole system in a relatively simple form. The multi-field coupling effect is generally added to the energy equation with relevant variable parameters. The energy equation mainly uses the Cahn–Hilliard equation, which can represent the evolution of the two-phase diffusion interface. For numerical simulation of fluid flow, the momentum conservation equation is generally used to describe the fluid flow, taking the wetting angle and interfacial tension into account. In the process of water driven oil, the dimensionless variable φ is transferred from − 1 to 1. The phase field equation can be decomposed into the following two equations:

where φ is the dimensionless phase field variable; u is the fluid velocity, m/s; γ is the mobility, m³ s/kg; λ is the mixing energy density; ε is the interface thickness parameter (a capillary width), m; Ψ is the dimensionless auxiliary variable (related to chemical potential).

The interface thickness parameter ε usually can be equal to the half of the maximum mesh element size in the region through which the interface passes. The mobility γ governs the diffusion-related time scale for the interface, and is determined by a mobility tuning parameter χ (m s/kg) and the interface thickness ε, γ = χε²; and χ usually can be equal to 1 m s/kg for most cases. The mixing energy density λ is related to the surface tension coefficient and interface thickness, as shown Eq. (6).

where σ is the interfacial tension coefficient, N/m.

Equations (6)–(11) illustrate that the key interconnected variables between the Navier–Stokes equations and the phase-field equations are density, dynamic viscosity, and surface tension.

where V_f1 and V_f1 mean the volume fractions of fluid 1 and 2, respectively; ρ_f1 and ρ_f1 mean the density of fluid 1 and 2, kg/m³, respectively; μ_f1 and μ_f1 mean the dynamic viscosity, of fluid 1 and 2, Pa/s, respectively.

2.3 Numerical setting

This paper uses the Computational Fluid Dynamics (CFD) module in simulation software. The mesh node and PARDISO solver are used to solve the model. As shown in Fig. 3, the computation domain has a size of 39.6 mm × 39.6 mm and contains randomly distributed spherical particles of different diameters (1.725 mm and 0.9 mm) to simulate the tight glutenite. The pore throat widths range from 0.275 to 1.1 mm. Additionally, a buffer zone with a width of 0.5 mm is set on the top of the model to ensure a steady inject velocity distribution. The contact angle and interfacial tension are 45° and of 20 mN/m, respectively. Other parameters of oil and water used in the simulation can be found in Table 1.

Computation domain and mesh

The domain is discretized by the mesh that has sufficient accuracy while maintaining the optimal number of elements, containing a total of 107,366 elements. The maximum element size is 2.12 mm; the minimum element size is 0.012 mm; the maximum element growth rate and the curvature factor are 1.3 and 0.3, respectively.

The initial pressure is cause by the gravity. The initial velocity is zero. The injection velocity is set to 10 mm/s at the upper boundary. The lower boundary is set as the fixed pressure equal to the value caused by the gravity. Two side boundaries are no-flux condition. During the solving process, the time step of the transient solver is set to be 0.02 s.

3.1 The role of gravity assistance in gas drainage

3.1.1 Effect of gravity on flow front

The objective of this section is to investigate the impact of gravity on the stabilization of the drive front (leading edge of the moving gas). Figure 4 presents the numerical simulation process of the GAGD process that takes into account the influence of gravity. The injected gas forms a gas roof in the upper part due to the density difference, and the drive front of the gas roof remains quite stable in the whole process. The finger-flow phenomenon is not significant, and the gas-drive front moves steadily downward under the influence of gas pressure, except for the residual oil trapped in some narrow pores as a result of capillary forces.

Temporal and spatial evolution of gas and oil volume fraction with the consideration of gravity. a t = 0.3 s; b t = 0.6 s; c t = 0.9 s; d t = 1.2 s; e t = 1.5 s; f t = 1.8 s; g t = 2.1 s; h t = 2.4 s (gas breakthrough); i t = 4 s (steady spatial distribution of fluids)

3.1.2 Effect of gravity on recovery efficiency

This part is intended to analyze the effect of gravity on enhanced recovery efficiency. As shown in Fig. 5. When gravity is incorporated into the gas drive process, the recovery efficiency increases linearly with the increase of CO₂ injection. The gas-drive front moves steadily downward with the aid of gravity, leading to higher recovery efficiency. The final recovery efficiency is not reached at the time of BT (gas breaks through the domain), but still increases slightly with subsequent CO₂ injection. Another intriguing finding is that, once the final recovery efficiency has stabilized, residual oil trapped in tiny pores at higher elevations will undergo flow and diffusion towards lower elevation pores. This occurrence may be attributed to the combined action of gravity and CO₂ pressure, which reduce the residual oil, overcoming the constraints imposed by viscosity and capillary forces.

Recovery efficiency evolution pattern during GAGD oil driving process

3.1.3 Jumping phenomenon of local pressure during GAGD

The GAGD process exhibits unstable pressure fluctuations during CO₂ transportation and filling in pore throats of different sizes (Armstrong et al. 2015). These fluctuations mostly occur during the transition of gas flow from small to large pores, as illustrated in Fig. 6. Once the pressure exceeds the capillary entry pressure of the pore throat (Berg et al. 2013), the pore body volume is rapidly filled by the incoming fluid (Haines et al. 1930). Figure 6a demonstrates the pressure jump phenomenon during 0.86–1.0 s, with the time step of 0.02 s as a transport trajectory. When the gas-drive front reaches the pore entrance, the pressure at the monitoring point shows a certain fluctuation and a slow decreasing trend. As the gas-drive front enters the pore space, the pressure exhibits a steep increase until the pore space is fully occupied, followed by a period of fluctuation and gradual decrease (Fig. 6a, b). Due to the change of pressure caused by the jump, the backward shift phenomenon of gas-drive front occurs at the smaller pore throats, as shown in Fig. 6c.

Local pressure jump phenomenon. a CO₂ transport trajectory in 0.86–1.0 s; b Monitoring point at the large pore space during the Haynes’ jump pressure change; c CO₂ transport trajectory at the right small pore space

3.2 Sensitivity analysis of factors that influence GAGD

3.2.1 Effect of heterogeneity on recovery efficiency

1. (1)

   Structure heterogeneity

The complex diagenetic process of glutenite leads to the spatial variation of its internal structure. We use the dimensionless variation coefficient of permeability to quantitatively characterize the structure heterogeneity, as shown in Eq. (12). The dimensionless variation coefficient of permeability, \(V_{k}\), represents the degree of dispersion or variation of the permeability within a layer relative to its average value. Generally, the dimensionless variation coefficient of permeability is greater than 0, and the larger value means the stronger heterogeneity.

where \(V_{k}\) is the dimensionless variation coefficient of permeability in the layer; \(k_{i}\) (unit: m²) is the permeability of the i^th sample in the layer; \(\overline{k}\)(unit: m²) is the average value of permeability of all samples in the layer; and n is the number of samples in the layer.

The structure heterogeneity is achieved by randomly increasing or decreasing the amount of large gravels. Figures 7a, b show the high and low heterogeneity models, respectively. In Fig. 7c, the domain is divided into four equal samples to calculate \(k_{i}\), and consequently the dimensionless variation coefficient of permeability, by the Darcy’s Law as shown in Eq. (13).

where \(u\) is the seepage velocity, m/s; \(K\) is the permeability, m²; \(\mu\) is the dynamic viscosity of the fluid, Pa s; and \(\nabla P\) is the pressure gradient, Pa/m.

Heterogeneous domain. a High structure heterogeneity; b Low structure heterogeneity; c Four samples used to calculate the dimensionless variation coefficient of permeability of the high heterogeneous domain

Generally, when \(V_{k}\) < 0.5, the degree of heterogeneity is weak; when \(V_{k}\) is between 0.5 and 0.7, the degree of heterogeneity is of medium level; and when \(V_{k}\) > 0.7, the degree of heterogeneity is strong. According to this criterion, we determine the three cases with weak, medium, and strong heterogeneity, and the \(V_{k}\) are equal to 0.49, 0.67, and 0.88, respectively.

Theoretically, lower inject velocity can achieve higher recovery efficiency. However, in the real engineering, the injection velocity cannot be reduced to the theoretical minimum due to the consideration of production cycle and cost return. Therefore, for the characteristics of strongly heterogeneous glutenite, we simulate two scenarios with low and high injection velocity. As shown in Fig. 8, at a low injection velocity, the gas drive front moves steadily downward, and the injected gas sweeps downward the residual oil in the pores. At this situation, gravity and capillary force dominate the flow process, and only a small amount of oil in the narrow pores is difficult to be driven out. Therefore, the GAGD process is relatively less affected by the heterogeneity at lower injection velocity.

Influence of structure heterogeneity on fluid distributions at lower injection velocity. a At BT in case of \(V_{k}\) = 0.49; b At BT in case of \(V_{k}\) = 0.67; c At BT in case of \(V_{k}\) = 0.88; d Final in case of \(V_{k}\) = 0.49; e Final in case of \(V_{k}\) = 0.67; and f Final in case of \(V_{k}\) = 0.88

As shown in Fig. 9, when the injection velocity is high, the stabilizing effect of gravity on the fluid level becomes worse, and the viscous force dominates at this situation. When the heterogeneity is weak, the injection of CO₂ can facilitate the formation of bypass oil at the narrow pore throats, which is caused by the large capillary force at the narrow pore throats. In the strong heterogeneity model (as shown in Figs. 9a, c), the pore throat distribution is very complex, and the complex pore throat distribution makes the injected CO₂ mainly move along the dominant flow channel. Strong finger-flow phenomenon occurs during the drive process, forming a large area of bypass oil on both sides. When some fine parallel dominant flow paths converge, they will wrap a piece of stagnant oil zone, which will lead to a further reduction in recovery efficiency. Therefore, the GAGD process is strongly influenced by the heterogeneity at high injection velocity.

Influence of structure heterogeneity on fluid distributions at high injection velocity. a At BT in case of \(V_{k}\) = 0.49; b At BT in case of \(V_{k}\) = 0.67; c At BT in case of \(V_{k}\) = 0.88; d Final in case of \(V_{k}\) = 0.49; e Final in case of \(V_{k}\) = 0.67; and f Final in case of \(V_{k}\) = 0.88

As shown in Fig. 10, when the injection velocity is low, different heterogeneous models all have higher oil recovery efficiency. With the increase of injection velocity, the recovery efficiency of each model has entered a rapid decline stage. Furthermore, its decline rate is intricately linked to its heterogeneity coefficient. The stronger the heterogeneity, the greater the decline rate. After the injection velocity exceeds 8 mm/s, the recovery efficiency decreases less with the increase of injection velocity in spite of some fluctuations in the case of strong heterogeneity.

1. (2)

   Wettability heterogeneity

Relationship between injection velocity and recovery efficiency under different structural heterogeneity

The nature of gravel types and mineral compositions result in the heterogeneity of wettability of glutenite. Oil and non-oil fluid has a competition in terms of wetting (holding contact with) the gravel surface. In the oil–gas–rock system, the smaller value of the contact angle means the greater oil-wetted area on the solid surface, so-called oleophilic. The petroleum industry standard specifies the oil–solid surface contact angle in the range of [75°, 105°) as neutral wetting and in the range of [105°, 180°) as oil wetting (Cao et al. 2021), which can be described by the oil–solid contact angle as shown in Fig. 11.

Schematic distribution of gravels of different wettability

The wettability heterogeneity index \(I_{\text{w}}\) is used to quantitatively describe the wettability heterogeneity, defined as the volume fraction of gas-wetted surfaces over the total pore surface, as shown in Eq. (14). The range of \(I_{\text{w}}\) varies from 1, indicating gas-wetted conditions, to 0 for neutral wettability, and to − 1 for completely oil-wetted conditions. (Su et al. 2018).

In this study, we randomly adjusted the contact angles of some of the wetted walls in the model to 135° (the blue edge in Fig. 12), while keeping the remaining walls at 45° unchanged, in order to simulate the wettability heterogeneity observed within glutenite.

Schematic distribution of gravels of different wettability

Influence of wettability heterogeneity on fluid distributions is shown in Fig. 13, with wettability heterogeneity index equal to 0.5, 0, and −0.5, respectively. When the pore medium at the smaller pore size is not oil-wetted after CO₂ injection, due to the role of interfacial tension and its large capillary force, it will make part of the oil remain in the pore space and cannot be driven out, resulting in lower recovery. Figure 14 presents a comparison of the displacement efficiency between gas-wetted and oil–gas mixed-wetted gravel wetting walls. When the pore medium is oil-wetted at the smaller pore size, the oil–gas interface will be concave to the oil phase during the displacement process. The formation of this concave to curved interface is caused by the pressure difference between the oil and gas two-phase non-mixed-phase fluid in the capillary. As the oil wets the capillary wall, the pressure in the non-wetted phase is greater than that in the wetted phase, and the result is that the oil–gas interface is concave to the oil phase (Figs. 14f–h). As the wetted wall is oil-wet at this time, the residual oil interfacial tension is small, and the capillary force to be overcome is relatively small. During the displacement process, oil exhibits continuous phase migration with enhanced fluidity, enabling it to overcome capillary forces and be discharged from narrow pore spaces, ultimately enhancing overall oil recovery efficiency.

Influence of wettability heterogeneity on fluid distributions. a At BT in case of \(I_{\text{w}}\) = 0.5; b At BT in case of \(I_{\text{w}}\) = 0; c At BT in case of \(I_{\text{w}}\) = − 0.5; d Final in case of \(I_{\text{w}}\) = 0.5; e Final in case of \(I_{\text{w}}\) = 0; and f Final in case of \(I_{\text{w}}\) = − 0.5

Fluid distribution near the oil-wetted wall during the drive-off process. a Oil-wetted diagram; b Early stage; c Mid - stage; d Late stage; e Gas-wetted diagram; f Early stage; g Mid - stage; h Late stage

Figure 15 shows the distribution of BT time and final recovery efficiency of GAGD process under different wettability index. According to the previous analysis, when the wettability index is greater than zero, the overall wettability of the model is gas-wetted, and the recovery efficiency is low. As the wetting index gradually decreases and the proportion of oil-wet walls increases, the overall wettability of the model gradually shifts towards oil-wet, resulting in an increase in recovery efficiency.

Distribution of BT time and final recovery at different wettability indexes

3.2.2 Effect of the capillary number and oil–gas viscosity ratio

The Capillary number (Ca), also known as the interfacial tension number, is a dimensionless quantity initially proposed by Taylor in 1934. The Capillary number is defined as the ratio of fluid viscosity to interfacial tension, and it reflects the balance between different forces during the displacement process of two phases in a porous medium. The expression is given by

where the Ca is the Capillary number, μ represents the viscosity of the injected CO₂ in Pa·s; v is the injection velocity of CO₂ in m/s; σ is the interfacial tension between the oil and gas phases in N/s.

The Viscosity ratio (M), defined as the ratio of oil viscosity to gas viscosity, is used to indicate the difference in viscosity between the two phases. It is a dimensionless quantity. The expression is given by

where \(\mu_{\text{O}}\) is the viscosity of the oil phase in Pa s; \(\mu_{\text{C}}\) is the viscosity of the gas phase in Pa s.

Figure 16 shows the fluid distribution during breakthrough for different Ca and M values, with Log Ca values of - 4.5, - 4.7, - 5.7, and Log M values of 1.7, 2.7, 3.0. The flow process is primarily governed by three forces: gravity, capillary force, and viscous force. When Log Ca ≤ - 5.7, with the assistance of gravity, gas drives the interface to stabilize and move downward, resulting in a high displacement efficiency. Only residual oil is present in some narrow pores, where the capillary force plays a dominant role. As Log Ca increases, the assistance of gravity weakens, and the capillary force continues to dominate. However, some injected gas bypasses certain narrow pores, leading to the formation of residual oil. The flow pattern deviates from the viscous fingering pattern. When Log Ca ≥ - 4.5, there is a significant fingering phenomenon during the displacement process, with gas rapidly breakthrough and forming large patches of residual oil before breakthrough. In this case, the viscous force plays a dominant role, and the flow pattern tends toward viscous flow. The leading interface of gas displacement exhibits finger-shaped advancement, and the continuous shape of large patches of residual oil is related to the flow characteristics of the two fluids.

Fluid distribution during displacement breakthrough under different Log M and Log Ca conditions

Increasing the oil–gas viscosity ratio leads to an increase in the mobility ratio between gas and oil, resulting in the gas exhibiting higher flowability than the crude oil. This reduces the stability of the gas displacement front and makes it susceptible to bypassing oil. In contrast, under low viscosity ratio conditions, the gas displacement front has higher stability, and the front interface moves downward, effectively sweeping the residual oil in the pore space with a displacement efficiency approaching 100%. Additionally, even after BT, higher final recovery efficiency is observed under low viscosity ratio conditions. This indicates that gas displacement under low viscosity ratio conditions can better expel the crude oil from the reservoir and improve the overall recovery efficiency. Figure 17b illustrates the relationship between recovery efficiency at BT and the final recovery efficiency under different viscosity ratio conditions. This further confirms the superior recovery performance under low viscosity ratio conditions.

Recovery efficiency under different capillary number and oil–gas viscosity ratio at BT. a Recovery efficiency under changing Ca; b Recovery efficiency under changing M

Based on the numerical simulations of the GAGD process discussed earlier, three displacement modes are mapped onto the Log Ca-Log M stability phase diagram. Due to the difficulty in defining precise flow regime boundaries in small-scale porous media, the light blue region in Fig. 18 is used as a transitional zone to account for uncertainties and prevent errors caused by them.. The numerical simulation results in this study mostly fall within the range of Log M < - 1.5 and Log Ca < - 4.5, which belong to the regions dominated by viscous forces or capillary forces. This consistency with the dominant force analysis in the displacement process supports the findings in the previous discussion.

Log Ca-Log M stability phase diagram, the light blue area is the boundary area

3.2.3 Effect of the Bond number on recovery efficiency

The Bond number (N_B) is a dimensionless quantity determined by the influence of surface tension, which characterizes the ratio of liquid gravity to surface tension. It is expressed as:

where ρ is the density of the liquid in kg/m³, g is the acceleration due to gravity in m/s², L is a characteristic length (typically the diameter of the radius of a capillary) in m, and σ is the surface tension of the liquid–gas or liquid–solid interface in N/m.

Figure 19 shows the fluid distribution at different Bond numbers (N_B) from left to right, with N_B values of 3.23×10⁻⁴, 5.14×10⁻⁴, and 7.48×10⁻⁴, respectively. When N_B is relatively low, the gravity-assisted drainage effect is not significant, and gravity plays a minor role in the displacement process, while capillary and viscous forces are the main driving forces. In this case, the injected CO₂ will advance along preferential flow paths, causing fingering phenomenon and resulting in relatively low recovery efficiency with continuous bypassed oil on either side. As N_B increases, the formation of bypassed oil during the displacement process decreases and is limited to smaller areas with smaller pores. The gravity-assisted drainage strengthens.

Influence of Bond number on fluid distributions. a At BT in case of \(N_{\text{B}}\) = 3.23×10⁻⁴; b At BT in case of \(N_{\text{B}}\) = 5.14×10⁻⁴; c At BT in case of \(N_{\text{B}}\) = 7.48×10⁻⁴; d Final in case of \(N_{\text{B}}\) = 3.23×10⁻⁴; e Final in case of \(N_{\text{B}}\) = 5.14×10⁻⁴; and f Final in case of \(N_{\text{B}}\) = 7.48×10⁻⁴

When N_B = 7.48×10⁻⁴, the displacement process is characterized by stable gravity-assisted drainage. Under the influence of gravity, the front interface of gas displacement moves steadily, sweeping the oil in the pores, and the spreading efficiency of CO₂ approaches 100%. The gravity-assisted drainage effect is significant, and it works together with capillary and viscous forces as dominant factors in the displacement process. Figure 20 displays the recovery efficiency and ultimate recovery efficiency at different pore sizes at BT.

Recovery efficiency and ultimate recovery efficiency under different pore sizes at BT

3.3 Relationship between dimensionless number and GAGD recovery

As shown in Fig. 21a, the relationship between the capillary number N_C and recovery factor under different Bond numbers is depicted. Based on the image, it can be concluded that there is no definitive relationship observed between the capillary number N_C and the recovery factor. Therefore, it is not possible to predict the recovery factor solely based on the capillary number N_C. However, when the Bond number is held constant, the recovery efficiency decreases with an increase in the capillary number.

Correlation of dimensionless numbers a The relationship between Capillary number N_C and oil recovery under different Bond numbers. b The relationship between Bond number and oil recovery

According to Fig. 21b, the fitted curve demonstrates a strong correlation between the Bond number and the recovery factor. When the Bond number remains constant, an increase in the recovery factor leads to a decrease in the effectiveness of capillary forces, while viscous forces gradually take over the dominant role in the flow process (as indicated by the blue box in Fig. 21b). As the Bond number increases, both capillary forces and viscous forces diminish, and gravity begins to play a more significant role in the flow process (as indicated by the green box in Fig. 21b). With the increased influence of gravity, the stability of the gas displacement front interface improves, reducing fingering and gas channeling phenomena. This expansion of the gas injection’s spreading range enhances the displacement efficiency and subsequently increases the recovery efficiency.

To find a better correlation, Rostami et al. (2010) proposed a combined dimensionless number that is directly proportional to the Bond number and viscosity ratio, while inversely proportional to the capillary number. Building upon their concept and considering heterogeneity factors, a new dimensionless number, N_Glu, is introduced for predicting the recovery factor in GAGD processes of heterogeneous glutenite reservoir. N_Glu is defined as follows:

where α, β and γ are correlation coefficients.

The dimensionless number N_Glu is obtained through non-linear regression analysis with three coefficients: 1.8, 0.13, and − 22. The formula for N_Glu is expressed as follows:

According to Fig. 22, it can be observed that the new dimensionless parameter N_Glu exhibits a strong correlation with the recovery efficiency, with a high coefficient of 0.93965. This correlation includes the heterogeneity coefficient V_K and can be utilized to predict the recovery efficiency during the exploitation of highly heterogeneous glutenite reservoir using the GAGD technique.

Relationship between N_Glu and oil recovery

As shown in Fig. 23, there is a strong correlation between the dimensionless number N_Glu and the recovery efficiency, with a correlation coefficient of 0.93965, which is significantly higher than other models (such as Fig. 23b). In addition to capillary number, viscosity ratio, and Bond number, the dimensionless number for pore scale also takes into account the strong heterogeneity of glutenite reservoirs. It can be used to predict the recovery efficiency of glutenite reservoir using the GAGD technique.

a Relationship between recovery efficiency and dimensionless number; b Correlation comparison

This study aims to reveal the role of gravity in GAGD and the impacts of heterogeneity on the oil recovery process by microscopic scale immiscible two-phase flow modelling. The mechanism of gravity in gas-oil drive is revealed in terms of stabilizing the drive front and increasing the recovery efficiency. The influences of structure heterogeneity and the wettability heterogeneity are analyzed. Main findings are concluded as follows.

1. (1)

   Gravity plays a role of forming a gas roof on the top, so the gas-drive front can move steadily downward under the push of the gas top, except for the residual oil left in some narrow pores due to capillary force.

2. (2)

   When CO₂ is transported and filled in different sized pore throats, unstable pressure fluctuations are observed, and consequently the backward shift phenomenon of gas-drive front occurs at the smaller pore throats.

3. (3)

   Structure heterogeneity generally has a negative influence on the oil recovery efficiency, though low injection velocity could decline such impact.

4. (4)

   Wettability heterogeneity has a significant influence on GAGD process. With the decrease of wettability index and the proportion of oil-wet walls increases, more oil can successfully overcome the capillary force to be recovered.

5. (5)

   The recovery efficiency is directly proportional to the Bond number and inversely proportional to the capillary number and viscosity ratio. A new dimensionless number, N_Glu, is proposed, incorporating the gravity number, Bond number, viscosity ratio, and heterogeneity coefficient. This aims to enhance the accuracy of predicting the recovery efficiency in highly heterogeneous glutenite reservoir using the GAGD technique.