Microfluidic investigation of CO₂ foam flow in a heterogeneous porous medium

Abstract

CO₂ foam has emerged as a promising alternative to CO₂ gas for mobility control in enhanced oil recovery (EOR) applications, yet the pore-scale dynamics and mechanisms governing foam flow and oil displacement in heterogeneous porous media remain relatively underexplored. This study investigates the pore-scale behavior of CO₂ foam flow and its efficacy for EOR in a microfluidic heterogeneous porous medium under ambient conditions. The microfluidic device, composed of parallel high- and low-permeable regions with a permeability ratio of 5.8, enabled direct visualization of foam generation, propagation, and oil displacement dynamics. A flow-focusing geometry was used to produce stable foam bubbles ranging from 30–270 μm, with size and texture governed by gas pressure and liquid flow rates. The foam morphology and transport characteristics were further analyzed as functions of the gas injection ratio (R_g), revealing that increasing R_g led to a higher gas areal fraction (F_g), increased gas trapping, and reduced foam velocity (V̄_f) and texture (n_f). Oil displacement experiments using SDS solution, CO₂ gas, and CO₂ foam showed distinct differences in performance. SDS solution and CO₂ foam formed relatively stable displacement fronts in the high-permeable zone, while CO₂ gas exhibited severe viscous fingering and early breakthrough. Foam flooding achieved significantly higher oil recovery (up to 100%) at faster rates compared to gas injection, owing to foam-induced pore blocking and vertical flow diversion into the low-permeable zone. These findings provide insights for optimizing foam-based EOR processes and highlight the value of microfluidics for resolving multiphase transport phenomena at the pore scale.

---

1 Introduction

Foams comprise gas bubbles dispersed in a surfactant solution,^1,2 forming a dynamic, multiphase system often described as a colloidal dispersion.³ However, foams differ from classical colloidal systems due to their athermal nature, with structural dynamics governed primarily by interfacial and hydrodynamic forces,⁴ enabling them to modify flow behavior in porous media. Foams have been widely applied in subsurface processes, including enhanced oil recovery (EOR),^5–7 soil and aquifer remediation,⁸ carbon dioxide sequestration,^9,10 hydraulic fracturing,¹¹ and hydrogen storage.^12,13 Foam improves displacement efficiency by mitigating gas channeling and enhancing sweep efficiency, compared to gas, water, or surfactant flooding alone, particularly in heterogeneous media.¹⁴

CO₂ gas injection is a well-established EOR technique,¹⁵ expected to be used for 1.6 million barrel of oil production per day by 2040.¹⁶ However, its poor mobility control due to the gas's inherently low density and viscosity drives viscous fingering, gravity override, and preferential flow along high-permeable zones.¹⁷ This unfavorable behavior is quantified by the mobility ratio (M_R) between the displacing and residing fluids, defined as:

with k_r and μ representing the relative permeability and viscosity of the displacing (D) and residing (R) fluids. M_R > 1 denotes conditions prone to instability and poor sweep efficiency. Such mobility contrast limits recovery (up to 30%¹⁸) in light oil systems and becomes more severe in heavy oil (high μ_R) reservoirs as M_R ≫ 1.¹⁹

Foam mitigates these mobility challenges by increasing the displacing fluid's viscosity (μ_D), thereby shifting M_R toward more favorable values and suppressing unstable flow.²⁰ The enhanced resistance in foam flow stems from interfacial viscosity at gas–liquid boundaries²¹ and capillary resistance from foam lamellae confined in pores.²² Foam can also accumulate in high-permeable regions, diverting flow into less-swept, low-permeable zones.⁶

Conventional core flooding has long been used to assess foam performance in porous media, wherein fluids are injected through rock cores to quantify oil displacement under field-relevant conditions.^23,24 While valuable for evaluating bulk displacement efficiency, these methods are inherently limited by the opacity of geological samples, which obscures the pore-scale mechanisms governing foam evolution and oil displacement.^25,26 To address this, 2D Hele–Shaw cells of different configurations, such as empty chambers,^20,27 bead packs,²⁸ and patterned plates,²⁹ have been used to offer visual access but are limited in mimicking realistic pore geometries. Advanced imaging tools, such as X-ray computed tomography (CT)³⁰ and nuclear magnetic resonance (NMR),³¹ have also been applied for detailed internal observations, but are constrained by high cost, resolution limits, and experimental complexity.³²

Microfluidic platforms have emerged as powerful tools for resolving pore-scale processes with exceptional spatial and temporal resolution.^{25,26,33–36} These lab-on-a-chip systems replicate complex porous architectures on transparent chips, enabling direct observation of foam generation, transport, and collapse under controlled flow conditions. Advantages such as optical clarity, low fluid consumption, high reproducibility, and rapid prototyping make microfluidics an increasingly indispensable asset for studying multiphase interactions, wettability alterations, and interfacial transport in porous media.^25,33,34,37

Recently, microfluidic experiments have played an important role in advancing foam characterization and optimization for EOR applications. Oil is known to destabilize foams by disrupting films, limiting their effectiveness in EOR. To address this, recent microfluidic studies have examined foam stability under varied conditions—including oil composition,^38,39 gas component,^40,41 brine salinity,^39,42 pressure and temperature,^41,43 and surfactant type and concentration.^44,45 Foam coarsening was observed to increase at higher gas ratios, higher temperatures, lower surfactant concentrations, or when using gas types with higher solubility.⁴¹ Bubble population dynamics were shown to depend on the interplay between foam generation and decay.⁴⁶ Moreover, Lopes et al.⁴⁷ recommended adding N₂ to CO₂ gas to decrease bubble coarsening through gas diffusion.

Foam texture—the size, shape, and distribution of gas bubbles—has emerged as a key parameter, besides velocity, governing flow resistance in porous media.^48–51 Foam texture is influenced by pore size distribution, gas fraction, surfactant type, and bubble generation and decay dynamics.^1,17,52 The transition between “wet” (low gas fraction) and “dry” (high gas fraction) foams is regulated by capillary pressure thresholds and the corresponding water saturation levels.⁵³ At high liquid saturation, foam bubbles are small and easily transported, while at reduced liquid contents and low flow velocities, bubbles tend to coarsen, decreasing foam strength. Microfluidic studies have visualized this transition and its influence on pressure drop and flow redistribution.^51,54 Foams with higher gas ratios tend to enhance foam viscosity and improve propagation and diversion efficiency.^42,45,49 Wang et al.⁵⁴ showed that increasing the gas ratio up to 0.8–0.9 increases bubble size, lamellae density, and apparent viscosity. However, some studies indicate that at high gas ratios, coarser foam texture may reduce viscosity.⁵¹

Microfluidic visualizations have consistently revealed that fluid displacement by foam is an efficient method to suppress viscous/capillary fingering, thereby improving sweep efficiency.^5,6,22,55 In homogeneous devices, only CO₂ foam generated a stable displacement front, outperforming liquid, gas, and supercritical CO₂ injections.⁵⁵ Similarly, air foam injections produced better sweep efficiency than air gas, primarily due to enhanced mobility control and delayed breakthrough.⁶ Conn et al. showed that foam could displace oil from both high- and low-permeability zones and intermediate fractures, achieving up to 80% oil recovery, in contrast to 50% and 20% for water and WAG (water-alternating-gas), respectively.⁵ Surfactant and surfactant-nanoparticle foams can enhance recovery compared to water injection, reaching 77% and 94%.⁵⁶ These findings are also supported by Guo et al., who observed foam-enhanced recovery gains exceeding 30% over water and gas flooding alone.^57,58

Foam injection also serves as an effective means of mitigating preferential flow pathways commonly encountered during gas injection in heterogeneous porous media.^5,59 Microfluidic experiments reveal that foam bubbles preferentially occupy and block high-permeability channels, promoting fluid redistribution into adjacent, previously unswept regions.^46,60 This selective occupation is governed by capillary pressure (P_c), which increases as the pore size decreases. When P_c exceeds the entry pressure for gas, the non-wetting phase is excluded from low-permeability pores, thereby encouraging aqueous phase invasion and improved sweep.²² This interplay between capillary thresholds and the foam structure supports foam's unique ability to propagate in heterogeneous porous zones.⁵⁸

Despite these advances, pore-scale displacement of oil by CO₂ foam remains insufficiently understood due to the complex interfacial interactions between oil, gas–liquid foam structures, and solid surfaces within porous matrices. A comprehensive investigation into how injection conditions and foam formulation affect displacement dynamics is critical for optimizing foam performance across EOR, aquifer remediation, and CO₂ sequestration scenarios. In this study, we employ a heterogeneous microfluidic platform to visualize and quantify CO₂ foam transport and oil displacement behaviors, focusing on the role of the gas fraction and displacing fluid type in controlling sweep efficiency and mobility stabilization.

2 Experimental

