Fluid entrapment during forced imbibition in a multidepth microfluidic chip with complex porous geometry

Abstract

Understanding and controlling fluid entrapment during forced imbibition in porous media is crucial for many natural and industrial applications. However, the microscale physics and macroscopic consequences of fluid entrapment in these geometric-confined porous media remain poorly understood. Here, we introduce a novel multidepth microfluidic chip, which can mitigate the depth confinement of traditional two-dimensional (2-D) microfluidic chips and mimic the wide pore size distribution as natural-occurring three-dimensional (3-D) porous media. Based on microfluidic experiments and direct numerical simulations, we observe the fluid-entrapment scenarios and elucidate the underlying complex interaction between geometric confinement, capillary number and wettability. Increasing depth variation can promote fluid entrapment, whereas increasing capillary number and contact angle yield the opposite effect, which seemingly contradicts conventional expectations in traditional 2-D microfluidic chips. The fluid-entrapment scenario in depth-variable microfluidic chips stems from microscopic interfacial phenomena, classified as snap-off and bypass events. We provide theoretical analyses of these pore-scale events and validate corresponding phase diagrams numerically. It is shown that increasing depth variation triggers snap-off and bypass events. Conversely, a higher capillary number suppresses snap-off events under strong imbibition, and an increased contact angle inhibits bypass events under imbibition. These macroscopic imbibition patterns in microfluidic porous media can be linked with these pore-scale events by improved dynamic pore-network models. Our findings bridge the understanding of forced imbibition between 2-D and 3-D porous media and provide design principles for newly engineered porous media with respect to their desired imbibition behaviours.

1. Introduction

Immiscible two-phase displacement in permeable media has drawn extensive research attention across various fields of science and engineering, such as microfluidic logic control (Cybulski, Garstecki & Grzybowski 2019), solvents separation in chromatography columns (Shalliker et al. 2007), formation of liquid-infused materials (Yeganeh et al. 2022), film wetting in nanofabrication (Anand et al. 2021), geological carbon sequestration (Szulczewski et al. 2012) and hydrocarbon recovery (Blunt 2017). However, complex interfacial behaviour and multiscale variability between fluids and porous structures pose great challenges to understanding and controlling multiphase flow dynamics and fluid entrapment. For instance, in geological carbon sequestration, viscous instability during CO₂ rising through brine-filled pore spaces and capillary instability during water film spreading in complex porous structures can significantly increase trapping efficiency by increasing the interfacial area (Huppert & Neufeld 2014). Unlike the drainage process where the fluid entrapment patterns in porous media with fractal fluid–fluid fronts are often observed, the imbibition pattern is typically expected to be stable and compact due to the favoured pore-scale cooperative pore-filling mechanism (Stokes et al. 1986; Cieplak & Robbins 1988, 1990; Holtzman & Segre 2015; Jung et al. 2016; Singh et al. 2019). However, the recent discoveries of incomplete imbibition processes in various structures make the above pore-scale imbibition results more subtle (Odier et al. 2017; Lei et al. 2022a; Wang et al. 2022). Clearly, there is a need for a thorough examination of imbibition mechanisms in diverse porous structures.

During the forced imbibition process, the macroscopic displacement patterns are governed by the competition between the capillary force and the viscous force, which is further complicated by flow rate, wettability and complex geometry (Rabbani et al. 2018; Hu et al. 2019; Mehmani et al. 2019; Primkulov et al. 2019, 2021; Singh et al. 2019; Lei et al. 2022a; Wang et al. 2022; Lei, Gong & Wang 2023a). Flow regimes in complete or incomplete displacements are controlled not only by the capillary number (the ratio of viscous to capillary force) (Stokes et al. 1986; Cieplak & Robbins 1988, 1990; Zhao, Macminn & Juanes 2016; Primkulov et al. 2021), but also by the local capillary pressure variation (Craster & Matar 2009; Hu et al. 2017; Gu, Liu & Wu 2021; Wang, Mehmani & Xu 2021). The capillary pressure$\; {P_c} = 2\gamma \,\textrm{cos}\,\theta /{R_h}$ is determined by interfacial tension (IFT) $\gamma $, contact angle $\theta $ and hydraulic radius of pores ${R_h}$. While most focus has been placed on the effects of fluid properties or injection conditions, the influence of 3-D geometrical features of porous structures on multiphase flow received less attention. To control and optimize fluid entrapment in applications, it is of crucial importance to understand the interplay between capillary number, contact angle and hydraulic diameter on forced imbibition in porous media with diverse geometric confinements.

Various natural and industrial porous media present different geometric constraints, as appear in the geological formation, vascular networks and membranes. By affecting pore-scale interfacial behaviour, such as the competition between the main meniscus and water film interface, geometric confinements of porous media lead to different fluid entrapments in complex porous structures (Lei et al. 2023c). Compared with the difficulty of visualization and quantification of multiphase flow in real-world 3-D porous media (Berg et al. 2013), microfluidic chips with complex porous geometry provide a convenient and flexible tool to visualize and quantify both microscopic flow dynamics and macroscopic fluid configuration (Sinton 2014). These microfluidic chips as ‘2-D models’ can realize size variation in horizontal dimensions but constant depth in the vertical dimension (Lei, Lu & Wang 2023b). However, there exist strong geometric confinements in these 2-D microfluidic chips with complex porous geometry that ${\lim _{W/H \gg 1}}{R_h} = {\lim _{W/H \gg 1}}(W/(W + H))H \approx H$, where W and H are the pore width and depth, respectively. Consequently, the hydraulic radius of designed larger pores will be confined by depth into the diameter range of 0∼H. This depth confinement will always result in a homogeneous porous media with a narrow 3-D pore size distribution whatever the heterogeneous design of the 2-D porous structure. The capillary pressure-saturation functions of porous media can be influenced by the pore size distribution directly (Rücker et al. 2020). Considering hydraulic-radius-dependent capillary pressure, the results based on 2-D microfluidic porous media may be not sufficient for fully understanding the microscopic capillary physics and macroscopic flow patterns during imbibition in 3-D porous media.

A few attempts have been made to add 3-D geometric features into microfluidic chips, such as 3-D printing (Ju, Gong & Zheng 2022) and glass-based isotropic hydrofluoric acid (HF) etching (Xu et al. 2017). Unfortunately, as yet, there has been no complete and rigorous fabrication technique reported for a microfluidic chip to mitigate this depth-confined capillarity limitation. The space resolution (tens of microns) and limited printing material (fusible and transparent) of current 3-D printing technology make it impossible to build a 3-D microfluidic chip with complex porous geometry (Yazdi et al. 2016). Glass-based isotropic HF etching can achieve a shallower throat if two designed neighbouring pores are close enough (Xu et al. 2017), while this special design is limited by the channel aspect ratio in wet-etched glass. Moreover, it is difficult to fabricate microfluidic chips with complex porous geometry on the glass due to the dense, complex and small patterns (Lei et al. 2023b).