In this study, we conducted three distinct sets of experiments designed to quantify and understand the structure, flow transport, and oil displacement performance of CO₂ foam in a heterogeneous microfluidic porous medium. First, we characterized foam bubble generation at the microfluidic flow-focusing junction at varying gas injection pressures (P_g) and liquid flow rates (Q_l), prior to entry into the microfluidic porous medium. Second, we investigated foam transport dynamics within the porous medium (without a residing fluid), examining the influence of the gas ratio (R_g) on fluid distribution, bubble velocity, and overall foam velocity. Third, we evaluated oil displacement performance using SDS solution, CO₂ gas, and CO₂ foam, with emphasis on how variations in Q_l, P_g, and R_g impact breakthrough behavior and recovery efficiency. Their respective results are presented in sections 3.1–3.3.

2.1 Microfluidic chip design and fabrication

A heterogeneous porous medium (27 800 μm × 12 400 μm) was designed, comprising two homogeneous regions with distinct porosities and permeabilities (see Fig. 1). All the pillars are circular and arranged in an orthorhombic lattice to represent rock grains. The high-permeable (HP) region contains pillars of diameter D_p = 210 μm, while the low-permeable (LP) region features larger pillars of D_p = 250 μm. The throat distances, D_th, are 60 μm and 21 μm, corresponding to porosities of 45% and 23% in the HP and LP regions, respectively. The geometric parameters—pillar size, pore body, and throat dimensions—were selected based on petrophysical data available in the literature for sandstone reservoirs^61–69 to closely mimic natural rock structures.

Fig. 1 Microfluidic device: the schematic of the microfluidic chip and fluid injection/production ports. The liquid phase is in blue, gas in white, oil in yellow, and solid pillar in gray. The definitions of pillar diameter (D_p), throat size (D_th), bubble size (D_b), and pore area (A_p) are provided. All scale bars represent 500 μm.

The total average permeability, K_avg, of the porous medium was estimated by injecting pure water at various flow rates (Appendix 5.1), measuring the pressure drop using a differential pressure transducer (Omega 0–50 psi), and applying Darcy's law. The average permeability of the medium (including both low and high permeable regions) was estimated to be 6.07 Darcy. As the permeability of each homogeneous layer is proportional to the porosity and hydraulic radius,⁷⁰ the permeability ratio between the layers was determined to be 5.8. The hydraulic radii were estimated based on wetted surface areas and pore volumes. Considering the permeability averaging rule for parallel layers, the permeabilities of the low- and high-permeable layers were estimated to be 1.8 and 10.1 Darcy, respectively (see Appendix 5.1 for detailed calculations).

Foam was generated in the same chip using a flow-focusing geometry placed upstream of the porous medium (Fig. 1). Here, a surfactant solution and CO₂ gas were co-injected and converged at a 60 μm-wide orifice to produce bubbles. A serpentine (winding) channel segment was incorporated between the foam generation region and the porous medium. The bubbles traveled through the winding channels connected to a drain channel, which remained open during stabilization to prevent unsteady foam entry into the porous structure. The winding channels served both as a temporary reservoir for the generated bubbles and as a hydraulic resistance element that isolates the foam generator from downstream pressure variations, thereby facilitating stable bubble formation and delivery. All flow and distribution channels, including those leading to the porous medium, were uniformly 200 μm wide.

The microfluidic pattern was fabricated on a 6′′ silicon wafer coated with a positive photoresist (AZ P4330) using maskless photolithography (MLA 150, Heidelberg Instruments), followed by uniform etching to a depth of 100 μm (ICPRIE – Oxford Estrelas) to create the microfluidic mold. A soft lithography approach was then used to cast and cure PDMS on a hydrophobized mold (treatedwith Trichloro(octadecyl)silane, Sigma Aldrich) to produce the chips.

All inlet and outlet ports (1.5 mm diameter) were punched using a coring tool (EMS Rapid-Core 1.5) to facilitate tubing connections. The PDMS slab was bonded to a glass slide to enclose the channels. Both surfaces were cleaned with isopropanol and exposed to oxygen plasma to functionalize them with silanol groups, enabling strong, permanent PDMS–glass bonding.

2.2 Materials and methods

The surfactant used to generate foam was 0.5 wt% sodium dodecyl sulphate (SDS), dissolved in Milli-Q water (resistivity > 20 MΩ cm). To enhance visual contrast with the oil phase and improve image processing, 1.0 wt% Nile Blue dye was added to the aqueous solution. The gas phase used was CO₂ gas with 99.99% purity (Linde). For oil displacement experiments, a paraffinic oil (ISOPAR M, ExxonMobil) with a viscosity of 3.2 cP was employed. We initially conducted exploratory experiments with oils of higher viscosities (10 cP and 35 cP); however, the higher injection pressures required to saturate and subsequently displace these oils exceeded the mechanical integrity of the PDMS–glass bonding, resulting in de-bonding and fluid leakage. To address this limitation in future studies, we plan to employ alternative microfluidic materials—such as fully glass-based or silicon–glass hybrid chips—with tailored porosity and permeability designs. These platforms offer greater pressure tolerance while maintaining optical accessibility, enabling systematic investigation of foam–oil displacement dynamics over a broader range of viscosity ratios.

The contact angles of the aqueous phase and the oil droplets on PDMS were measured in the presence of air to be 86° and 18°, respectively (Appendix 5.3). On bare glass, the contact angles were 32° for the aqueous phase and 12° for the oil, indicating that the porous medium—composed of a glass base and PDMS features—was overall more hydrophobic (or oil-wet). The wettability of the microfluidic chip remained unchanged after exposure to both the SDS solution and ISOPAR M oil, indicating no appreciable alteration from surfactant adsorption or oil contact. In addition, the interfacial tension (IFT), σ, between various fluids was measured using pendant drop tensiometry.^71,72 The results showed the following values for the SDS solution and oil used: σ_SDS,air = 36 mN m⁻¹, σ_oil,air = 16 mN m⁻¹, σ_SDS,CO2 = 38 mN m⁻¹, σ_oil,CO2 = 18 mN m⁻¹, and σ_SDS,oil = 17 mN m⁻¹ (see Appendix 5.4).

The stability of foam bubbles could be evaluated using the IFT values to calculate the entering (E), bridging (B), and spreading (S) coefficients following established criteria.^73,74 According to the results in Table 1, ISOPAR M oil can enter the gas–liquid interface (E > 0), rupture it (B > 0), and spread along it (S > 0), leading to foam destabilization.

Table 1 The entering (E), bridging (B), and spreading (S) coefficients of the SDS solution–CO₂–ISOPAR M oil system, estimated using the interfacial tensions (σ) obtained in Fig. 14

Coefficient	Definition	Value
Entering (E)	σ SDS/CO2 + σSDS/oil − σoil/CO2	37.1
Bridging (B)		1400.2
Spreading (S)	σ SDS/CO2 − σSDS/oil − σoil/CO2	2.5

---

The blue aqueous solution and ISOPAR M oil were injected using a syringe pump (Chemyx) operated in constant flow rate (Q_l) mode, and the gas flow was controlled via a pressure-regulated mass flow controller (ALICAT) operated in constant pressure mode (P_g), producing a corresponding gas flow rate (Q_g). Microfluidic valves (Upchurch Scientific) were used to open or close the outlet and drain ports as needed.

2.3 Image processing

The microfluidic chip was mounted and secured on an inverted microscope (Zeiss Axio Observer 7 Materials) connected to a high-speed camera (HSC, Photron Nova S9) for experimental visualization. The captured videos were processed and analyzed using Fiji (ImageJ) and MATLAB. In all image analyses in this work, a 2D flow condition was assumed.

Different magnifications (2.5×, 5×, and 10×) and frame rates (125–15 000 fps) were used depending on the resolution required for each specific analysis. Images were converted to grayscale, and fast Fourier transform (FFT) band-pass filtering^75,76 was applied for shading correction. Thresholding was then used to isolate region of interest, such as bubble areas or fluid saturation. Each pixel was assigned to a single fluid phase.

Bubble sizes were calculated assuming a circular geometry, using the relation , where D_b is the bubble diameter and A_b is the bubble area measured from the 2D images (see Fig. 1). The main oil displacement experiments were carried out at a low magnification of 2.5× to maximize the field of view, which limited resolution and contrast, particularly in the low-permeable area.

3 Results and discussion

In this study, we evaluated the dynamics of CO₂ foam flow and assess its potential for EOR application in porous media. Therefore, the results are presented in three parts: a) characterization of bubbles generated at the flow-focusing foam generator, b) analysis of foam size and velocity in the porous medium, and c) oil displacement by SDS solution, CO₂ gas, and CO₂ foam in the porous medium.

3.1 Foam bubble size characterization

To generate different foam textures using the flow-focusing geometry, gas was injected at constant pressure while varying the liquid flow rate. Two different gas pressures—P_g = 22 and 26 psi—were applied, and the liquid flow rates (Q_l) ranged from 25–300 μL min⁻¹ and 30–300 μL min⁻¹, respectively. All experiments in this part were captured at a frame rate of 5000 fps.

During these experiments, the outlet port remained open, which allowed the bubbles to travel through the winding channels and into the porous medium, while the drain port was closed. The resulting bubble trains remained stable within the region of interest (ROI shown in Fig. 2 snapshots)–flow-focusing foam generation area–and throughout the observation period. Each experiment was repeated 3 to 4 times to ensure reproducibility.

Fig. 2 Foam bubble characterization: (a) variation of average bubble diameter D̄_b with the liquid flow rate (Q_l) at P_g = 22 psi and P_g = 26 psi. (b) Effect of the gas-to-liquid flow rate ratio (Q_g/Q_l) on the dimensionless bubble diameter (D̄_b/W_orifice), where W_orifice = 60 μm is the flow-focusing geometry. Equations of the fitted power-law trends (dashed lines) are shown, where x = Q_g/Q_l and y = D̄_b/W_orifice, respectively. The R² values for the fitting curves are 0.98 and 0.97 for P_g = 22 psi and P_g = 26 psi cases, respectively. Representative experimental snapshots of produced bubble trains are labeled (I–V), with border colors matching the corresponding data points.

To calculate the average diameter of the bubbles for each condition, all bubble diameters (D_b) in a single representative frame were averaged and reported as D̄_b. Fig. 2(a) shows the variation of D̄_b as a function of Q_l for both gas pressures, indicating that D̄_b decreases with increasing Q_l at constant gas pressure.

Fig. 2(b) presents the relationship between the normalized bubble size, D̄_b/W_orifice, and the flow rate ratio, Q_g/Q_l, at both gas pressures, along with representative experimental snapshots. Here, D̄_b is the average bubble diameter obtained in Fig. 2(a) and W_orifice = 60 μm is the constant width of the orifice in the flow-focusing generator.

At a constant gas flow rate—Q_g = 0.13 mL min⁻¹ and 0.19 mL min⁻¹ at gas pressures P_g = 22 psi and P_g = 26 psi, respectively—increasing Q_l resulted in smaller bubbles. Higher liquid flow rates impose greater shear stress on the (constant-pressure) gas thread, increasing the bubble pinch-off frequency at the expense of bubble size.^77–79 Moreover, higher gas pressure generated larger bubbles. The bubble diameters ranged from 30–270 μm.

Previous studies using microfluidic flow-focusing bubble generators with different dimensions have proposed a power-law relation, D̄_b/W_orifice = a(Q_g/Q_l)^b, to predict the bubble size.^78,80 Similarly, power-law fits to our data (with R² values of 0.98 and 0.97) in Fig. 2(b) confirm that the bubble size can be reliably predicted under these injection conditions. The constants a and b appear to depend on the flow-focusing geometry, surface wettability, and fluid properties.^78,80 Both fitted curves exhibit similar power-law exponents of 0.6, which aligns with the reported range (0.3 < b < 1.0) in the literature.⁸⁰ This power-law exponent (b) is primarily determined by the flow focusing geometry,⁸⁰ while the coefficient a appears to be more sensitive to variations in the injection gas pressure.

3.2 Foam flow in the microfluidic porous medium

The flow of generated bubbles in the porous medium was captured to investigate the distribution of fluid phases in terms of 2D areal fractions—CO₂ gas (F_g), both trapped and mobile—and the bubble velocities. From this point forward, all experiments were conducted at a constant gas pressure of P_g = 22 psi. Liquid flow rates were varied from 10 to 400 μL min⁻¹, resulting in gas ratios (R_g) ranging from 0.25 to 0.93. The gas ratio R_g is defined as the ratio of the gas flow rate (Q_g) to the total flow rate (including both gas and liquid injection rates):During these experiments, the drain line was closed, while the outlet remained open to permit flow in the porous medium (see Fig. 1). To distinguish and track the bubbles, we selected the image magnification of 5× and the frame rate of 15 000 fps. Due to camera memory limitation (2 second videos at this frame rate), videos were recorded when the fluid distribution was established in the porous medium. The middle snapshot of each video, e.g., Fig. 2(a-I-left), was analyzed to calculate the distribution fractions. During the period of the captured videos, fluid saturation barely changed, as confirmed by the representative videos in the supplementary information. Each experiment was repeated 3–4 times for reproducibility.

3.2.1 Effect of the gas ratio on fluid distribution fractions. In the sample raw image in Fig. 3(a-I-left), the gas, aqueous, and solid phases are shown in white, blue, and light gray, respectively. In the processed image of Fig. 3(a-I-middle), the areal fraction occupied by the gas, F_g, is shown in yellow and defined as: F_g = (pore area occupied by the gas phase)/(total pore area).

Fig. 3 Foam distribution in the microfluidic porous medium: (a-I) the raw image (left) processed to segregate the areal gas fraction, F_g, shown in yellow (middle), and differentiate the trapped and moving gas fractions shown in cyan and purple, respectively (right) at R_g = 0.57. (a-II) and (a-III) indicate trapped and moving F_g at R_g = 0.30 and R_g = 0.87, respectively. The corresponding data points are labeled as in (c). (b) The variation of gas fraction, F_g, in the porous medium as a function of gas ratio (R_g). (c) Effect of gas fraction on gas trapping in the porous medium.

As revealed in Fig. 3(b), increasing R_g from 0.25 to 0.93 resulted in higher gas fractions in the porous medium. This is expected because with increased gas ratios, less water volume is introduced into the device, leading to greater F_g. The data indicate that between 50% and 93% of the pore area is occupied by the gas phase using different tested ranges of wet to dry foams. A linear fit with a coefficient of determination R² = 0.74 indicates a similar increasing trend in this range of gas ratios.

Each experiment revealed the coexistence of mobile and trapped gas bubbles: some bubbles flowed through the porous structure due to hydrodynamic pressure drop, while others became immobilized, trapped in the pore structures. Thus, we categorized the gas phase into trapped and moving components, and defined their fractions as follows: the trapped (or moving) gas fraction is the ratio of the area occupied by trapped (or moving) gas to the total area occupied by the gas phase.

In Fig. 3(a-I-right), the yellow gas region is divided into purple (moving) and cyan (trapped) regions (to distinguish the moving and trapped bubbles). As shown in Fig. 3(c), increasing R_g resulted in a larger proportion of trapped gas. At R_g = 0.25, the trapped gas fraction is 0.23, which increases to 0.91 at R_g = 0.93.

Beyond R_g = 0.76, the trapped F_g does not change significantly and reaches a relative plateau value close to 0.90. Fig. 3(a-I-right), (a-II), and (a-III) correspond to the labeled data points in Fig. 3(c). As the foam propagates and accumulates in porous media, certain pathways become increasingly blocked, leading to local increases in flow resistance. To sustain fluid transport, flow is redirected through less obstructed routes. In all cases, a central pathway—the horizontal channel within the high-permeability zone closest to the main inlet flow channel—remained unblocked and allowed bubble movement. This persistent open channel likely represents the least resistive path due to its direct alignment with the inlet, enabling preferential flow. This path contributes to ≈10% of the pore volume and accounts for the observed plateau.

Even at the highest moving gas fraction of 0.77 at R_g = 0.25, no moving bubbles were detected in the low-permeable zone. This implies that although the foam bubbles overcome the critical capillary pressure to enter the pore areas in both zones, mobile bubbles exclusively traveled through the high-permeable zone, which offered the path of least resistance.

Across all gas ratios, D̄_b/D_th > 1 in both zones, indicating that bubbles generated by the flow focusing geometry needed to shrink or break to enter the pore throats of both zones. This ratio is smaller (0.5–4.7) for the high-permeable zone than in the low-permeable zone (1.5–13.3). As discussed, foam bubbles in the microfluidic porous media are either trapped or mobile. The average sizes of moving bubbles, D̄_mb, and trapped bubble areas Ā_tb were obtained by averaging the corresponding detected bubbles from the same representative frames used in Fig. 3(a–c).

The comparison between D̄_mb and D_th, in addition to the visualized results, reveals that regardless of R_g, the bubbles flow only through the high-permeable zone as the less resistive path and they do not need to shrink much (maximum D̄_mb/D_th is around 2). In contrast, entering the low-permeable zone would require bubbles to shrink by a factor of 4–8. This preferential foam flow aligns with previous observations that foams tend to migrate through a larger pore pathway.⁸¹

Finally, a comparison of trapped bubble areas (Ā_tb) with local pore areas (A_p) showed that trapped bubbles conformed closely to the pore size, with A_tb/A_p ∼ 1 in both high-permeable and low-permeable zones. Being restricted by the pillar walls, this bubble alignment phenomenon has been observed⁸¹ and even used for spatial patterning.⁸²

The solubility and diffusivity of CO₂ in water at atmospheric pressure and ambient temperature (14.7 psi and 25 °C) are relatively low.^83–86 Under these conditions, the estimated Peclet numbers in our system are considerably high (>1600), indicating that convective transport strongly dominates over diffusion. In addition, the observed foam bubble sizes remain relatively constant during flow, further supporting that CO₂ dissolution in water via diffusion is negligible under our experimental conditions.

3.2.2 Effect of the gas ratio on velocities of bubbles and foam. To estimate the velocity of the moving bubbles, V_mb, at different R_g values, the recorded videos were converted to binary images and enhanced for contrast (see Fig. 4(a)). The wrMTrck plugin in Fiji (ImageJ) was used to identify and track bubble movements.