In this paper, we combine advanced microfluidic experiments, numerical simulations and theoretical analysis to elucidate the complex interplay between capillary number, wettability and geometric confinement on multiphase flow dynamics and fluid entrapment. We first present a novel method for designing and fabricating multidepth microfluidic chips (§ 2.1) and the corresponding microfluidic platform (§ 2.2). Then, we introduce the fundamentals of direct numerical simulations with OpenFOAM (§ 2.3) and an improved dynamic pore-network model (§ 2.4). By combining microfluidic experiments and direct numerical simulations (§§ 3.1 and 3.2), we determine the underlying flow mechanism and purify the complex displacement modes as the basic pore-scale events. Theoretical models of these pore-scale events, such as snap-off and bypass, were proposed and numerically verified (§ 3.3). Finally, by inputting these pore-scale events into the improved pore-network model, these microscopic pore-scale events were linked directly to the macroscopic displacement patterns (§ 3.4). Our findings highlight the interplay of depth confinement, capillary number and contact angle on both microscopic flow dynamics and macroscopic forced imbibition patterns. These insights not only bridge the understanding of fluid entrapment in various 2-D and 3-D porous media but also open the door for the potential of 3-D microfluidic design for striking macroscopic imbibition patterns.

2. Materials and methods

2.1. Design and fabrication of multidepth microfluidic chips

In this study, the sequential photolithography and multietching processes are combined for fabricating silicon-based multidepth microfluidic chips. Among these processing steps, mask design is the crucial step. Figure 1 shows the whole process from the 2-D image generation to the area segmentation for labelling depth information. Here, a high-resolution 2-D porous structure, aiming for a target pore size distribution resembling natural 3-D porous media, is generated using a random generation algorithm known as the quartet structure generation set (Wang et al. 2007; Lei et al. 2020). To simplify the microfabrication processes, the target depth of each etched pattern is configured to the same value. Here, this value can be chosen as the characteristic diameter of the desired pore size distribution (40 μm). Therefore, the depth value of pores in the microfluidic chips can be 40, 80 or 120 μm.

Figure 1. Labelling different depth regions according to the width distribution of the 2-D geometry, which is regenerated based on the statistic information of natural porous media. (a) The 2-D diameter distribution calculated by the maximal ball algorithm. (b) The basic rule of dividing the depth regions according to the width range of the 2-D geometry. (c) Corresponding designed depth distribution on the 2-D geometry to achieve the target pore size distribution in multidepth microfluidic chips. (d) Different etching layers in sequential photolithography and multietching processes.

We utilized the maximal ball algorithm (Blunt 2017) to identify the 2-D width distribution of the porous structure, as shown in figure 1(a). Statistical information about the width distributions (figure 1b) can be acquired by calculating the ratios of these pore areas. Considering the impact of depth confinement on the hydraulic diameter’s distribution, we divided the aforementioned 2-D diameter-labelled regions into three categories: the first-depth region, the second-depth region and the third-depth region, as shown in figure 1(b,c). To achieve depth variations in pores and ensure mismatch tolerance during the fabrication, we designed the etching range to shrink to the inside of the pores with the etching times increase. As a result, the first mask for the first etching is the original 2-D image containing the all-depth region, the second mask for the second etching is extracted from the first mask that contains the second-depth and third-depth regions, the third mask for the third etching is extracted from the second mask that only contains the third-depth region. Subsequently, the corresponding 2-D geometry can be categorized into different depth regions (figure 1c), enabling the fabrication of depth-variable microfluidic chips through sequential photolithography and multietching methods (figure 1d).

It is worth noting that directly dividing the width distribution into three regions according to 2-D diameter values will result in discontinuities in the 3-D hydraulic diameter distribution due to $2{R_{{h_i}}} = W{H_i}/(W + {H_i})\; \ne \; 2{R_{{h_j}}} = W{H_j}/(W + {H_j})$, where ${H_i}$ and ${H_j}$ are different depth values, i and j are different depth labels. To address this issue and maintain a continuous pore size distribution, it is essential to include an overlap width range between two different depth regions, where Area II was designed to make up the gap caused by the first-depth region and the second-depth region that $1/W{2_{min}} + 1/2H = 1/W{1_{max}} + 1/H$ and Area IV was designed to make up the gap caused by the second-depth region and third-depth region that $1/W{3_{min}} + 1/3H = 1/W{2_{max}} + 1/2H$, where W1_max, W2_max and W3_max represent the largest width for first-depth, second-depth and third-depth regions, respectively; W2_min and W3_min represent the smallest width for second-depth and third-depth regions, respectively, as shown in figure 1(b). In the overlap width regions II and IV, we randomly select 50 % of the pore area labelled with shallower depth and the other pore area labelled with deeper depth.

Based on the obtained pattern, the porous structures were created from a silicon wafer using microfabrication techniques, including multiple uses of standard photography, and inductively coupled plasma-deep reactive ion etching. Figure 2(a i–iii) shows the corresponding sequential photolithography and multietching processes. With an increase in etching cycles, the geometric confinement in the depth direction is gradually weakened, which is evident in pore-network features that the larger 3-D pores are increased clearly (figure 2b i–iii). To get equal wettability on all surfaces within the microfluidic chip, all etched silicon wafers underwent a one-hour heating process at 1000 °C in the presence of oxygen to create a uniform 100 nm silicon dioxide layer via thermal oxidation (Chomsurin & Werth 2003). Finally, a Schott Borofloat 270 glass wafer with inlet and outlet ports was anodically bonded to the above thermally oxidized silicon. Figure 2(c i–iii) present the final morphology of uniform-depth and depth-variable microfluidic chips observed by the stereo microscope.

Figure 2. Sequence photolithography and multietching processes for depth-variable microfluidic chips with complex porous geometry. (a) Photolithography regions for first, second and third etching processes, where black, blue and red pixels represent pores etched sequentially. (b) Different pore-network features are generated by the SNOW algorithm (Gostick 2017) based on different microfluidic chips with uniform depth, two different depths and three different depths. Grey colour areas are small pores, blue colour areas are middle pores and red colour areas are large pores. (c) The corresponding images from (i) to (iii) are uniform-depth, two-different-depth and three-different-depth microfluidic chips.

The depth confinements of these microfluidic chips can be directly reflected by the 3-D pore size distribution, as shown in figure 3. Here, the depth variation factor is defined as $\alpha = {d_{max}}/{d_{min}}$, ${d_{max}}$ and ${d_{min}}$ are the maximum and minimum depth in the microfluidic chips. Smaller depth variation $\alpha $ will result in strong depth confinement (figure 3). Figure 3(b) presents the strong depth confinement ($\alpha = 1$) in pore size distribution, which is commonly observed in traditional 2-D microfluidic experiments. Our novel multidepth design in microfluidic chips ($\alpha = 2,3$) can release the depth confinement and even realize the wide 3-D pore size distribution as 3-D porous media (figure 3d).