Fig. 4 Effect of the gas ratio on bubble velocity: (a) the representative snapshots of moving bubbles (white) in the microfluidic porous media at gas ratios R_g = 0.25 (I), 0.63 (II), and 0.93 (III). The circles are the pillars of the high-permeable zone. (b) The variation of 〈V̄_mb〉 (left axis) and D̄_mb (right axis) as a function of gas ratio, R_g. The corresponding labeled snapshots in (a) are labeled similary in (b).

The displacements of all the bubbles' centroids were tracked and measured throughout the image sequences. The velocities of moving bubbles were then averaged to obtain a single mean value, 〈V̄_mb〉. Bubbles were classified as mobile or trapped based on their measured velocity. A threshold of 1 pixel per frame was used to distinguish the two populations: bubbles with velocities below this value were considered stationary (trapped), while those with higher velocities were classified as mobile.

Fig. 4(a and b-I to III) show snapshots of three representative moving foam bubbles at different R_g values along with the corresponding average velocity results. As shown, 〈V̄_mb〉 generally falls in the range of 75–95 mm s⁻¹, except at R_g = 0.93, where 〈V̄_mb〉 = 130 mm s⁻¹. In contrast to uniform bubbles generated at the flow-focusing junction (Fig. 2(a)), moving bubbles within the porous medium displayed less uniformity and were more influenced by the local porous structure. As shown in Fig. 4(b), the mean size of the moving bubbles, D̄_mb, remains relatively unchanged (0.074–0.091 mm) for R_g ≤ 0.76, but increases at higher gas ratios. Bubbles at R_g ≤ 0.76 are well separated and more circular (Fig. 4(a-I and II)), whereas R_g > 0.76 caused connected large bubbles (Fig. 4(a-III)). Moreover, bubble velocity appears to have a weak dependence on the gas ratio, as no clear correlation is observed in 〈V̄_mb〉 with increasing R_g from 0.25 to 0.93 (see Fig. 4(b)).

In addition to hydrodynamic pressure, the Jamin effect plays a significant role in governing foam flow dynamics within porous media.^87–89 When a bubble approaches a pore constriction, its leading edge compresses while the trailing edge remains relatively undeformed, producing a curvature difference that results in additional capillary resistance.⁸⁸ Two factors help reduce this resistance: (1) smaller bubbles require less energy to deform, weakening the Jamin effect, and (2) higher gas entrapment increases the pressure gradient across the porous medium, enhancing foam bubble velocity. These competing factors may explain the lack of a consistent trend in 〈V̄_mb〉 versus R_g.

The average foam velocity (V̄_f), incorporating both moving and trapped bubbles, was estimated as:

where F_g,mb and F_g,tb are the fractions of moving and trapped bubbles, obtained from Fig. 3(c); V̄_f denotes the average foam velocity; 〈V̄_tb〉 represents the trapped foam velocity, which is assumed to be zero (as its traveling distance is less than 1 pixel between consecutive frames); 〈V̄_mb〉 is the velocity of moving bubbles as obtained in Fig. 4(b).

As revealed in Fig. 5(a-right axis), the overall foam velocity V̄_f decreases from approximately 73 mm s⁻¹ to 7 mm s⁻¹, with increasing gas ratio. Comparing trends in V̄_f and F_g,mb reveals that bubble entrapment has a stronger influence on overall foam velocity than on the velocity of individual moving bubbles (Fig. 4(b)). In other words, a less mobile gas fraction results in slower foam propagation. It is important to emphasize that V̄_f accounts for all bubbles, while 〈V̄_mb〉 only represents the moving ones.

Fig. 5 Foam velocity: variation of overall foam velocity (V̄_f) with the gas ratio (R_g) and the effects of (a) foam trapping (F_g,mb), (b) foam texture (n_f), and (c) foam size (D̄_tf). The right axis of all plots are assigned to the overall foam velocity (V̄_f); the left axes display the respective independent variables.

Furthermore, foam velocity variation with R_g is compared to foam texture—defined as the number of bubbles per unit area (n_f)—and to the mean bubble diameter of the same foam (D̄_f) in Fig. 5(b) and (c), respectively. These parameters account for both moving and trapped bubbles. The results show that overall, larger bubbles (or lower foam texture) are associated with slower foam movement. The mean bubble diameter D̄_f ranges from 740 μm to 1000 μm, while n_f ranges from 70 to 104 mm⁻². This inverse relationship reflects that the increased gas ratio leads to larger bubbles and lower foam texture, which reduces foam mobility. These trends are clearly seen in the comparison between Fig. 5(b) and (c).

In summary, a higher foam gas ratio results in a greater fraction of the porous medium being occupied by CO₂ gas, much of which is trapped rather than flowing. Due to capillary barriers, mobile bubbles consistently flow through the high-permeable zone, while trapped bubbles accumulate in both high- and low-permeable zones depending on R_g. This increased entrapment results in larger average bubble size and lower foam texture, thereby reducing the overall foam velocity.

3.3 Oil displacement experiments by gas, aqueous solution, and foam

In this section, we investigated the enhanced oil recovery (EOR) performance of three displacement strategies by (1) CO₂ gas, (2) SDS aqueous solution, and (3) CO₂ foam. All experiments were conducted using a 2.5× magnification and 125–250 fps, which provided a field of view of 4096 μm × 4096 μm. The displacement was monitored near the center of the porous medium, away from the outlet, to minimize capillary end effects.^55,90 Sudden reductions in capillary pressure near the outlet may lead to local accumulation of wetting fluid and distortion of fluid distribution analysis.

ISOPAR M oil was injected from the outlet at a constant flow rate of 100 μl min⁻¹. The drain port remained open to allow the oil phase to exit, while the gas and liquid injection ports were closed (using Idex shut-off valves) to prevent oil from entering the foam generation area. Once the microfluidic porous medium was fully saturated with oil, the respective injection port(s) were opened depending on the displacement fluid. Continued oil injection maintained a back pressure to prevent entry of unstable displacing fluid into the porous space.

The displacing fluid flowed out through the drain port until stable injection conditions were reached, i.e., constant gas pressure (for CO₂ gas injections), constant liquid flow rate (for SDS solution injections), or steady bubble regime (for foam injections). Once stabilized, the outlet was opened to atmospheric pressure, the drain port was closed to redirect the displacing fluid into the porous medium, and data acquisition began.

The tested conditions included five CO₂ gas pressures (15–28 psi), five SDS solution flow rates (5–100 μm min⁻¹), and six foam gas ratios (R_g = 0.57–0.98 at P_g = 22 psi). Each experiment was repeated 3 to 4 times to ensure reproducibility.

The image sequences in Fig. 6 illustrate oil displacement by the three fluids: (a) SDS aqueous solution, (b) CO₂ gas, and (c) CO₂ foam. In all cases, displacement initiated in the high-permeable (HP) zone (labeled “HP” in Fig. 6). The SDS solution (Q_l = 25 μl min⁻¹) forms a relatively stable displacement front and exhibits a delayed breakthrough at T_bt = 1.33 s, as shown in Fig. 6(a). Although the high-permeable (HP) zone was quickly swept (within 2 s), the solution fails to enter the low-permeable (LP) zone, even after extended injection (up to 78.89 s).

Fig. 6 Oil displacement dynamics: image sequences of oil displacement by (a) SDS solution (25 μl min⁻¹), (b) CO₂ gas (P_g = 22 psi), and (c) CO₂ foam generated with SDS solution (R_g = 0.87 at P_g = 22 psi). Displacement progress is from left to right; colored regions represent displacing fluids, and black regions include oil and pillars. Red frames indicate breakthrough times (T_bt). Raw images at T_bt are shown on the right. “HP” and “LP” refer to high- and low-permeable zones. The raw images of the breakthrough times are shown on the right side of each row.

In contrast, CO₂ gas injection (at P_g = 22 psi) led to viscous fingering in the direction of flow (left to right), producing an early breakthrough at T_bt = 0.02 s (Fig. 6(b)). Although this unstable gas–oil interface delayed complete oil displacement in the HP zone (24.19 s), gas was able to penetrate the LP zone, achieving nearly complete oil recovery after a prolonged injection period of 161.89 s. The improved sweep may be attributed to a larger pressure, estimated using Darcy's law (along with the chip geometry, fluid properties,^91,92 and operating conditions) to be 0.23 psi for SDS solution vs. 7.30 psi for CO₂ gas (see Appendix 5.1). Viscous fingering was less severe in the LP zone due to its significantly lower permeability, ≈6× lower than that of the HP zone.

Fig. 6(c) shows the oil displacement by CO₂ foam (R_g = 0.87), which achieved a relatively stable front (in the HP zone) similar to SDS solution, but with deeper penetration (in the LP zone) like gas. CO₂ foam broke through at T_bt = 1.67 s and recovered oil from the LP zone in just 30.02 s, about 5× faster than the pure gas injection. This better performance stemmed from foam lamellae trapping gas bubbles, thereby creating additional flow resistance and suppressing viscous fingering. Consequently, CO₂ foam exhibits both delayed breakthrough and enhanced oil recovery efficiency, combining the benefits of stable propagation (in the HP zone) and deeper penetration (in the LP area).

To elucidate the fluid distributions during the microfluidic EOR process, oil displacement experiments under the same conditions as in Fig. 6 were repeated at a higher magnification (5×). The resulting final fluid distributions (after the displacement) are shown in Fig. 7. Consistent with previous observations, the SDS solution (blue) failed to displace oil (gray) in the LP zone (Fig. 7(a)). In contrast, both CO₂ gas and foam effectively swept oil from the field of view (Fig. 7(b) and (c)). Fig. 7(c) and (d), acquired at 5× and 10× magnifications, respectively, reveal that the high- and low-permeable zones are mainly displaced by foam and liquid, respectively. This behavior may be attributed to the fact that the non-wetting gas phase overcomes opposing capillary forces in the high-permeable zone more readily. Similar trends have been reported in the literature for heterogeneous porous media, where lower-permeability regions tend to become liquid-rich, while injected foam bubbles preferentially fill the higher-permeability zones.^5,22

Fig. 7 Fluid distribution at oil displacement: magnified (5×) snapshots showing fluid configurations after oil displacement by (a) SDS solution (at t = 9.7 s) showing no penetration of blue SDS solution into the (upper) low-permeable zone, (b) CO₂ gas (at t = 94.3 s) showing penetration of the CO₂ gas into both the high and low-permeable zones, and (c) CO₂ foam (at t = 88.2 s) showing penetration of CO₂ foam into the high-permeable zone and SDS solution into the low-permeable zone. (d) Higher magnification (10×) from the same experiment as (c).

Breakthrough time (T_bt) is a key performance metric in EOR field applications, as it correlates with sweep efficiency and operational feasibility for a given strategy.^93–95 To ensure cost-effectiveness, the volume of produced gas (gas–oil ratio) and water (water cut), relative to oil, must remain sufficiently low.

Fig. 8 presents the T_bt and oil recovery factor at breakthrough (RF_bt) for all the displacements.

Fig. 8 Oil displacement efficiency at breakthrough: the time of the breakthrough of the displacing fluid (T_bt) and the corresponding oil recovery factor (RF_bt) for (a) SDS aqueous solution, (b) CO₂ gas, and (c) CO₂ foam generated with SDS solution. Each data point represents 3–4 repeated experiments.

For the SDS aqueous solution, T_bt decreases from 6.6 s to 1.0 s as the injection flow rate (Q_l) increases from 5 μl min⁻¹ to 100 μl min⁻¹ (Fig. 8(a)). This is expected, as higher flow rates result in greater hydrodynamic pressure drops, leading to faster fluid propagation and consequently earlier breakthrough. However, a higher Q_l, with a shorter T_bt, reduces sweep efficiency, with an RF_bt as low as 40% at Q_l = 100 μL min⁻¹. The maximum oil recovery at breakthrough for SDS solution injection is 52%. Contact angle measurements confirm that the microfluidic chip is hydrophobic. Even at the lowest flow rate (Q_l), the displacement front remains stable, and no lateral capillary fingering is observed. This suggests that oil displacement by the SDS solution is primarily facilitated by desirable viscosity contrast (M < 1) rather than capillary forces. A comparative analysis of the displacement behavior using pure water versus SDS solution would provide further insight into the role of surfactants and capillary effects. However, such a comparison lies beyond the scope of the present study.

Injection of CO₂ gas accelerates the breakthrough times from 0.18 s down to 0.02 s when P_g = 15–28 psi (Fig. 8(b)). In all gas injection experiments, even at the lowest injection pressure P_g, gas fingers through the oil toward the outlet, as shown in Fig. 9(a). Viscous fingering reduced the sweep efficiency, with RF_bt being reduced to 12–18%, significantly lower than that of SDS solution and CO₂ foam with stable displacement fronts. At P_g = 15 psi, the oil recovery factor by CO₂ gas is low despite a delayed breakthrough, as shown in Fig. 9(b). This is because the resulting pressure drop (≈0.3 psig) is too small to effectively push the gas through the porous medium, slowing down the displacement process. Observations show that at this pressure, CO₂ cannot enter the LP zone, likely due to its higher entry capillary pressure. Capillary pressure is defined as P_c = P_non-wet − P_wet = 2γκ_c cos(θ), where γ is the interfacial tension, θ is the contact angle, and is the average curvature of the porous medium; r₁ and r₂ are the two principal pore radii or radii of curvature of the meniscus. Here, and ; D_th is the diameter of the pore throat (60 and 21 μm for the high- and low-permeable zones, respectively), and H is the height of the microfluidic chip (100 μm). Assuming a constant contact angle (cosθ = 1) for both regions, and using the measured pore throat radii, the estimated capillary entry pressures (for gas displacing oil) are approximately 960 Pa for the high-permeability region and 2074 Pa for the low-permeability region. According to Fig. 9(b), it seems that P_g = 15 psi is insufficient to overcome this 1114 Pa pressure difference between the two zones and unable to invade the LP zone. However, when the gas pressure increases to 20 psi, CO₂ gas begins to penetrate the LP zone, suggesting that the gas critical entry pressure (P_oil + P_c) lies between 15 and 20 psi.

Fig. 9 Fluid–fluid interactions: (a) viscous fingering of CO₂ gas (red frames) through the oil phase within the HP zone at P_g = 15 psi (left) and P_g = 28 psi (right). (b) CO₂ cannot overcome the capillary pressure of the low-permeable zone (top part) to enter and sweep the oil at P_g = 15 psi but at P_g = 20 psi (right), CO₂ enters the low permeable zone. (c) Foam destabilization at the contact of oil under the effect of gas ratios of 0.98 (left) and 0.57 (right) within the HP zone. (d) Effect of the presence of gas on the displacement front of SDS solution (left) and CO₂ foam injection within the HP zone (right). The liquid flow rate in both cases is 22 μl min⁻¹.

Foam flooding significantly improves recovery. As R_g increases from 0.57 to 0.98, T_bt decreases from 1.2 s to 0.54 s, while RF_bt ranges from 44% to 70% (Fig. 8(c)). Foam injection at P_g = 22 psi improves the breakthrough time (T_bt = 0.03 s) by over an order of magnitude compared to gas and enhances recovery by 27% or more. Our results reveal that replacing CO₂ gas injection with CO₂ foam at 0.57 ≤ R_g ≤ 0.98 is beneficial by stabilizing the displacing front, delaying the breakthrough time, and increasing the oil sweep efficiency.

In contrast to the SDS solution, increasing Q_l in foam flooding (i.e., lowering R_g) improved oil displacement with longer T_bt and higher oil RF_bt. Two factors may contribute to this trend. First, contact between the foam and oil can destabilize the foam, leading to segregated gas slugs, mild viscous fingering, and earlier breakthrough (Fig. 6(c), 1.67 s). As shown in Table 1, ISOPAR M oil exhibits positive entering (E > 0), bridging (B > 0), and spreading (S > 0) coefficients, indicating that the oil can penetrate the gas–liquid interface, rupture the foam lamellae, and spread along the interfacial film. These characteristics suggest a strong tendency for foam destabilization. At the foam–oil front, this leads to lamella rupture, gas slug release, and the onset of gas fingering, as illustrated in Fig. 9(c). The destabilization reduces foam's ability to maintain a stable displacement front, resulting in shorter breakthrough times (T_bt) and lower oil recovery at breakthrough (oil RF_bt), as shown in Fig. 8(c). This effect can become more pronounced at higher R_g due to the increased gas content and reduced availability of the foam stabilizer (e.g., surfactant). Moreover, as discussed in section 1, foam flow in porous media can occur in either the high R_g (high-quality) or low R_g (low-quality) regimes. At a constant gas flow rate (22 psi here), the low R_g regime is relatively insensitive to variations in the liquid flow rate, resulting in minimal changes in pressure drop. In contrast, upon transitioning into the high R_g regime, even small reductions in the liquid flow rate lead to significant decreases in pressure drop.⁵³ Numerous studies have observed shear-thinning behavior in the high R_g regime.^51,96,97 A lower foam viscosity in this regime causes shorter T_bt, potentially contributing to the decreasing trend observed in Fig. 8(c) beyond R_g = 0.87. Second, gas injection pressure enhances the hydrodynamic driving force, accelerating the displacement front while maintaining stability and improving recovery (Fig. 9(d)). Under the conditions tested, foam destabilization likely dominates the observed trend in Fig. 8(c), while the added driving force explains the superior performance of foam over SDS solution, with higher recovery achieved despite shorter T_bt.

Fig. 10 displays the dynamic performance of oil production (recovery factor, RF, versus time) by the three representative displacing fluids: SDS solution (100 μl min⁻¹, blue), CO₂ gas (22 psi, black), and CO₂ foam (R_g = 0.57, red). The oil RF at each time was calculated by the ratio of the displaced to the initial oil-saturated area. Note that the SDS curve follows the bottom time axis, while the gas and foam curves correspond to the top axis. The horizontal gray lines indicate the recovery limits for the high-permeable zone (∼70%) and full displacement (100%).