Figure 3. Comparing pore size distributions in natural-occurring 3-D porous media and microfluidic chips. From (a) target pore size distribution to corresponding pore distributions of (b) traditional uniform-depth ($\alpha = 1$), (c) two-different-depth ($\alpha = 2$) and (d) three-different-depth ($\alpha = 3$) microfluidic chips. The 3-D microstructure of natural rock is regenerated by microcomputed tomography, and the local microstructure of microchips is scanned by the scanning electron microscope. The pore size distributions of natural rock and microfluidic chips are calculated by the maximal ball algorithm (Blunt 2017).

2.2. Experimental set-up and procedures

The experimental set-up used in this study has been introduced in our previous works (Lei et al. 2020, 2022b). The microfluidic chip was cleaned by acetone, alcohol, 4 vol. % hydrochloric acid and deionized water sequentially, and was dried by N₂, then fixed by a homemade holder under the fluorescence microscope (Nikon, SMZ18) with a high-speed camera (Nikon, DS-Ri2). The flow rate was controlled using a feedback-controlled pump (Fluigent, MFCS-EZ) with a precise flow meter (Fluigent, Flow-unit). The visualized displacement processes can be quantified by image processing and data analysis. During displacement experiments in the microfluidic porous media, we employed a spatial resolution of 1.92 μm pixel⁻¹ and a time resolution of 1 s frame⁻¹ to meet the requirements of a quasistatic experiment and analysis of the phenomena from the micron pore scale to centimetre porous media scale.

2.3. Direct numerical simulation methods

Direct numerical simulation within a computational fluid dynamics framework was conducted to investigate the impact of depth variation, capillary number and wettability on the multiphase flow dynamics and fluid entrapment during imbibition in porous media (Ferrari et al. 2015; Zheng et al. 2021). Gravity was neglected in this study because Bond number $Bd = \rho g{r^2}/\gamma \ll 1$, where $\rho $ is the fluid density, g is the gravitational acceleration, r is the characteristic pore size and $\gamma $ is fluid–fluid IFT. The imbibition processes can be captured by solving the Navier–Stokes (N–S) equations in the pore space and employing the volume-of-fluid (VOF) method to track the evolution of the fluid–fluid interface (Ferrari & Lunati 2013). The VOF method provides a simple and economical way to track the fluid–fluid interface, which can better ensure mass conservation. Generally, it is more flexible and efficient than other two-phase flow models in dealing with problems involving highly complex interface configurations. However, it is also worth noting that VOF simulations can be computationally expensive for problems involving complex geometry.

In our simulations, the direct solution of the governing equations was carried out using the open-source software OpenFOAM (interFoam solver). The wetting phase (water) and non-wetting phase (air or decane) can be treated as incompressible Newtonian fluids under a low flow rate. The flow system of experiments and simulations is under normal temperature and pressure conditions. The governing equations for the isothermal incompressible two-phase immiscible Newtonian fluids are

(2.1)\begin{gather}\frac{{\partial (\rho \boldsymbol{u})}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{uu}) ={-} \boldsymbol{\nabla }p + \boldsymbol{\nabla }\boldsymbol{\cdot }{\bf (}2\mu \boldsymbol{E}) + {\boldsymbol{f}_s},\end{gather}

(2.2)\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0,\end{gather}

where t designates time, the density $\rho $ and the dynamic viscosity $\mu $ are the function of volume fraction $\delta $ that $\rho (\delta ) = (1 - \delta ){\rho _{nw}} + \delta {\rho _w}$ and $\mu (\delta ) = (1 - \delta ){\mu _{nw}} + \delta {\mu _w}$. The volume fraction $\delta $ is introduced as an indicator function to represent the volume fraction of wetting fluid within flow cells, where $\delta = 1$ and $\delta = 0$ denote the wetting phase cell and non-wetting phase cell, respectively, while $0 < \delta < 1$ denotes the interfacial cell. The subscript nw stands for the non-wetting phase and the subscript w stands for the wetting phase. Here $\boldsymbol{E} = {\textstyle{1 \over 2}}(\boldsymbol{\nabla }\boldsymbol{u} + {\nabla ^T}\boldsymbol{u})$ is the rate-of-strain tensor, u is the vector of the fluid velocity, p is the pressure. Here ${\boldsymbol{f}_s}$ is the volumetric surface tension force and describes the effect of the Laplace pressure. These equations are combined with a transport equation for the volume fraction $\delta $,

(2.3)\begin{equation}\frac{{\partial \delta }}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\delta \boldsymbol{u}) = 0,\end{equation}

which can be regarded as a mass-conservation equation for the incompressible wetting fluid.

According to the continuum surface force model (Brackbill, Kothe & Zemach 1992), the volumetric surface tension force ${\boldsymbol{f}_s}$ can be calculated as

(2.4)\begin{equation}{\boldsymbol{f}_s} = \gamma \kappa \boldsymbol{\nabla }\delta ,\end{equation}

where $\kappa $ is the interface curvature. The interface curvature $\kappa $ is obtained as

(2.5)\begin{equation}\kappa ={-} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{n} ={-} \boldsymbol{\nabla }\boldsymbol{\cdot }\left( {\frac{{\boldsymbol{\nabla }\delta }}{{|{\boldsymbol{\nabla }\delta } |}}} \right).\end{equation}

When fluids flow through a porous medium, their access to space is restricted – they are confined within the pore space defined by the solid matrix surface. The fluid–wall interaction adopts no-slip boundary conditions. The wettability of the fluid–fluid–solid system is considered by imposing an additional constraint at the intersection between the interface and solid surface (Ferrari & Lunati 2013),

(2.6)\begin{equation}\boldsymbol{n}{|_{x \in wall}} = {\boldsymbol{n}_{\boldsymbol{s}}}\,\textrm{cos}\,\theta + {\boldsymbol{t}_{\boldsymbol{s}}}\,\textrm{sin}\,\theta ,\end{equation}

where $\theta $ is the static contact angle of the invading fluid, $\boldsymbol{n}$ is the normal cell vector of the fluid–fluid interface. Here ${\boldsymbol{n}_{\boldsymbol{s}}}$ is the normal cell pointing into the solid and ${\boldsymbol{t}_{\boldsymbol{s}}}$ is the cell tangent to solid points into the wetting phase.

The advection equation (2.3) is solved explicitly using the velocity from the previous time step. A nonlinear convective term is added in (2.3) to keep the interface sharp and create a shock that balances numerical diffusion. These yield

(2.7)\begin{equation}\frac{{\partial \delta }}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\delta \boldsymbol{u}) + \boldsymbol{\nabla }\boldsymbol{\cdot }[\delta (1 - \delta ){\boldsymbol{u}_r}] = 0,\end{equation}

where ${\boldsymbol{u}_r}$ is a suitable compression velocity, which does not affect the solution outside the interface region because of the factor $\delta (1 - \delta )$. An appropriate discretization and a suitable choice of the parameters ensure mass-conservation of the phase fraction and guarantee a bounded and unbiased solution (Berberović et al. 2009; Ferrari & Lunati 2013). This is achieved by using the multidimensional universal limiter with explicit solution solver, which is implemented in interFoam, to compute the compression velocity, ${\boldsymbol{u}_r} = \min [{C_f}|\boldsymbol{u}|,\;\max (|\boldsymbol{u}|)]\boldsymbol{n}$, which acts perpendicular to the interface due to the presence of the interface normal $\boldsymbol{n}$ (Ferrari & Lunati 2013). Here $\textrm{max}(|\boldsymbol{u} |)$ is the maximum velocity in the interface region and ${C_f}$ is a user-defined parameter that controls the strength of the compression. In our simulations, we have used ${C_f} = 1$, which guarantees a conservative compression and ensures a convergent numerical scheme: the compression term vanishes and (2.3) is recovered if the interface thickness tends to zero (Weller 2008; Berberović et al. 2009).

A rectangular buffer region with a length of 40 μm is set up at the inlet and outlet of the porous medium. At the initial state, the wetting phase fills the buffer region at the inlet, and the non-wetting phase fills other volumes. The coupled continuity and momentum equations were solved using the pressure implicit with the splitting of the operator method. At each computing step, the time step $\Delta t$ is automatically adjusted based on the Courant–Friedrichs–Lewy condition, where the maximum Courant number Co is set as 0.5 (Zhou et al. 2022). To manage the substantial computational cost associated with solving the N–S equations together with the VOF method, we focused our simulations on a reduced 3-D domain, specifically a 1/4 × 1/4 portion of the original microfluidic chip. This approach allowed us to qualitatively capture the structural complexity and the dynamics of multiphase flow processes as observed in the whole microfluidic chips. The calculation domain only contains the pore space of porous media, where an unstructured grid was used. It is worth noting that even with the current local porous structure, the calculation amount of direct numerical simulation is huge. Each case in this article approximately requires 224 cores to be calculated for more than two months, which means that it is nearly impossible to calculate the displacement in whole microfluidic porous media because the porous structure size increases 16 times, and displacement time increases four times. The discretization and meshing methods and verifications of the numerical model for capturing correct multiphase flow dynamics are presented in Appendix A. Our simulations were conducted at the High-Performance Computing Center of Tsinghua University, China.

2.4. Dynamic pore-network models

Considering the huge calculation and time consumption of direct numerical simulations for macroscopic porous media, pore-network models are intuitive and computationally inexpensive to link the pore-scale flow behaviour and macroscopic multiphase flow patterns. To simplify calculations and present the main flow dynamics, the porous geometry in the pore-network model is approximated by a pores-throats network, including series pore-throat channels and parallel pore-doublets, and the flow follows the fully developed Poiseuille flow law. The dynamic pore-network model was applied to consider the competition between viscous forces and capillary forces on displacement processes. Considering the importance of pore-scale events, such as bypass and snap-off (Singh et al. 2022), we will improve the dynamic pore-network model to link the pore-scale events and macroscopic flow patterns in porous media. Details of improved dynamic pore-network models are introduced in Appendix B.

3. Results and discussion

In this section, forced imbibition experiments were conducted in different microfluidic chips first. These microfluidic chips have the same 2-D morphology, except that the depth variation ranges from uniform depth ($\alpha = 1$), two different depths ($\alpha = 2$) to three different depths ($\alpha = 3$). During imbibition in porous media, the trapped non-wetting phase can be classified into three types: (i) isolated ganglia, which are trapped by the instability of water film flow along the corner (arc-menisci), also called corner flow; (ii) clusters, which occupy several pores bypassed by the main meniscus (the main terminal meniscus); and (iii) unswept non-wetting area, which occupies a continuous region skirted by the macroscopic preferential flow (preferential flow pathway). Through a combination of microfluidic experiments and direct numerical simulations, it was evident that depth variation, capillary number and wettability can significantly affect the non-wetting phase entrapment, which stems from microscopic flow dynamics including the instability of this corner flow and the velocity mismatch of the main meniscus in pore-doublets. Then, pore-scale theoretical models and corresponding phase diagrams were proposed and verified based on series pore-throat channels or pore doublets. Finally, to emphasize the importance and generality of these pore-scale events, these pore-scale events were inputted into an improved dynamic pore-network model, which can reproduce the above macroscopic flow patterns in microfluidic porous media.

3.1. Fluid entrapment during forced imbibition in microfluidic porous media

3.1.1. The effect of capillary number on the fluid entrapment scenario

The effect of capillary number on fluid entrapment and imbibition stability was explored in different microfluidic chips under strongly water-wet conditions ($\theta \approx 30^\circ $). The air was used as the non-wetting phase which was presaturated in the microfluidic chip before the waterflooding. Dyed water (100 ppm fluorescein sodium salt, Macklin) was used as the aqueous phase. The IFT $\gamma $ between the two phases is 73 mN m⁻¹. The dynamic viscosity of the wetting phase (water) ${\mu _w}$ is 1.0 mPa s, and the dynamic viscosity of the non-wetting phase (air) ${\mu _{nw}}$ is 0.018 mPa s. The water flow rate q was adjusted for different capillary numbers, $Ca = ({\mu _w}{u_i})/\gamma $, where ${u_i} = q/S$ is the characteristic velocity of the inlet cross-section of porous structure, which is located at the connection between the buffer region and the porous structure, $S = w\,\textrm{d}\phi $ was the area of the inlet cross-section, w = 6 mm is the width of the microfluidic porous media, d = 40 μm is the depth of buffer region and $\phi = 0.38$ is the 2-D porosity of the inlet cross-section, which is equal to the one-dimensional porosity of the inlet cross-section in the width direction. The contact angle of water on the thermal oxidation silicon wafer in the air is $\theta = 32 \pm 3^\circ $, which is measured by a contact-angle-measuring instrument (DSA25, Krüss).