Fig. 10 Oil recovery as a function of time, displaced by SDS solution (blue, 100 μl min⁻¹), CO₂ gas (black, P_g = 22 psi), and CO₂ foam (red, R_g = 0.57 at P_g = 22 psi). Time axes differ: bottom axis for SDS, top axis for gas and foam. The horizontal gray lines correspond to the oil RF of the high permeable zone (≈ 70%) and complete recovery (100%).

Consistent with the observations in Fig. 6, the SDS solution sweeps the HP zone smoothly, forming a stable displacement front, and reaches a plateau at ≈70%, the high-permeable RF. However, it fails to displace oil from the LP zone, likely due to insufficient injection pressure to overcome the zone's capillary entry pressure. In contrast, both CO₂ gas and foam reach the HP zone plateau and proceed to displace oil from the LP zone—after ≈40 s and 15 s, respectively. The gas pressure of 22 psi in both cases is sufficient to overcome the capillary threshold in the LP region.

Furthermore, CO₂ foam achieves complete recovery in just 31 s, whereas gas alone requires over 150 s—roughly five times longer. This enhanced performance is attributed to two main factors: (1) foam suppresses viscous fingering, in contrast to pure gas (Fig. 11(a and b)), and (2) foam lamellae may block pores in the HP zone, raising the pressure drop in the porous medium and redirecting flow into the LP zone (Fig. 11(b, right)). In addition, while foam initially sweeps the HP zone, it is the liquid phase that ultimately penetrates and displaces oil from the LP region.

Fig. 11 Displacing fluid protrusions: (a) fingering of CO₂ gas into the low-permeable zone at P_g = 22 psi. (b) Oil displacement by CO₂ foam showing a more stable horizontal front (left), along with vertical protrusions resulting from the foam-blocked high-permeable zone into the low-permeable zone (right). All scale bars represent 500 μm.

4 Conclusions

We experimentally investigate pore-scale CO₂ foam dynamics and oil displacement behavior in a microfluidic heterogeneous porous medium, designed to mimic dual-permeability reservoir structures. The microfluidic device was fabricated from PDMS and featured two parallel homogeneous zones—high- and low-permeable regions—characterized by a permeability ratio of 5.8 and porosities of 45% and 23%, respectively. Oil-displacement experiments were conducted under ambient conditions using ISOPAR M oil (as the displaced fluid) and three displacing agents: SDS aqueous solution, CO₂ gas, and CO₂ foam.

Stable foam bubbles were generated via a flow-focusing geometry with an orifice width W_orifice = 60 μm. Bubble sizes ranged from 30 μm to 270 μm and were sensitive to both gas pressure (P_g = 22 and 26 psi) and the liquid flow rate (Q_l = 25–300 μL min⁻¹). Larger foam bubbles formed at lower Q_l and higher P_g. The bubble diameter normalized by the orifice width (D̄_b/W_orifice) followed a power-law trend against the gas-to-liquid flow rate ratio (Q_g/Q_l), with an exponent of 0.6 under constant P_g, aligning with literature models of foam formation in microfluidic systems.

Foam's structure and velocity were characterized by analyzing the gas areal fraction (F_g), foam texture (n_f), and velocity when foam flowed in the microfluidic porous medium at a constant gas pressure of 22 psi and varied liquid flow rates of 10 to 400 μL min⁻¹, producing gas ratios (R_g) of 0.25–0.93. As R_g increased, F_g rose linearly from 0.50 to 0.94, while n_f decreased, indicating coarser foam. Subdivision of F_g into trapped and mobile components revealed that higher R_g resulted in significantly greater gas trapping, with the trapped gas fraction (F_g,tb) increasing from 0.10 to 0.90. Spatial analysis showed that trapped bubbles exclusively occupied the low-permeable zone, while both trapped and mobile bubbles were present in the high-permeable region, which served as the primary flow pathway. Despite the changes in the foam structure, the average diameter of moving bubbles (D̄_mb) remained relatively stable (0.074–0.091 mm) for R_g ≤ 0.76. Their velocities (〈V̄_mb〉) ranged from 75 to 95 mm s⁻¹ without a clear trend. However, the overall foam velocity (V̄_f), estimated using the relative proportions of mobile and trapped gas, decreased with R_g, reflecting the growing impact of gas entrapment on flow resistance.

Oil displacement experiments provided a comparative evaluation of recovery performance across different displacing fluids. Both SDS solution and CO₂ foam formed stable displacement fronts in the high-permeable zone, whereas CO₂ gas led to pronounced viscous fingering and very early breakthrough. For SDS solution, increasing Q_l from 5 to 100 μL min⁻¹ reduced T_bt from 6.6 to 1.0 s, but also decreased RF_bt from 51.7% to 35.9%. CO₂ gas yielded T_bt values between 0.02 and 0.18 s as P_g increased from 15 to 28 psi; however, RF_bt remained consistently low (12–18%) due to rapid fingering and inefficient sweep. Foam flooding, by contrast, demonstrated delayed breakthrough and improved recovery: T_bt ranged from 1.52 to 0.54 s and RF_bt from 69.1% to 44.2% as R_g increased from 0.57 to 0.98.

Although higher R_g improved bubble generation, it also promoted foam destabilization in contact with oil, likely due to the reduced surfactant availability and shear-thinning rheology. Nevertheless, CO₂ foam injection significantly outperformed CO₂ gas in recovery dynamics. Complete oil recovery (100%) was achieved within 31 s using CO₂ foam, compared to over 150 s with CO₂ gas—demonstrating the foam's ability to overcome capillary entry pressures in the low-permeable zone more efficiently. This improvement was attributed to foam lamellae that blocked high-permeability pathways, increased resistance, and redirected flow vertically into unswept low-permeable regions. Such mechanisms were absent in gas-only injection, which relied solely on pressure-driven flow and suffered from bypassing.

These findings highlight that CO₂ foam's powerful potential for enhancing mobility control, sweep efficiency, and displacement depth in heterogeneous porous media. A promising recovery strategy could involve a sequential injection: first deploying SDS solution to establish a stable front in high-permeable channels, followed by CO₂ foam to divert flow into low-permeable areas. Operational parameters—including surfactant concentration, gas/liquid flow rates, foam quality (R_g), and injection pressure—must be optimized based on desired outcomes, technical constraints, and reservoir conditions.

Future work should explore the effect of reservoir-relevant factors such as elevated temperature and pressure, brine salinity, crude oil properties, and surface mineralogy on foam stability and flow behavior. These factors are essential for scaling up microfluidic insights to practical EOR, aquifer remediation, and carbon sequestration applications.

Author contributions

Nikoo Moradpour: writing original draft, visualization, software, methodology, investigation, formal analysis. Peichun Amy Tsai: conceptualization, methodology, writing, review & editing, validation, visualization, supervision, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information: two videos are provided in the supplementary information (SI): Video: Fig2_Foam bubble characterization: bubble regimes generated by the flow-focusing geometry at varied gas pressures (P_g) and liquid flow rates (Q_l). This is a representative video for the results in Fig. 2. Video: Fig4_Effect of gas ratio on bubble velocity: (left) trapped and moving bubbles in a representative porous medium at R_g = 0.25 and (right) moving bubbles (white) in the microfluidic porous medium at the gas ratios R_g = 0.25 (top), R_g = 0.63 (middle), and R_g = 0.93 (bottom). This is a representative video for the results in Fig. 4. See DOI: https://doi.org/10.1039/D5LC00544B.

Data analysis was conducted using ImageJ Fiji and MATLAB.

The data supporting the findings of this study are available upon reasonable request.

5 Appendix

5.1 Permeability of the porous medium

The average permeability, k_avg, is estimated by injecting pure water into the microfluidic chip and recording the pressure drop, ΔP, across the porous medium. The water injection flow rate is varied in the range of Q = 25–600 μl min⁻¹. Fig. 12(a) shows the result of Q/A versus ΔP/μL. L and A are the length and cross-sectional area of the porous medium, and μ is the water viscosity. According to Darcy's law, the slope of the linear fit in this plot represents the permeability of the porous medium, k_avg = 6.07 × 10⁻¹² m² or 6.07 Darcy.

Fig. 12 Permeability estimation of the porous medium using Darcy's law. (a) The flow rate per cross-sectional area (Q/A) is plotted against the measured pressure drop divided by the product of medium length and water viscosity, ΔP/Lμ. The slope yields the permeability. (b) (Left) the top view schematic of the circular pillar configuration (gray) with the pore area (green) between them. (Right) the 3D schematics of the pore volume and wetted surface area used to determine the hydraulic radius of the pore. (c) Schematic of fluid flow into two parallel layers with low and high permeability values. The pressure drop, ΔP, length, L, and the viscosity of the flowing fluid, μ, in the layers are equal while their cross sectional area, A, and liquid flow rate, Q, are different.

We used the Kozeny–Carman equation to find the ratio of the permeability of the low- and high- permeable layers using their hydraulic diameter, D_h, and porosity, ϕ.^70,98D_h is defined as the pore volume/wetted area, indicated in Fig. 12(b).