Under strong imbibition ($\theta \approx 30^\circ $), we compared experimental and simulated displacement patterns in various microfluidic chips ($\alpha = 1,2,3$) at different capillary numbers (Ca = 2.6 × 10⁻⁴, 2.6 × 10⁻⁵), as shown in figure 4. At higher Ca = 2.6 × 10⁻⁴, the observed patterns exhibit a cooperative pore-filling mechanism driven by the continuous fluid–fluid main meniscus within several narrow preferential flow pathways. Upon reducing the capillary number to a lower level (Ca = 2.6 × 10⁻⁵), for the uniform-depth microfluidic chip ($\alpha = 1$), the pore bodies are predominantly occupied either entirely by the wetting phase or the non-wetting phase (figure 4a iv). In the case of $\alpha = 1$, decreasing Ca results in an increase in the invading width, a finding that aligns qualitatively with the results reported by Stokes et al. (1986). However, under lower capillary number (Ca = 2.6 × 10⁻⁵), varying from uniform-depth microfluidic chip to three-different-depth microfluidic chip destabilizes the immiscible two-phase flow (figure 4a iv–vi). For the three-different-depth microfluidic chip ($\alpha = 3$), nearly every pore body entraps isolated non-wetting ganglia, while the wetting phase occupies pore edges and small pores by sequentially coating the perimeters of neighbouring walls (figure 5a,b). What is particularly intriguing is that the fluid entrapment within depth-variable microfluidic chips contradicts the invasion dynamics reported by Stokes et al. (1986), where decreasing Ca leads to lower water saturation S_w (figure 4c,d), and numerous isolated non-wetting fluid ganglia are trapped by the water flow in the three-different-depth microfluidic chip (figure 4a vi,b vi). Consequently, an increase in depth variation $\alpha $ appears to transform the main entrapment morphology from preferential flow-induced unswept areas into unstable bypass-induced clusters or corner flow-induced isolated ganglion types under strong imbibition conditions (figure 4a iv–vi).

Figure 4. Capillary number effect on the invasion morphologies at the breakthrough stage when the wetting fluid reaches the outlet in experiments (a) and simulations (b). Two different capillary numbers were used in the experiments and simulations, with higher Ca = 2.6 × 10⁻⁴ (a i–iii, b i–iii) and lower Ca = 2.6 × 10⁻⁵ (a iv–vi, b iv–vi) under the same contact angle $\theta \approx 30^\circ $. The water saturation S_w at the breakthrough stage of microfluidic experiments (c) and simulations (d) was influenced by capillary number Ca and depth variation $\alpha $.

Figure 5. The imbibition process in the three-different-depth microfluidic experiments and corresponding direct numerical simulations under capillary number Ca = 2.6 × 10⁻⁵ and contact angle $\theta \approx 30^\circ $. (a) The global view of water displacing air in a three-different-depth microfluidic chip before the breakthrough stages. (b) Enlarged view of corresponding snapshots during forced strongly imbibition illustrating the water film flow and snap-off phenomena. (c) The simulated invasion morphologies until the invading fluid reaches the outlet of the calculation domain, where the calculated area is the same as the enlarged structure. The pore volume of microfluidic porous media is used to represent the dimensionless time of multiphase displacement.

It can be concluded that under higher Ca conditions, the main meniscus velocity far exceeds the corner flow velocity, so corner flow is submerged, and conversely, the water film flow will be triggered to trap the non-wetting phase. The water saturation S_w of microfluidic experiments and simulations at the breakthrough stage also presents that the invasion dynamics reported by Stokes et al. (1986) were primarily applicable to traditional uniform-depth microfluidic chips ($\alpha = 1$), where decreasing Ca leads to a more complete displacement. Conversely, in the case of multidepth microfluidic chips ($\alpha = 2,3$), decreasing Ca results in a pronounced fluid entrapment with a more incomplete displacement (figure 4c,d). This fluid entrapment in multidepth microfluidic chips is rarely observed by traditional 2-D microfluidic experiments (Jung et al. 2016; Singh et al. 2019; Cha et al. 2021) and 2-D pore-scale simulations (Holtzman & Segre 2015; Hu et al. 2019; Wolf, Siebert & Surmas 2020), but is frequently encountered in 3-D porous media (Chaudhary et al. 2013; Mehmani et al. 2019).

By comparing experimental invading processes and corresponding numerical simulations results in the three-different-depth microfluidic chip at the same local region (figure 4a vi,b vi), we find that the numerical method accurately captures the invading pathway and fluid entrapment behaviour, despite using only a 1/4 × 1/4 portion of the microfluidic structures due to computational limitations (figure 5). Figure 5 shows that the patterns evolve into thinner films forming patterns with branches and frequent reconnections. The thin water film forms and propagates throughout the porous media, which is far speedier than the main meniscus flow. The simulated imbibition morphologies are qualitatively consistent with the experimental results that increasing depth variation $\alpha $ promotes isolated ganglia and small cluster formation by the unstable corner flow, as shown in figure 4(b i–vi). Our numerical results effectively reflect the main trapping characteristics, determined by the competition between corner flow and the main meniscus, during fluid displacement within the microfluidic chips.

3.1.2. The effect of wettability on the fluid entrapment scenario

To investigate the impact of wettability on the entrapment of non-wetting fluids, we varied the wettability from strongly water-wet ($\theta \approx 30^\circ $) to weakly water-wet ($\theta \approx 60^\circ $). The experimental conditions remained consistent, except that dyed decane (100 ppm Nile red) was used as the non-aqueous phase to adjust the wettability to weakly water-wet conditions. The IFT between the two phases is 49 mN m⁻¹. The viscosity of the non-wetting phase (decane) ${\mu _{nw}}$ is 0.838 mPa s. The contact angle of water on the thermal oxidation silicon wafer immersed in decane is $\theta = 61 \pm 5^\circ $, which is measured by a contact-angle-measuring instrument (DSA25, Krüss). The viscous instability can also be ignored due to the viscosity ratio ($M = {\mu _w}/{\mu _{nw}} = 1.2$) being much larger than 10⁻² (Lenormand, Touboul & Zarcone 1988; Primkulov et al. 2021).

As wettability transitions from strongly water-wet to weakly water-wet conditions in depth-variable microfluidic chips ($\alpha = 2,3$), the non-wetting fluid entrapment will transform from isolated ganglia to clusters or unswept areas because the corner flow was overwhelmed (figure 6a ii,iii). Simulation results also support this trend, indicating that an increase in the contact angle $\theta $ stabilizes the imbibition interface and transfers the corner flow instability-induced snap-off to stable main meniscus flow (figure 6b ii,iii). Conversely, in uniform-depth microfluidic chips ($\alpha = 1$), decreasing the contact angle $\theta $ can promote main meniscus flow advances in small channels, resulting in the emergence of more small branches from the preferential flow pathway, as shown in both experimental results (figure 6a i,iv) and simulation results (figure 6b i,iv), which are also consistent with previous studies (Holtzman & Segre 2015; Jung et al. 2016).

Figure 6. Wettability effect on the invasion morphologies at the breakthrough stage when the wetting fluid reaches the outlet in experiments (a) and simulations (b). Two different contact angles were used in the experiments and simulations, with a higher contact angle $\theta \approx 60^\circ $ (a i–iii,b i–iii) and a lower contact angle $\theta \approx 30^\circ $ (a iv–iv,b iv–vi) under the same capillary number (Ca = 2.6 × 10⁻⁵). The water saturation S_w at the breakthrough stage of microfluidic experiments (c) and direct numerical simulations (d) was influenced by contact angle $\theta $ and depth variation $\alpha $.

The water saturation S_w of microfluidic experiments and simulations at the breakthrough stage also presents that for the uniform-depth microfluidic chip ($\alpha = 1$), the water saturation S_w increases as decreasing contact angle $\theta $, which results in a more complete displacement; in contrast, for depth-variable microfluidic chips ($\alpha = 2,3$), decreasing the contact angle $\theta $ causes a reduction in water saturation, indicating a strengthening of fluid entrapment in this scenario (figure 6c,d).

3.2. Pore-scale theoretical models and phase diagrams

Considering that local direct numerical simulations exhibit a similar trend to that observed in microfluidic chips (figures 4–6), the pore-scale mechanism of the trapped ganglia and clusters formation during the forced imbibition process in uniform-depth and depth-variable microfluidic chips plays an important role in fluid entrapment. The complex porous structure can be dissected into a series of idealized series pore-throat channels and parallel pore doublets (figure 7). Corresponding pore-scale trapping events can be classified as snap-off and bypass modes. It is worth noting that the morphology of preferential flow pathways in the microfluidic chips appears to be nearly consistent (figures 4 and 6), potentially controlled by fluid properties, such as viscous instability (Zhao et al. 2016) and elastic flow instability (Browne et al. 2023) and the spatial correlation of pore spaces (Siena et al. 2019). Here we ignored this complex formation of preferential flow pathways and mainly focused on the formation of isolated ganglia and clusters, which is dominated by capillary forces at the microscopic scale.

Figure 7. Schematic diagram of the geometry and liquid–liquid interface of (a i–iii) the series pore-throat channel and (b i–iii) pore doublets. (a i,b i) Top view of channels. (a ii,b ii) The main view of channels. (a iii,b iii) The side view of channels. (a iv,b iv) Multiphase distribution and fluid–fluid interface. The blue line is the fluid–fluid interface.

To gain physical insight into isolated ganglia formation, the snap-off criterion was proposed under strongly water-wet conditions. We develop a critical capillary number $C{a_c}$ that incorporates the contribution of the contact angle $\theta $ and quantifies the relative importance of main meniscus flow and corner flow in the series pore-throat channel. Furthermore, the depth-variation effect on the snap-off phenomenon was discussed. For the microscopic mechanism of cluster formation, the bypass criterion was proposed to separate stable and unstable (bypass) displacement patterns. We develop a critical contact angle ${\theta _c}$, which incorporates the contribution of the geometry factor $\varPsi $ and quantifies the competition relation of different main meniscus flow in parallel large and small channels.

3.3.1. Snap-off criteria in series pore-throat channels

The criterion for the occurrence of snap-off in a series pore-throat channel is determined by the instability of the water flow, which is controlled by three conditions: (i) wettability condition, $\theta < ({\rm \pi} - \varphi )/2$, where $\varphi $ is the corner angle and the contact angle in our microfluidic chips for corner flow should satisfy $\theta < 45^\circ $, which can maintain the existence of corner flow in the angular channel (Concus & Finn 1969); (ii) geometry condition (Roof 1970; Lenormand, Zarcone & Sarr 1983), $(1/{W_p} + 1/{D_p}) - (1 - \textrm{tan}\,\theta )/\textrm{min}({w_t},{d_t}) < 0$, ${W_p}$ and ${D_p}$ are the width and depth of the pore, respectively, ${w_t}$ and ${d_t}$ are the width and depth of the throat, respectively (figure 7a). It can cause the unstable interface of the water film at the throat by $\Delta {P_c} = {p_{c - p}} - {p_{c - {t_{cir}}}} < 0$, where ${p_{c - p}}$ is the capillary pressure at the main meniscus in the pore and ${p_{c - {t_{cir}}}}$ is critical capillary pressure when the corner flow interface merges at the throat; (iii) capillary number condition, $Ca < C{a_c}$, where $C{a_c}$ is the critical capillary number, which can determine the domination role of the precursor corner flow. The detailed derivation of the snap-off criterion is presented in Appendix C.

Snap-off events only occur under strongly water-wet conditions ($\theta < 45^\circ $ for rectangular channels), but the occurrence of snap-off was also controlled by the pore geometry and capillary number, as shown in figure 4. The snap-off criterion depends on the relative value of these geometry-dependent capillary pressures, when the following criterion is satisfied, snap-off will occur in the series pore-throat channel:

(3.1)\begin{equation}\Delta {P_c} = 2\gamma \,\textrm{cos}\,\theta \left[ {\left( {\frac{1}{{{W_p}}} + \frac{1}{{{D_p}}}} \right) - \frac{{1 - \textrm{tan}\,\theta }}{{\min ({w_t},{d_t})}}} \right] < 0.\end{equation}

In depth-variable microfluidic chips, a larger pore width will have a larger pore depth, thus the depth of the pore channel will be several times that of the throat, ${D_p} = \alpha {d_t}$. Therefore, compared with uniform-depth microfluidic chips ($\alpha = 1$), depth-variable microfluidic chips can satisfy the geometry conditions (3.1) more easily, described as follows:

(3.2)\begin{equation}\Delta {P_c}(\alpha \ge 1) - \Delta {P_c}(\alpha = 1) = \left( {\frac{1}{\alpha } - 1} \right)\frac{1}{{{d_t}}} \le 0.\end{equation}

It can be concluded that increasing depth variation $\alpha $ will promote the occurrence of snap-off, as shown in figure 4. However, the above geometric criteria are only valid under lower capillary numbers where the corner flow velocity is much higher than the main meniscus flow velocity. By comparing the mean flow velocity of corner flow and main meniscus flow, the critical capillary number $C{a_c}$ can be determined by

(3.3)\begin{equation}C{a_c} = 2a{(6K\beta \tau )^{ - 1/2}},\end{equation}

where $\beta = \mathcal{C}/[2{\hat{g}_w}{(\textrm{cos}\,\theta - \textrm{sin}\,\theta )^4}]$ is the dimensionless flow resistance as a function of contact angle $\theta $, where ${\hat{g}_w}$ is the dimensionless conductance of the wetting phase proposed by Patzek & Kristensen (2001). Other parameters are constant, such as a is a parameter introduced by the trial function, and the parameter K is calculated by a variational method. For a tube with rectangular corners, a and K are 0.59 and 1.447 (Ransohoff & Radke 1988), respectively. Here $\tau = {10^4}$ is the dimensionless time for snap-off phenomena based on previous experimental observation (Gauglitz, St. Laurent & Radkle 1987).