In the above equation, high and lower subscripts denote the properties for the high-permeable and low-permeable zones.

Given that the average permeability of two parallel layers (see Fig. 12(c)) can be determined using Darcy's law, the relationship is expressed as follows:

Q_total = Q_high + Q_low, (5)

Knowing that A_low = A_high = A_total/2, eqn (6) is reduced to:Using k_avg, eqn (4), and eqn (7), k_low and k_high are calculated to be 1.8 and 10.1 Darcy, respectively.

5.2 Darcy's law

Fluid flow through porous media can be described using Darcy's law:⁹⁹where Q = fluid flow rate, m³ s⁻¹, K = permeability of the medium, m², A = open cross-sectional area, m², ΔP = pressure drop across the medium, Pa, L = length of the medium, m, and μ = viscosity of the fluid, Pa·s.

To estimate ΔP for oil displacement by SDS solution, Q_l = 25 μl min⁻¹, A = 6.3 × 10⁻⁷ m², L = 1.4 m, K = 6.07 × 10⁻¹² m², and μ is very close to 1 cP at an SDS concentration of 0.5 wt%.^91,92 Therefore, ΔP = 0.23 psi.

5.3 Contact-angle measurement results

Droplets of SDS solution (0.5 wt%) and ISOPAR M oil were placed on clean surfaces of bare glass slides and PDMS slabs. A high-resolution camera captured snapshots of the liquid droplets on the solid substrates. Fiji ImageJ software was used to find the left and right contact angles (inner side as shown in Fig. 13). The values were averaged to determine the reported ones.

Fig. 13 The contact angles of SDS solution and ISOPAR M oil on bare glass and PDMS. The contact angles are measured by the sessile drop technique. The right and left contact angles are averaged and the mean values are listed in the table.

5.4 IFT measurement results

Interfacial tension (IFT) between the fluids was obtained using the pendant drop technique. For gas–liquid IFT measurements, a transparent cuvette was filled with the liquid phase and the gas bubbles were released using a syringe and a U-shaped needle. Oil droplets were released in the SDS solution similar to that above. The captured images were processed for drop shape analysis^71,72 and the attained IFT values are listed in Fig. 14.

Fig. 14 The interfacial tension, σ, measured between different fluids (air, CO₂, SDS solution, and ISOPAR M oil) using the pendant drop method.

5.5 Oil displacement by CO₂ foam

The raw and binarized temporal sequence of oil displacement by CO₂ foam is displayed in Fig. 15. The foam is injected at P_g = 22 psi with R_g = 0.87. This figure represents the same experiment shown in Fig. 6(c).

Fig. 15 Oil displacement dynamics: image sequences of oil displacement by CO₂ foam by SDS solution at R_g = 0.87 and P_g = 22 psi. Displacement progress is from left to right; the upper row represents raw images and the lower row illustrates the processed images where colored regions represent displacing fluids and white regions include oil and pillars. “HP” and “LP” refer to high- and low-permeable zones. The magnification is 2.5×. The scale bar represent 1000 μm.

Acknowledgements

We gratefully acknowledge the support from the Canada First Research Excellence Fund (CFREF), the Future Energy System (FES T02-P05 CCUS projects) at the University of Alberta, and the Canada Foundation for Innovation (CFI 34546). P. A. T. holds a Canada Research Chair (CRC) in Fluids and Interfaces and gratefully acknowledges funding from the Natural Sciences and Engineering Research Council of Canada (NSERC), in particular the NSERC Canada Research Chairs Program (CRC 233147) and Discovery Grant (RGPIN-2020-05511).

Notes and references

1. D. L. Weaire and S. Hutzler, The physics of foams, Oxford University Press, 1999.
2. W. Drenckhan and S. Hutzler, Adv. Colloid Interface Sci., 2015, 224, 1–16.
3. L. L. Schramm, Emulsions, foams, suspensions, and aerosols: microscience and applications, John Wiley & Sons, 2014.
4. B. Dollet and C. Raufaste, C. R. Phys., 2014, 15, 731–747.
5. C. A. Conn, K. Ma, G. J. Hirasaki and S. L. Biswal, Lab Chip, 2014, 14, 3968–3977.
6. K. Ma, R. Liontas, C. A. Conn, G. J. Hirasaki and S. L. Biswal, Soft Matter, 2012, 8, 10669–10675.
7. P. Nguyen, H. Fadaei and D. Sinton, Energy Fuels, 2014, 28, 6221–6227.
8. A. Karthick, B. Roy and P. Chattopadhyay, J. Environ. Manage., 2019, 243, 187–205.
9. S. S. Datta, I. Battiato, M. A. Fernø, R. Juanes, S. Parsa, V. Prigiobbe, E. Santanach-Carreras, W. Song, S. L. Biswal and D. Sinton, Lab Chip, 2023, 23, 1358–1375.
10. W. R. Rossen, R. Farajzadeh, G. J. Hirasaki and M. Amirmoshiri, SPE Improved Oil Recovery Conference?, 2022, p. D021S014R001.
11. W. Wanniarachchi, P. Ranjith, M. Perera, T. Rathnaweera, D. Zhang and C. Zhang, Eng. Fract. Mech., 2018, 194, 117–135.
12. T. Lu, Z. Li and L. Du, J. Chem. Eng., 2024, 481, 148575.
13. K. R. Chaturvedi, S. Bajpai, J. Trivedi and T. Sharma, Energy Fuels, 2022, 36, 5036–5046.
14. P. Liu, X. Zhang, Y. Wu and X. Li, J. Pet. Sci. Eng., 2017, 150, 208–216.
15. Z.-x. Liu, Y. Liang, Q. Wang, Y.-j. Guo, M. Gao, Z.-b. Wang and W.-l. Liu, J. Pet. Sci. Eng., 2020, 193, 107449.
16. C. McGlade, G. Sondak and M. Han, Whatever happened to enhanced oil recovery?, International energy agency technical report, 2018.
17. N. Moradpour, J. Yang and P. A. Tsai, Curr. Opin. Colloid Interface Sci., 2024, 101845.
18. T. Muther, H. A. Qureshi, F. I. Syed, H. Aziz, A. Siyal, A. K. Dahaghi and S. Negahban, J. Pet. Explor. Prod. Technol., 2022, 12, 1463–1488.
19. R. Santos, W. Loh, A. Bannwart and O. Trevisan, Braz. J. Chem. Eng., 2014, 31, 571–590.
20. N. Moradpour, R. Azadi and P. A. Tsai, Colloids Surf., A, 2025, 705, 135533.
21. H. Manikantan and T. M. Squires, J. Fluid Mech., 2020, 892, P1.
22. S. Xiao, Y. Zeng, E. D. Vavra, P. He, M. Puerto, G. J. Hirasaki and S. L. Biswal, Langmuir, 2018, 34, 739–749.
23. R. Phukan, S. B. Gogoi and P. Tiwari, Colloids Surf., A, 2020, 597, 124799.
24. K. R. Chaturvedi and T. Sharma, Chem. Eng. Sci., 2021, 235, 116484.
25. V. A. Lifton, Lab Chip, 2016, 16, 1777–1796.
26. A. Anbari, H. T. Chien, S. S. Datta, W. Deng, D. A. Weitz and J. Fan, Small, 2018, 14, 1–15.
27. K. Osei-Bonsu, N. Shokri and P. Grassia, J. Colloid Interface Sci., 2016, 462, 288–296.
28. M. J. Shojaei, K. Osei-Bonsu, S. Richman, P. Grassia and N. Shokri, Ind. Eng. Chem. Res., 2018, 58, 1068–1074.
29. B. Géraud, S. A. Jones, I. Cantat, B. Dollet and Y. Méheust, Water Resour. Res., 2016, 52, 773–790.
30. M. Simjoo, Y. Dong, A. Andrianov, M. Talanana and P. Zitha, Ind. Eng. Chem. Res., 2013, 52, 6221–6233.
31. G. Shi, K. Tang, F. Wang, Q. Luo, L. Bai, P. Sun and D. Zhu, Energy Fuels, 2020, 35, 465–472.
32. S. Li, Q. Wang, K. Zhang and Z. Li, Fuel, 2020, 263, 116648.
33. D. Sinton, Lab Chip, 2014, 14, 3127–3134.
34. J. Yang and P. A. Tsai, Chem. Eng. Sci., 2024, 300, 120543.
35. T. H. M. Ho and P. A. Tsai, Lab Chip, 2020, 20(20), 3806–3814.
36. T. H. M. Ho, D. Sameoto and P. A. Tsai, Chem. Eng. Res. Des., 2021, 174, 116–126.
37. J. Yang, M. Saadat, I. Azizov, M. Dudek, G. Øye and P. A. Tsai, Lab Chip, 2022, 22(24), 4974–4983.
38. J. Zhao, F. Torabi and J. Yang, Fuel, 2021, 287, 119443.
39. N. Yekeen, M. A. Manan, A. K. Idris, A. M. Samin and A. R. Risal, J. Pet. Sci. Eng., 2017, 159, 115–134.
40. M. Lv, Z. Liu, C. Ji, L. Jia and Y. Jiang, Ind. Eng. Chem. Res., 2018, 57, 15172–15180.
41. W. Yu, X. Zhou and M. Y. Kanj, Langmuir, 2022, 38, 2895–2905.
42. A. Gizzatov, S. Pierobon, Z. AlYousef, G. Jian, X. Fan, A. Abedini and A. I. Abdel-Fattah, Sci. Rep., 2021, 11, 3360.
43. Z. Xu, Z. Li, S. Cui, B. Li, D. Chen, Q. Zhang, L. Zheng and M. M. Husein, J. Pet. Sci. Eng., 2022, 211, 110170.
44. Q. Lv, R. Zheng, T. Zhou, X. Guo, W. Wang, J. Li and Z. Liu, Fuel, 2022, 330, 125533.
45. G. Jian, A. Gizzatov, M. Kawelah, Z. AlYousef and A. I. Abdel-Fattah, Appl. Energy, 2021, 292, 116815.
46. B. Benali, T. L. Føyen, Z. P. Alcorn, M. Haugen, J. Gauteplass, A. R. Kovscek and M. A. Fernø, Int. J. Greenhouse Gas Control, 2022, 114, 103607.
47. L. F. Lopes, J. M. Façanha, L. Maqueira, F. R. Ribeiro and A. Pérez-Gramatges, J. Pet. Sci. Eng., 2021, 207, 109141.
48. W. Yan, C. A. Miller and G. J. Hirasaki, Colloids Surf., A, 2006, 282, 348–359.
49. S. Jones, N. Getrouw and S. Vincent-Bonnieu, Soft Matter, 2018, 14, 3497–3503.
50. A. Kovscek, T. Patzek and C. Radke, Chem. Eng. Sci., 1995, 50, 3783–3799.
51. N. Nazari and A. R. Kovscek, Lab Chip, 2022, 22, 3489–3498.
52. W. R. Rossen, Colloids Surf., A, 2003, 225, 1–24.
53. J. Alvarez, H. Rivas and W. Rossen, SPE J., 2001, 6, 325–333.
54. Y. Wang, M. Puerto, C. Bai, K. Ma, K. Mateen, G. J. Hirasaki and S. L. Biswal, Ind. Eng. Chem. Res., 2025, 64, 2995–3003.
55. X. Zheng, N. Mahabadi, T. S. Yun and J. Jang, J. Geophys. Res.: Solid Earth, 2017, 122, 1634–1647.
56. Q. Sun, Z. Li, S. Li, L. Jiang, J. Wang and P. Wang, Energy Fuels, 2014, 28, 2384–2394.
57. F. Guo, S. A. Aryana, Y. Wang, J. F. McLaughlin and K. Coddington, Int. J. Greenhouse Gas Control, 2019, 87, 134–141.
58. F. Guo and S. A. Aryana, J. Pet. Sci. Eng., 2018, 164, 155–163.
59. D. Li, G. Xin and S. Ren, ACS Omega, 2022, 7, 36503–36509.
60. S. Jones, N. Getrouw and S. Vincent-Bonnieu, Soft Matter, 2018, 14, 3490–3496.
61. A. E. Peksa, K.-H. A. Wolf and P. L. Zitha, Mar. Pet. Geol., 2015, 67, 701–719.
62. B. A. Abuamarah and B. S. Nabawy, J. Nat. Gas Sci. Eng., 2021, 88, 103807.
63. P.-R. Thomson, M. Jefferd, B. L. Clark, D. Chiarella, T. M. Mitchell and S. Hier-Majumder, Mar. Pet. Geol., 2020, 122, 104614.
64. P.-R. Thomson, R. Ellis, D. Chiarella and S. Hier-Majumder, Front. Earth Sci., 2020, 8, 246.
65. M. Kashif, Y. Cao, G. Yuan, M. Asif, K. Javed, J. N. Mendez, D. Khan and L. Miruo, Pet. Sci., 2019, 16, 981–1000.
66. S. Ma, G. Jin, R. Sy and H. Kesserwan, SPE Kingdom of Saudi Arabia Annual Technical Symposium and Exhibition, 2018, p. SPE-192220.
67. R. Kareem, P. Cubillas, J. Gluyas, L. Bowen, S. Hillier and H. C. Greenwell, J. Pet. Sci. Eng., 2017, 149, 436–455.
68. M. Wetzel, T. Kempka and M. Kühn, Minerals, 2021, 11, 151.
69. H. Safari, B. J. Balcom and A. Afrough, Comput. Geosci., 2021, 156, 104895.
70. E. J. Peters, Advanced Petrophysics: Geology, porosity, absolute permeability, heterogeneity, and geostatistics, Greenleaf Book Group, 2012, vol. 1.
71. J. D. Berry, M. J. Neeson, R. R. Dagastine, D. Y. Chan and R. F. Tabor, J. Colloid Interface Sci., 2015, 454, 226–237.
72. A. Esfandiarian, A. Maghsoudian, A. Davarpanah, Y. Tamsilian and S. Kord, J. Pet. Sci. Eng., 2022, 215, 110607.
73. N. D. Denkov, K. G. Marinova and S. S. Tcholakova, Adv. Colloid Interface Sci., 2014, 206, 57–67.
74. R. Aveyard, B. Binks, P. Fletcher, T. Peck and C. Rutherford, Adv. Colloid Interface Sci., 1994, 48, 93–120.
75. J. Broeke, J. M. M. Pérez and J. Pascau, Image processing with ImageJ, Packt Publishing Ltd, 2015.
76. T. Ferreira and W. Rasband, ImageJ user guide, ImageJ/Fiji, 2012.
77. P. Garstecki, I. Gitlin, W. DiLuzio, G. M. Whitesides, E. Kumacheva and H. A. Stone, Appl. Phys. Lett., 2004, 85, 2649–2651.
78. N. Dietrich, S. Poncin, N. Midoux and H. Z. Li, Langmuir, 2008, 24, 13904–13911.
79. T. Fu, D. Funfschilling, Y. Ma and H. Z. Li, Microfluid. Nanofluid., 2010, 8, 467–475.
80. S. Cleve, A. Lassus, C. Diddens, B. Van Elburg, E. Gaud, S. Cherkaoui, M. Versluis, T. Segers and G. Lajoinie, J. Fluid Mech., 2023, 972, A27.
81. B. Géraud, Y. Méheust, I. Cantat and B. Dollet, Phys. Rev. Lett., 2017, 118, 098003.
82. Z. Huang, M. Su, Q. Yang, Z. Li, S. Chen, Y. Li, X. Zhou, F. Li and Y. Song, Nat. Commun., 2017, 8, 1–9.
83. Y. Zhang, W. Geng, M. Chen, X. Xu, L. Jiang and Y. Song, Energy Fuels, 2023, 37, 19695–19703.
84. H. M. Polat, F. M. Coelho, T. J. Vlugt, L. F. Mercier Franco, I. N. Tsimpanogiannis and O. A. Moultos, J. Chem. Eng. Data, 2024, 69, 3296–3329.
85. J. J. Carroll, J. D. Slupsky and A. E. Mather, J. Phys. Chem. Ref. Data, 1991, 20, 1201–1209.
86. L. W. Diamond and N. N. Akinfiev, Fluid Phase Equilib., 2003, 208, 265–290.
87. S. Liu, X. Dou, Q. Zeng and J. Liu, J. Pet. Sci. Eng., 2021, 196, 107635.
88. S. Ansari and D. S. Nobes, Fluids, 2022, 7, 333.
89. F. Mo, Z. Du, X. Peng, Y. Tang and H. Sun, J. Nat. Gas Sci. Eng., 2017, 44, 314–327.
90. D. D. Huang and M. M. Honarpour, J. Pet. Sci. Eng., 1998, 19, 103–117.
91. T. P. Niraula, S. K. Shah, S. K. Chatterjee and A. Bhattarai, Karbala Int. J. Mod. Sci., 2018, 4, 26–34.
92. Y. Chai, X. Li, J. Geng, J. Pan, Y. Huang and D. Jing, Colloid Polym. Sci., 2019, 297, 1025–1035.
93. H. B. Todd and J. G. Evans, SPE Rocky Mountain Petroleum Technology Conference/Low-Permeability Reservoirs Symposium, 2016, p. SPE-180270.
94. L. Xue, P. Liu and Y. Zhang, Water, 2023, 15, 1342.
95. L. Zhao, T. Liu, F. Lu, J. Sun and Y. Yao, ACS Omega, 2024, 9, 28505–28519.
96. W. Osterloh and M. Jante Jr, SPE Improved Oil Recovery Conference?, 1992, p. SPE-24179.
97. R. Gajbhiye and S. Kam, Chem. Eng. Sci., 2011, 66, 1536–1549.
98. P. C. Carman, Chem. Eng. Res. Des., 1997, 75, S32–S48.
99. T. Ahmed, Reservoir engineering handbook, Gulf professional publishing, 2018.

---