When $Ca > C{a_c}$, the corner flow from the moving meniscus will not be fully developed, thus snap-off phenomena will be suppressed. The theoretical results illustrate how $C{a_c}$ varies as a function of contact angle $\theta $ and the snap-off criterion under different contact angles $\theta $ are presented in figure 8. Figure 8(a) illustrates that the analytically predicted $C{a_c}$ agrees well with the numerical results in distinguishing between the stable and snap-off for all combinations of the capillary number Ca and contact angle $\theta $. In the simulations, a series pore-throat channel with a square cross-section was used. The geometry condition $1/{W_p} + 1/{D_p} - (1 - \textrm{tan}\,\theta )/\textrm{min}({w_t},{d_t}) < 0$ is satisfied in this pore-throat channel. A typical evaluation of two-phase distributions in the pore-throat geometry at different Ca under strongly water-wet conditions is present in figure 8(b). The competition between the main meniscus and precursor corner flow is also presented well in figure 8(b), which numerically supports our theoretical analysis in (3.3).

Figure 8. Increasing the capillary number Ca over the critical value Ca_c will suppress the snap-off event. (a) Semilog simulation plot showing the relationship between Ca and θ. The solid line indicates the analytically predicted Ca_c (3.3). The symbols are the results obtained by the direct numerical simulations, where filled and open symbols represent the snap-off and stable displacement, respectively. (b) Time evaluation of two-phase distributions in the pore-throat geometries by the direct numerical simulation at the different Ca but $\theta = 30^\circ $. The widths of the pore channel and the throat channel are 120 and 60 μm, respectively, and the depths of the pore channel and the throat channel are 80 and 40 μm, respectively.

An important finding from our simulations is that under strongly water-wet conditions ($\theta < 45^\circ $), the impact of the contact angle $\theta $ on the critical capillary number $C{a_c}$ remains relatively constant, showing a weakly linear relationship, with $C{a_c}$ typically falling within the range of around 10⁻³–10⁻⁴. The trend is qualitatively consistent but the value is slightly higher than our microfluidic experiments (figure 4) and the previous porous media studies summarized by Hu et al. (2018). This slight discrepancy may be attributed to the geometric disparities between the simplified pore-throat channel used in our simulations and the complex porous structures.

3.3.2. Bypass criteria in pore-doublet channels

Bypass events are generated by the velocity mismatch of the main meniscus flow in pore doublets, as shown in figure 7(b). Our analysis focuses solely on scenarios in which fluid–fluid interfaces within these pore doublets are connected. Here the smaller channel as an interior edge will influence the interfacial stability, which plays an important role in capillarity-driven flow. This pore-scale model aligns with our microfluidic experiments, where shallow channels surround the deeper channel in a multidepth microfluidic chip designed and the merging of fluid–fluid interfaces occurs frequently, otherwise, the previous pore network model for separate main meniscus should be considered as well as the competition between viscous force and capillary force (Chatzis & Dullien 1983; Gu et al. 2021). Moreover, we assume that the capillary force in the large channel is not affected by the small channel, but the critical value of the capillary force in the small channel is affected by the large channel, this assumption is also validated by the simulation results (figure 9b). The detailed derivation of the bypass criterion is presented in Appendix D.

Figure 9. Increasing the contact angle $\theta $ over the critical value ${\theta _c}$ will inhibit the bypass event. (a) Simulation plot showing the relationship between $\theta $ and $\varPsi $. The solid line indicates the analytically predicted critical $\textrm{cos}\,{\theta _c}$, which is a solution of the bypass criterion predicted by (3.6). The filled and open symbols represent the bypass and stable displacement, respectively. The widths of the large and small channels are set as 120, 80 μm (square); 200, 60 μm (circle); 160, 100 μm (triangle); 160, 140 μm (rhombus), respectively. The depth of the small channel is set as 40 μm, and the large channel is set as 80 μm (red) or 120 μm (blue). Blue samples will have larger $\varPsi $ than red samples. (b) Time evaluation of two-phase distributions in the paralleled large and small channels by the direct numerical simulation under different $\theta $ but $\varPsi = 0.67$. The widths of the large and small channels are set as 120 μm and 80 μm. The depths of the large and small channels are set as 80 and 40 μm.

If the capillary force difference in the large channel and small channel $\Delta {P_c} < 0$, the interface in the small channel will invade preferentially and the bypass phenomenon will occur,

(3.4)\begin{gather}\Delta {P_c} = 2\gamma \left[ {\textrm{cos}\,\theta \left( {\frac{1}{{{W_l}}} + \frac{1}{{{D_l}}}} \right) - \left( {\textrm{cos}\,\theta \frac{1}{{{d_s}}} - \frac{{1 - \textrm{cos}\,\theta }}{{2{w_s}}}} \right)} \right] < 0,\end{gather}

(3.5)\begin{gather}\textrm{cos}\,\theta > \textrm{cos}\,{\theta _c} = \frac{1}{{2{w_s}\left( {\dfrac{1}{{{d_s}}} - \dfrac{1}{{{W_l}}} - \dfrac{1}{{{D_l}}}} \right) + 1}},\end{gather}

where ${\theta _c}$ is the critical contact angle for bypass phenomena, using the dimensionless geometry factor of the parallel channels $\varPsi = 2{w_s}(1/{d_s} - 1/{W_l} - 1/{D_l})$ gives rise to the bypass criterion,

(3.6)\begin{equation}\textrm{cos}\,{\theta _c} = \frac{1}{{\varPsi + 1}}.\end{equation}

In our microfluidic experiments, the large channel has a larger width ${W_l}$, thus the depth of the large channel will be several times that of the throat, ${D_l} = \alpha {d_s}$. We find the larger depth variation $\alpha $ will lead to a higher $\varPsi $, which will result in a higher critical contact angle ${\theta _c}$,

(3.7)\begin{equation}\varPsi (\alpha \ge 1) - \varPsi (\alpha = 1) = {w_s}\frac{1}{{{d_s}}}\left( {1 - \frac{1}{\alpha }} \right) \ge 0.\end{equation}

It can be concluded that increasing depth variation $\alpha $ will promote the occurrence of the bypass event. The theoretical results illustrate how ${\theta _c}$ varies as a function of $\varPsi $ and the bypass criterion under different $\varPsi $ is presented in figure 9. Figure 9(a) illustrates that the analytically predicted ${\theta _c}$ agrees well with the numerical results in distinguishing between the stable and unstable front (bypass) for all combinations of contact angle $\theta $ and geometry factor $\varPsi $. In the simulations, pore doublets (parallel large and small channels) with a square cross-section were accepted. A typical time evaluation of two-phase distributions in the pore doublets geometry at different $\theta $ is present in figure 9(b). The competition of the main meniscus in parallel small and large channels numerically supports our theoretical analysis in (3.6) and (3.7), as shown in figure 9.

We can find that smaller $\theta $ will result in the main meniscus flow in the small channel far exceeding this in the large channel, thus the bypass mode occurs more easily in smaller $\theta $. This pore-scale phenomenon explains that decreasing contact angle leads to more small channel branches and clusters in porous media. The critical contact angle ${\theta _c}$ will decrease with increasing geometry factor $\varPsi $, which is also consistent with the flow patterns in porous media that increasing depth variation $\alpha $ can trigger the bypass mode, as shown in figure 6(a i–iii).

3.3. Linking pore-scale events to macroscopic imbibition patterns

To demonstrate the dominance of these pore-scale events (snap-off and bypass) and give a quantitative comparison between experiments and simulations, we simulated the macroscopic imbibition patterns in microfluidic porous media by the improved dynamic pore-network model, where the snap-off and bypass behaviours are incorporated in the model. The simulated macroscopic imbibition patterns and quantification analysis of wetting phase saturation are presented in figure 10. Since the pore network model can only present the saturation of each pore but not the specific fluid distribution state in the pore, we use the distribution of water saturation in each pore to represent the multiphase distribution in microfluidic porous media, as shown in figure 10(a–i). Purple represents completely occupied by water, yellow represents occupied by non-aqueous phase and the colour between purple and yellow means aqueous and non-aqueous phases occupy the same pore at the same time. The morphology of imbibition patterns simulated by the improved dynamic pore-network model is similar to that of microfluidic experiments, where increasing depth variation $\alpha $ can promote fluid entrapment, whereas increasing capillary number Ca and contact angle $\theta $ yield the opposite effect (figure 10a–i). The water saturation variation S_w affected by capillary number Ca or contact angle $\theta $ also fits with the experimental results and the difference in the water saturation value S_w between the dynamic pore-network model and microfluidic experiments is less than 5 % (figure 10j,k).

Figure 10. Linking pore-scale flow behaviour to macroscopic imbibition patterns by the improved dynamic pore-network model. (a–f) The invasion morphologies and local water saturation distribution at the breakthrough stage: higher Ca = 2.6 × 10⁻⁴ (a–c) and lower Ca = 2.6 × 10⁻⁵ (d–f) under the same contact angle θ = 30°. (d–i) The invasion morphologies and local water saturation distribution at the breakthrough stage under two different wettability: lower θ = 30° (d–f) and higher θ = 60° (g–i) under the same capillary number Ca = 2.6 × 10⁻⁵. (j,k) The quantitative comparison between microfluidic experiments and dynamic pore-network model simulations about water saturation S_w in microfluidic porous media affected by (j) capillary number and (k) contact angle at the breakthrough stage.

3.4. Implications and discussion

According to microfluidic experiments, theoretical analysis and simulations, we can conclude: (i) increasing depth variation $\alpha $ results in incomplete displacement by triggering snap-off or bypass or their combinations; (ii) higher capillary number Ca can lead to complete displacement by suppressing the snap-off event when $\theta \le 45^\circ $, and increased contact angle $\theta $ can lead to complete displacement by inhibiting the bypass event when $\theta < 90^\circ $. These microscopic pore-scale conclusions can explain the complex capillary phenomena of microfluidic experiments, which has been demonstrated by the improved pore-network model, which establishes a bridge between our understanding of interfacial behaviour at the single pore scale and the formation of displacement patterns in porous media. However, it should be noted that the capillary forces are dominant, and the viscous instability of the fluid–fluid interface is ignored in this work, so our microscopic flow dynamic mechanisms and macroscopic fluid entrapment should be revisited when the viscosity ratio $M = {\mu _w}/{\mu _{nw}}$ is smaller than 10⁻² (Lenormand et al. 1988; Primkulov et al. 2021). Thin films from the moving meniscus will lead to different entrapment scenarios considering the Saffman–Taylor instability (Levaché & Bartolo 2014; Odier et al. 2017). The coupling of viscous instability and strong imbibition will present more complex entrapment patterns even in uniform-depth microfluidic chips (Zhao et al. 2016). In addition, the surface roughness will influence the precursor film flow behaviour to trigger snap-off phenomena (Sun, Mehmani & Torres-Verdín 2021). Preferential flow conditions can also promote imbibition instability and trigger fluid entrapment (Lei et al. 2022a). The Haines jump will lead to a clear initial effect, especially under strong drainage conditions, which may result in different snap-off events (McKinley & Renardy 2011; He et al. 2017). All those problems will be analysed to modify those models in future work.

4. Conclusions

In this study, we have introduced an innovative fabrication technique for multidepth microfluidic chips with complex porous geometry. Compared with traditional uniform-depth microfluidic chips, these multidepth microfluidic chips can release the geometric confinement to mimic the wide pore size distribution as naturally 3-D porous structures. The geometric confinement-induced capillarity limitations differ in uniform-depth, two-different-depth and three-different-depth microfluidic chips. These depth-variable microfluidic chips could be an advanced tool for adjusting geometric confinements of porous media and reflecting different fluid dynamics and fluid entrapment during forced imbibition in different geometric-confined porous media.

The experiments of forced imbibition have elucidated the complex interplay of depth variation $\alpha $, capillary number Ca and wettability $\theta $ on flow behaviours and fluid entrapment. Different from the imbibition behaviour in traditional uniform-depth microfluidic chips ($\alpha = 1$), increasing depth variation $\alpha $ to release the geometric confinement can promote fluid entrapment during forced imbibition experiments, whereas increasing capillary number Ca and contact angle $\theta $ yield the opposite effect. These microfluidic experiments and corresponding direct numerical simulations enabled us to purify the complex capillarity and classify the basic microscopic transport modes as snap-off and bypass phenomena. Increasing depth variation $\alpha $ in microfluidic chips can enhance the instability of corner flow or the main meniscus. However, the higher capillary number Ca will suppress the corner flow instability-induced snap-off phenomena, and the larger contact angle $\theta $ will weaken the main meniscus flow instability-induced bypass phenomena. We proposed a theoretical model to elucidate the combined effect of geometric confinement, capillary number and wettability on the snap-off and bypass phenomena. These phase diagrams predicted by the theoretical model describe the pore-scale imbibition behaviour transition from stable flow to snap-off or bypass modes, which were numerically verified. By incorporating these pore-scale events into pore-scale flow behaviours, we establish the link between the microscopic flow dynamics and macroscopic flow patterns in porous media.

Our study not only sheds light on the coupling effects of 3-D pore structure, capillary number and wettability on non-wetting phase entrapment in porous media but also delves into the dynamics of the microscopic main meniscus and water film, offering insights into macroscopic flow patterns. This analysis provides the understanding of CO₂ storage in the deep stratum for hydrocarbon reservoirs, contaminant transport in underground aquifers and hydrocarbon recovery processes, the optimization of multiphase behaviour during liquid-infused material fabrication or operational guidelines for the 3-D geometry design of microfluidic chips for controlling multiphase flow patterns.