# Impact of wetting films on the stability of two-phase flow in porous media: A pore-doublet perspective

Nathan Bernard, Cyprien Soulaine, Sophie Roman

► To cite this version:

Nathan Bernard, Cyprien Soulaine, Sophie Roman. Impact of wetting films on the stability of two-phase flow in porous media: A pore-doublet perspective. Advances in Water Resources, 2025, 204, pp.105060. 10.1016/j.advwatres.2025.105060. hal-05202641

HAL Id: hal-05202641

https://hal.science/hal-05202641v1

Submitted on 7 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Impact of wetting films on the stability of two-phase flow in porous media: A pore-doublet perspective

Nathan Bernard<sup>1,0</sup> *, Cyprien Soulaine<sup>1,0</sup>, Sophie Roman<sup>1,0</sup>

Institut des Sciences de la Terre d'Orléans, ISTO, UMR 7327, Université d'Orléans, CNRS, BRGM, OSUC, 45071, France

## ARTICLE INFO

### Keywords:

- Porous media

- Two-phase flow

- Pore-doublet model

- Microfluidics

## ABSTRACT

Two-phase flow stability in porous media depends on a complex interplay between surface tension effects and viscous forces. We investigate the role of wetting films on the stability of two-phase flow during drainage. Such films are ubiquitous in porous media and are deposited at the solid walls and corners after the passage of the meniscus. We derive a novel dynamic pore-doublet model with identical channels that includes the effects of lubricating films. Direct comparisons with microfluidic experiments highlight that neglecting wetting films in pore-doublets analysis leads to significant prediction errors. We demonstrate that pore-doublet models with wetting films reproduce stability phase diagrams similar to the so-called Lenormand diagrams. At high viscosity ratios, our model predicts phase transition between stable, viscous, and capillary fingering regimes that agree very well with Lenormand phase diagrams obtained using complex porous structures. Our study demonstrates that using simple systems, such as pore doublets with lubricating films, allows us to capture the non-linear response of two-phase flow stability similar to those observed in complex pore networks.

## 1. Introduction

Multiphase flow in porous media is of utmost importance in many industrial applications including carbon capture and storage (Kumar et al., 2020), distillation columns (Soulaine et al., 2014), soil remediation (Delshad et al., 1996), enhanced oil recovery (Moore and Slobod, 1955) or liquid composite molding (Lundström et al., 2008). Two-phase flows rely on complex coupled interactions between the two immiscible fluids and the porous matrix. The wettability of the solid plays an important role in the fluid distribution: the wetting fluid tends to flow in contact with the solid, along the walls of the pore or in the roughness of the porous matrix while the non-wetting fluid flows in the pore center minimizing the surface of contact with the solid (Blunt, 2017; Gredicak et al., 2025). Understanding the fluids' behavior is essential for predicting the residual saturation after injection.

In porous media, two-phase flow is inherently complex, primarily due to capillary effects (Norouzisadeh et al., 2024) and viscous coupling between the phases (Roman et al., 2020; Li et al., 2005). This complex system exhibits a non-linear behavior from which instabilities arise and significantly alter the fluid flow. A typical example is the Saffman-Taylor instability (Saffman and Taylor, 1958; Tabeling et al., 1987), which appears when a more viscous fluid is displaced by a less viscous one, leading to the formation of the well-known viscous fingering patterns. In porous media and neglecting gravity effects, Lenormand et al. (1988) identified three regimes of drainage depending on two

dimensionless numbers, namely the capillary number, $Ca$ (the ratio of viscous forces and capillary forces), and the viscosity ratio, $M$ (viscosity of the injected phase over the viscosity of the defending phase). For $M > 10$ and $Ca > 10^{-1}$, the viscous forces in the injected fluid dominate the other forces and the displacement is stable. In this regime, the invasion pattern presents a uniform front and the displacement efficiency is high, meaning that only very little wetting fluid is trapped behind the interface. For $M < 10^{-3}$ and $Ca > 10^{-8}$, the viscous forces in the displaced fluid dominate, leading to a finger-like pattern, due to Saffman-Taylor instability (Saffman and Taylor, 1958). This regime has poor displacement efficiency, and an important volume of wetting fluid is left behind the invasion front. For $M > 10^{-4}$ and $Ca < 10^{-6}$, the capillary forces dominate, leading to capillary instability (Lenormand and Zarcone, 1985) in which the pores with the lower capillary entry pressure are invaded preferentially, resulting in a seemingly random invasion pattern and a high residual saturation.

These three regimes – stable front, viscous fingering, and capillary fingering – have been studied extensively at the pore scale (Dullien, 2012). The two-phase flow patterns have been observed experimentally using microfluidic devices that replicate the structure of pores (Lenormand and Zarcone, 1985; Lenormand et al., 1988; Lenormand, 1990; Zhang et al., 2011; Zhao et al., 2016; Chen et al., 2023) or using core experiments combined with microtomography imaging (Chen et al., 2018). Pore-scale simulations have been proposed to investigate the

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: nathan.bernard@cnrs-orleans.fr (N. Bernard).</footnote>
Nomenclature

$\alpha$ Growth rate of the channel width or radius.

$\beta$ Numerical constant for the calculation of the rectangular channel permeability

$\chi_i$ Dry area in channel $i$

$\Delta P_i$ Pressure difference between the beginning and the end of channel $i$

$\Delta P_{cap}$ Capillary pressure between fluid I and II.

$\Delta P_{visc,j}$ Viscous drop in fluid $j$

$\Delta z$ Position difference between the two menisci

$\epsilon$ Aspect ratio of the rectangular channels

$\kappa_i$ Curvature of the interface in channel $i$

$\lambda_i$ Slip length in channel $i$

$\mu_I$ Viscosity of the wetting fluid

$\mu_{II}$ Viscosity of the non-wetting fluid

$\nabla P_{i,j}$ Pressure gradient in fluid $j$ for channel $i$

$\phi$ Parameter related to the viscous drop in the channels

$\psi_i$ Wet fraction in the asymptotic case $Ca \to 0$

$\sigma$ Interfacial tension

$\tau$ Initial difference in the menisci position imposed for the simulations

$\theta_i$ Contact angle measure through the wetting phase in channel $i$

$\zeta$ Dimensionless geometrical factor

$A_i$ Total cross section of channel $i$

$A_{i,j}$ Cross section of the channel $i$ occupied by fluid $j$

$c$ Numerical constant used for the conductance in rectangular channel

$Ca$ Capillary number of the whole system

$Ca_i$ Local capillary number in channel $i$

$D_0$ Constant depth of the rectangular channel

$e_i$ Film thickness in the channel $i$ for the cylindrical channel

$e_{avr,i}$ Average thickness of the film in the rectangular case, calculated from the wetted area in the channel $i$

$e_{r,j}$ Thickness of the films at the center of the walls for the rectangular channel $i$

$g_{i,j}$ Conductance of the channel $i$ to the fluid $j$

$k_{i,j}$ Permeability of the channel $i$ to the fluid $j$

$L$ Total length of the channel

$M$ Viscosity ratio

$Q_i$ Flow rate of non-wetting fluid in channel $i$.

$Q_p$ Total flow rate of non-wetting fluid injected in the inlet

$r$ Radius of the cylindrical channel

$r_0$ Initial radius of the cylindrical channel

$R_{h,i}$ Hydraulic radius of the rectangular channel

$t'$ Dimensionless time

$t'_{end}$ Dimensionless simulated time

$U_i$ Velocity of the meniscus in channel $i$

$u_{i,j}$ Velocity of fluid $j$ inside of channel $i$ at a given time and position

$w$ Width of the channel

$w_0$ Initial width of the channel

$w_f$ Width at the end of the channel in the rectangular case

$z_i$ Distance from the inlet to the position of the meniscus in channel $i$

conditions leading to the emergence of flow instabilities. Some of them rely on regular pore-network modeling (Lenormand et al., 1988; Cottin et al., 2010). With the continuous improvement of computational resources, more recent studies directly track the fluid-fluid interface in the pore-space by solving the Navier-Stokes equations (Ferrari and Lunati, 2014). More recent studies simulate a drainage in a 3D

reconstruction of the rock pore structure obtained with microtomography (Tsuji et al., 2016). All these studies recover the three different flow regimes. The regime boundaries, however, vary from study to study and depend strongly on the pore structure.

Recently Al-Housseiny et al. (2014) and Mansouri-Boroujeni et al. (2023) used pore doublets to investigate the stability of two-phase flow in porous media. Pore doublets represent the most basic constitutive element of porous media, i.e. a solid grain sandwiched in between two parallel pores (Moore and Slobod, 1955). This theoretical tool was initially designed to assess the oil recovery mechanisms in a petroleum reservoir during waterflooding (Rose and Witherspoon, 1956; Chatzis and Dullien, 1983). Dynamic pore-doublets models are made of coupled ordinary differential equations solving for the fluid-fluid interface position in the two parallel pores that can be of various sizes, shapes, and lengths (Dullien, 2012; Protière et al., 2010). Some models also account for inertia effects (Mansouri-Boroujeni et al., 2023; Moebius and Or, 2014; Sorbie et al., 1995). Modern usage includes the investigation of water infiltration in soils (Moebius and Or, 2012), the assessment of CO₂ storage capacity in aquifers (Mansouri-Boroujeni et al., 2023), and the modeling of liquid composite molding (Aaboud et al., 2019).

The knowledge of the meniscus positions informs on the two-phase flow stability in the pore doublet. In principle and by extrapolation, it should inform on the stability of larger-scale porous media as well. Al-Housseiny et al. (2014) studied drainage in a symmetrical pore doublet. They demonstrated that for two identical straight channels, the system stability depends only on the viscosity ratio in agreement with Saffman-Taylor instability. They managed to recover a transition between stable and unstable flow that depends on the viscosity ratio and the capillary number by using channels whose radius increases linearly with the distance to the inlet. The dependency in $Ca$ is due to the cross-section variations of the channels, as it changes the interface curvature when the meniscus moves, thus making the interfacial tension play a role in the pore-doublet invasion, according to Young-Laplace's law. Nevertheless, the phase diagram that they derived shows some discrepancies with the regime observed by Lenormand and others.

We posit that the transition between stable front, viscous fingering and capillary fingering can be properly captured by pore doublets if the wetting layers are considered. Wetting films are a crucial aspect of two-phase flow in porous media, yet they are often not included in pore-doublet models for the sake of simplification (Blunt, 2017; Dullien, 2012; Soulaine et al., 2021). During drainage, they deposit behind the front meniscus along the solid walls (Taylor, 1961). Their thickness depends on the balance of viscous and capillary forces (Bretherton, 1961; Aussillous and Quéré, 2000; Zhao et al., 2018; Balestra et al., 2018). Those layers play a role in the fluid dynamics via increased phase connectivity (Tuller and Or, 2001) and snap-off (Roman et al., 2017). During drainage, they also act as lubricants via viscous coupling at the fluid-fluid interface (Gong et al., 2023; Li et al., 2005; Yiotis et al., 2007; Rose, 2000; Roman et al., 2020; Picchi and Battiato, 2018). Shams et al. (2018) calculated the flow conductance in channels of various shapes including the effect of wetting layers in corners for
various viscosity ratios. The resulting conductance is reminiscent of flow with partial slip conditions at the wall (Pasquier et al., 2017), in which the slip length depends on the repartition of the fluids inside the channel (i.e. Ca) and the mutual shear stress (i.e. M). In this work, we develop a new pore-doublet model that accounts for the film effects at the wall via an effective slip length. We use the model to investigate two-phase flow stability in symmetrical pore doublets before confronting the results with microfluidic experiments.

The paper is organized as follows. First (Section 2), we introduce a new theory of pore doublet for drainage and imbibition including wetting layers at the solid surface. Then, in Section 3, we describe the microfluidic experiments used to probe emerging flow patterns during drainage and to verify the theory. Finally (Section 4), we discuss the role played by wetting films on the development of instabilities in porous media. We close with conclusions and perspectives (Section 5).

## 2. Theory of dynamic pore doublets with wetting films

In this section, we describe the mathematics of two-phase flow in pore-doublets and show how they are used to assess stability criteria. In particular, we introduce the equations controlling the interface dynamics when wetting layers are involved in both circular and rectangular sections.

### 2.1. Fundamental principles of pore doublets

Pore doublets offer a simplified geometrical representation of the main features of porous media. It consists of a channel that splits into two branches, 1 and 2, upstream of a solid obstacle. The two branches rejoin each other downstream after a certain distance. Fig. 1 illustrates a pore doublet with identical branches. In this particular example introduced by Al-Housseiny et al. (2014), the branch radius varies linearly along the channel mid-axis. The cross-section of each branch, noted $A_i$, with $i = (1,2)$ is thus a function of z, the distance along the channel. We have,

$$A_i = f(z), \qquad (1)$$

where f depends on the geometry of the channels.

During drainage in pore doublets, a non-wetting fluid, hereafter noted II is injected in a system saturated with a wetting fluid, denoted I. A meniscus is formed and moves from the injection channel to the bifurcation. According to the flow conditions, the meniscus splits into two menisci that invade the child branches simultaneously or at different rates, or invade only one of the two branches (Mansouri-Boroujeni et al., 2023). The position of the menisci, $z_i(t)$ are considered at the center of the channel, see Fig. 1 (Al-Housseiny et al., 2014; Moebius and Or, 2014). The knowledge of the meniscus position to the inlet, $z_i(t)$ with $i = 1,2$, in both branches and particularly the position difference, $\Delta z = z_1(t) - z_2(t)$, characterizes the interface front stability at larger scales. The case $\Delta z \to 0$, corresponds to a stable regime while a divergence of $\Delta z$ is characteristic of an unstable regime.

Flow in pore doublets is ruled by two fundamental principles:

* The flow rate is conserved throughout the system,

$$Q_p = Q_1(t) + Q_2(t), \qquad (2)$$

where $Q_p$ is the flow rate of the injected fluid in the parent channel, considered as constant in time, and $Q_i$ ($i = 1,2$) are the flow rates in the two child branches,

* The total pressure difference between the inlet and outlet is the same in the two branches,

$$\Delta P(t) = \Delta P_1(t) = \Delta P_2(t), \qquad (3)$$

where $\Delta P_i(t)$ is the pressure drop in branch $i = 1,2$.

These two principles allow the derivation of coupled ordinary differential equations (ODE) for the evolution of the meniscus positions in both branches.

In each branch, $i = 1,2$, the flow rate of fluid II is the same at any points upstream of the meniscus, with,

$$Q_i(t) = A_{i,II}(t, z) u_{i,II}(t, z), \qquad (4)$$

where $A_{i,II}(t, z)$ is the surface area occupied by the fluid II in a cross-section, and $u_{i,II}(t, z)$ is the speed of fluid II averaged on the section $A_{i,II}$ of channel $i$. Because of the wetting films left behind during the drainage of the channel or ahead of the meniscus during imbibition, $A_{i,j}(t, z_i)$ is not necessarily equal to the area of the channel cross-section. During a drainage, $A_{i,II}(t, z)$ is smaller than the channel cross-section due to wetting layers (see Fig. 1). $A_{i,j}(t, z_i)$ is also considered to be time dependent to account for thickening or thinning of the films, based on the evolution of the capillary numbers in each channel. One of the two differences between ours and the other dynamic pore-doublet models proposed in the literature is the differentiation between the surface area upstream and downstream of the meniscus due to the existence of wetting films (Protière et al., 2010; Moebius and Or, 2012; Al-Housseiny et al., 2014; Mansouri-Boroujeni et al., 2023).

The section-averaged flow velocity, $u_{i,j}(t, z)$ is driven by the pressure gradient, $\partial_z P_{,II}$, upstream the meniscus, and by $\partial_z P_{,I}$ downstream the meniscus. It is obtained by averaging the Stokes momentum equation over a channel cross-section considering either a no-slip condition if there is no film (i.e. if $z \in [z_i(t); L]$ for a drainage) or a partial slip condition in the presence of wetting layers (i.e. if $z \in [0; z_i(t)]$ for a drainage). The partial slip condition results from the viscous shear stress between the wetting film and the non-wetting fluid. A general formulation reads (Bacri et al., 1990),

$$u_{i,II}(z) = - \frac{g_{i,II}(z) A_{i,II}(t, z)}{\mu_{II}} \frac{\partial P_{,II}}{\partial z}, \quad \text{if } z \in [0; z_i(t)], \qquad (5)$$

$$u_{i,I}(z) = - \frac{g_{i,I}(z) A_{i,I}(t, z)}{\mu_{I}} \frac{\partial P_{,I}}{\partial z}, \quad \text{if } z \in [z_i(t); L], \qquad (6)$$

where $g_{i,j}(z) = k_{i,j}(z)/A_{i,j}(t, z)$ is the channel conductance, and $k_{i,j}(z)$ is the channel apparent permeability. $\mu_I$ and $\mu_{II}$ are respectively the dynamic viscosity of fluid I and fluid II. Unlike other dynamic pore-doublet models (Moore and Slobod, 1955; Moebius and Or, 2012; Al-Housseiny and Stone, 2013; Mansouri-Boroujeni et al., 2023), Eqs. (5)–(6) account for the impact of the wetting films at the wall if they are present. Indeed, the apparent permeability and conductance not only depend on the section shape but are a function of the film thickness and viscosity ratio as well. It constitutes the second main difference with other models. We will see in the next sections different conductances for various cross-section shapes in the presence of wetting films.

In pore-doublet models, the section-averaged velocity of the bulk of the fluid j is estimated using the speed of the meniscus, $\frac{dz_i}{dt}$. We have,

$$u_{i,j}(t, z_i) = U_i(t) = \frac{dz_i(t)}{dt}. \qquad (7)$$

The pressure drop in each branch is obtained using,

$$\Delta P_i(t) = \Delta P_{visc,II}(z_i(t)) + \Delta P_{cap}(z_i(t)) + \Delta P_{visc,I}(z_i(t)), \qquad (8)$$

where $\Delta P_{visc,II}(z_i(t))$ and $\Delta P_{visc,I}(z_i(t))$ are respectively the viscous pressure drop between the inlet and the meniscus location and between the meniscus and the outlet. $\Delta P_{cap}(z_i(t))$ is the pressure drop due to the curvature of the meniscus. The capillary pressure in channel $i$, $\Delta P_{cap}$, is due to the meniscus curvature and is estimated by Young–Laplace’s law. We neglect the dynamic aspect of the curvature in this study, as we restrict ourselves to low capillary numbers. We then have,

$$\begin{aligned} \Delta P_{visc,I}(z_i) &= \int_{z_i(t)}^{L} \frac{\mu_I u_{i,I}(z)}{g_{i,I}(z) A_{i,I}(t, z)} dz \\ &= \mu_I A_{i,I}(t, z_i) \frac{dz_i}{dt} \int_{z_i(t)}^{L} \frac{1}{g_{i,I}(z) A_{i,I}(t, z)^2} dz, \end{aligned} \qquad (9)$$
<figure><img src="image_3.png" /><figcaption>Fig. 1. Schematic of the pore doublet used in this study. The channels are either cylindrical or rectangular. For this particular example, $\alpha$ is taken to be positive (diverging channels), but $\alpha$ can also be negative, representing converging channels.</figcaption></figure>

$$ \Delta P_{visc,II}(z_i) = \int_0^{z_i(t)} \frac{\mu_{II} u_{i,II}(z)}{g_{i,II}(z) A_{i,II}(t, z)} dz \\ = \mu_{II} A_{i,II}(t, z_i) \frac{dz_i}{dt} \int_0^{z_i(t)} \frac{1}{g_{i,II}(z) A_{i,II}(t, z)^2} dz, \quad (10) $$

$$ \Delta P_{cap} = \sigma \cos(\theta_i) \kappa_i(t, z_i), \quad (11) $$

where $\kappa(t, z_i)$ is the curvature of the interface in channel i, $\sigma$ is the surface tension, and $\theta_i$ is the contact angle between fluid II and the channel's walls. The viscous pressure drops are obtained using Eqs. (5)–(6) that relate the pressure gradient to the mean velocity, and because, $U_i(t) = Q_i(t)/A_{i,j}(t, z_i)$ and $Q_j$ being constant along the branch length at a time t, it can be pulled out the integral. Finally, we have,

$$ \Delta P_i(t) = \frac{dz_i}{dt} \left( \mu_{II} A_{i,II}(t, z_i) \int_0^{z_i(t)} \frac{dz}{A_{i,II}^2(t, z) g_{i,II}(z)} + \mu_I A_{i,II}(t, z_i) \int_{z_i(t)}^{L} \frac{dz}{A_{i,II}^2(t, z) g_{i,II}(z)} \right) + \sigma \cos(\theta_i) \kappa_i(z_i). \quad (12) $$

The ODEs governing the evolution of the meniscus location in both branches are obtained by combining Eqs. (4)–(7)–(12) with the two pore-doublet fundamental principles (i.e. conservation of the flow rate at the bifurcation, and same pressure drop in both branches). Hereafter, we form non-dimensionalized ODEs using $z' = z/L$, $t' = t/t_c$, $P' = P_l/P_c$, $A'_{i,j} = A_{i,j}/A_c$, $k'_{i,j} = k_{i,j}/k_c$, $g'_{i,j} = g_{i,j}A_c/k_c$, $A'_{i,j} = A_{i,j}/A_c$, $\kappa'_{i,j} = \kappa_i/\kappa_c$. The reference variables are noted with subscript c. $A_c$, $k_c$, and $\kappa_c$ depend on the channel geometry and are given in the following sections for cylindrical and rectangular channels. We define the characteristic time as $t_c = A_c L/Q_p$, the reference pressure as $P_c = L\mu_1 Q_p/A_c k_c$, the viscosity ratio as $M = \mu_{II}/\mu_l$, and the capillary number as $Ca^* = \frac{\mu_1 L Q_p}{\alpha k_c A_c \kappa_c} = Ca \zeta^{-1}$ with $Ca = \frac{\mu_1 Q_p}{\sigma A_c}$ and $\zeta^{-1} = \frac{L}{k_c \kappa_c}$ a geometrical factor.

The dimensionless conservation of the injected flow rate reads,

$$ 1 = A'_{1,II} \frac{dz'_{1}}{dt'} + A'_{2,II} \frac{dz'_{2}}{dt'}. \quad (13) $$

Because the two channels share the same inlet and outlet (see Fig. 1), the total pressure difference in channel 1 is the same as in channel

2. That means that from Eq. (12), we have

$$ \frac{dz'_{1}}{dt'} = \frac{\phi(z'_{2}) + \frac{\zeta}{Ca}(\cos(\theta_{2})\kappa'(t', z'_{2}) - \cos(\theta_{1})\kappa'(t', z'_{1}))}{A'_{1,II}(t', z'_{2})(\phi(z'_{1}) + \phi(z'_{2}))}, \quad (14) $$

$$ \frac{dz'_{2}}{dt'} = \frac{\phi(z'_{1}) - \frac{\zeta}{Ca}(\cos(\theta_{2})\kappa'(t', z'_{2}) - \cos(\theta_{1})\kappa'(t', z'_{1}))}{A'_{2,II}(t', z'_{1})(\phi(z'_{1}) + \phi(z'_{2}))}, \quad (15) $$

with

$$ \phi(t', z'_{i}) = M \int_{0}^{z'_{i}(t')} \frac{dz}{A'_{i,II}(t', z)^2 g'_{i,II}(z)} + \frac{A'_{i,II}(t', z')}{A'_{i,II}(t', z') \int_{z'_{i}(t')}^{L} \frac{dz}{A'_{i,II}(t', z)^2 g'_{i,II}(z)}}, \quad (16) $$

Eqs. (14) and (15) form a generic system of coupled ODEs that solves for the meniscus location in both branches regardless of length and channel geometry, as long as the channels fall within the realm of pore-doublet theory. In the next two parts, we calculate the conductance $g_{i,j}$, the functions $\phi(z)$, the section area $A_{i,j}(z)$, and the meniscus curvature, $\kappa(z)$ for a drainage ($\theta < 90^\circ$) with wetting films for cylindrical and rectangular channels as illustrated in Fig. 1.

## 2.2. Converging/diverging cylindrical channel with wetting layers

In this section, we derive the $\phi_i$ functions in the case of drainage in converging or diverging cylindrical channels. The channels are identical, and their radius is described by $r(z) = r_0 + \alpha z$ (see Fig. 1). No assumption is made on the sign of $\alpha$, as long as it is small enough to fit in the lubrication theory, which is the case for $|\alpha| \ll 1$ (Al-Housseiny et al., 2014). The reference variables for this geometry are $A_c = \pi r_0^2$, $k_c = r_0^2/8$, $g_c = 1/8\pi$, and $\kappa_c = 1/r_0$. The capillary number is, therefore, defined as $Ca^* = \frac{8\mu_1 L Q_p}{\sigma \pi r_0^2} = Ca \zeta^{-1}$, with $\zeta^{-1} = \frac{8L}{\pi r_0}$.

In non-dimensionalized form, the channel radius reads,

$$ r'(z') = 1 + z' \alpha', \quad (17) $$

where $\alpha' = \alpha \frac{L}{r_0}$. The cross-section is,

$$ A'_{i}(z') = r'(z')^{2}, \quad (18) $$
and the conductance reads,

$$g'_{i} = 1. \qquad (19)$$

The latter is obtained from the permeability of the cylindrical channel $k(z) = r^2(z)/8$.

The wetting fluid layer, also known as Bretherton film (Bretherton, 1961; Taylor, 1961) has a thickness $e_i(z)$ that can be approximated using the results of Balestra et al. (2018),

$$\frac{e'_{i}(t', z')}{r'(z')} \approx \frac{P(M)(3Ca_{i}(t'))^{2/3}}{1 + P(M)Q(M)(3Ca_{i}(t'))^{2/3}}, \qquad (20)$$

where $P(M)$ and $Q(M)$ are given in Appendix A, and $Ca_i(t') = Ca \frac{U_i(t')}{U_i(t') + U_2(t')}$ is the capillary number of channel i. Eq. (20) is a generalization of the law derived by Aussillous and Quéré (2000), which assumes the non-wetting fluid to be inviscid. It agrees with the equation of Aussillous and Quéré (2000) for $M < 1$, and predicts an increase of the film's thickness with the viscosity ratio until it reaches a plateau for $M > 10^3$ where the film is $2^{2/3}$ the thickness of the film in the inviscid case. This law is valid for capillary numbers below 1. In practice, $e_i$ is small compared with the channel radius (e.g. if $Ca_i = 10^{-3}$, then $\frac{e_i}{r_0} \approx 10^{-2}$).

The wetting films reduce the space available for the flow of fluid II in the channel. We introduce $\chi_i(z_i)$, as the fraction of the cross-section occupied by the non-wetting fluid in the center of the channel. For a drainage, we have

$$\chi_i(t, z) = \begin{cases} \frac{A'_{\text{II}}}{A'_i} = \left(1 - \frac{e'_i(t')}{r'(z)}\right)^2 & \text{for } z \in [0; z'_i], \\ 0 & \text{for } z \in [z'_i; 1]. \end{cases} \qquad (21)$$

The wetted area is then defined as $1 - \chi_i(t, z)$. The surface area and conductance downstream of the meniscus reads,

$$A'_{i,II} = A'_i \text{ and } g'_{i,I} = g'_i, \qquad (22)$$

while in the presence of wetting films upstream the meniscus we have,

$$A'_{i,II} = \chi_i A'_i \text{ and } g'_{i,II} \approx g'_{i,1} \left( 1 + 4 \frac{\lambda_i(t')}{r'(z')} \right), \qquad (23)$$

with $\lambda_i = e'_i(t')M$. The equation for $g'_{i,II}$ is obtained by averaging the velocity profile obtained by solving Stokes equations in a cylindrical channel with a uniform wetting film that is only driven by mutual shear stress (Bacri et al., 1990) and assuming that the film thickness is small compared to the cylinder radius. Mathematically, Eq. (23) can also be obtained by considering a partial slip condition at the fluid-fluid interface with the effective slip length $\lambda_i$ describing the lubrication caused by the wetting film.

Eq. (23) shows that the surface area of contact between the two fluids as well as the viscosity ratio play a role in the dynamics of the non-wetting fluid, highlighting the viscous coupling between the two fluids (Ehrlich, 1993; Li et al., 2005). When $M$ is very small, the wetting layers behave like a wall, and the partial slip becomes a no-slip condition at the interface between the fluid. On the opposite, when $M$ is very high, the wetting layer behaves like a lubricant, decreasing the viscous drop. This notion of slip is used later on to build the conductance in a rectangular wetted channel by analogy with the cylinder.

The curvature of the meniscus reads,

$$\kappa'(z'_i) = \frac{2}{r'(z'_i) - e'_i(t')}. \qquad (24)$$

Finally after integration (see Appendix B), we obtain, to the first order in $e_i$,

$$\phi(t', z') = \frac{r'(z')^2}{(r'(z') - e'_i(t'))^2} \times \frac{1 - \frac{r'(z')^3}{r'(1)^3}}{3\alpha' r'(z')^3} + \frac{M}{\alpha' K^3} \left( -\ln r'(z') - K \left( 1 - \frac{1}{r'(z')^2} \right) - \frac{K^2}{2} \left( 1 - \frac{1}{r'(z')^2} \right) + \ln \left| \frac{r'(z') + 2\lambda_i}{1 + 2\lambda_i} \right| \right), \qquad (25)$$

with $K = -4(M - 1)e'_i(t')$.

## 2.3. Converging/diverging rectangular channel with wetting layers

In the last section, we saw how the wetting layers act as lubricants in cylindrical channels. Porous media and microfluidics devices, however, rarely have cylindrical cross-sections and often present corners and crevices in which fluid accumulates leading to thicker wetting layers in the corners. In this section, we consider the drainage in rectangular channels whose width vary linearly from the inlet to the outlet. The two channels are identical: their width is described by $w(z) = w_0 + \alpha z$ (see Fig. 1), and their depth is constant $D(z) = D_0$. The aspect ratio between the channel length and width is $\epsilon(z) = w(z)/D(z)$, and the reference variables are $w_c = w_0$, $D_c = D_0$, $A_c = w_0 D_0$, $k_c = w_0^2/12$, $g_c = \frac{w_0}{6D_0}$, and $\kappa_c = 2/w_0$. The capillary number is defined as $Ca^* = \frac{6\mu L Q \rho}{\sigma D_0 w_0^2} = Ca \zeta^{-1}$ with $\zeta^{-1} = 6L/w_0$. In our pore doublet, the cross-section at the inlet of the channel is square and $D_0 = w_0$. Because our channel depth is constant, $D' = 1$. In non-dimensionalized form, the channel width reads,

$$w'(z') = 1 + z'\alpha', \qquad (26)$$

where $\alpha' = \alpha \frac{L}{w_0}$. The cross-section is,

$$A'_{i}(z) = w'(z')D', \qquad (27)$$

and the conductance reads,

$$g'_{i}(z) = \frac{1}{w'(z)} \left( 1 - \frac{\beta}{w'(z)} \right). \qquad (28)$$

The latter is obtained because the permeability of a rectangular channel is $k(z) \approx \frac{D_0^2}{12} \left(1 - \beta \frac{D_0}{w(z)}\right)$ with $\beta = 6(2/\pi)^5$. This equation comes from the first term of the Fourier series giving the solution for the velocity profile. This approximation is valid for $D_0 < w(z)$, and leads to a deviation from the exact solution of less than 10% for $D_0/w(z) < 0.7$ (Al-Housseiny et al., 2014).

In rectangular cross-sections, the wetting film is no longer axisymmetric. It is thin along the walls and gets thicker near the corners, as illustrated in Fig. 1. At the center of the walls, the thickness is given by Magnini et al. (2022),

$$\frac{e'_{r,i}(t', z')}{R'_{h}(z')} = \frac{0.7Ca_{i}^{2/3}(t')}{1 + 3.35Ca_{i}^{2/3}(t')}, \qquad (29)$$

where $R'_{h}(z') = \frac{w'(z')D'}{w'(z') + D'}$ is the hydraulic radius. This estimation of the film thickness is no longer true at the corners where a much larger volume of fluid accumulates (Blunt, 2017; Dullien, 2012). We assess the reduction in the cross-section occupied by the non-wetting fluid by considering an effective uniform film whose mean thickness is given by

$$\frac{e'_{\text{avr},i}(t', z')}{R'_{h}(z')} = \frac{(1 - \chi_i(t', z'))}{2}, \qquad (30)$$

where $1 - \chi_i(t', z')$ is the cross-sectional fraction occupied by the wetting fluid. Magnini et al. (2022) used numerical simulations to propose the
semi-analytical law,

$$ \chi_i(t', z) = \begin{cases} \frac{A'_{i,II}}{A'_i} = 1 - \psi(z) (1 - (1 - \chi_\infty)) \\ \exp \left( a_1 \left( Ca_i(t') \left( 1 + \frac{\epsilon(z')^2}{\epsilon_i^2} \right)^{a_2} \right) \right) & \text{for } z \in [0; z'_i], \\ 0 & \text{for } z \in [z'_i; 1], \end{cases} \quad (31) $$

where $\psi(z)$ is the wet fraction in the asymptotic case $Ca \to 0$, given in Appendix A, from Wong et al. (1995), with $\epsilon_i = 5.5$, $\chi_\infty = 0.37$, $a_1 = -0.24$ and $a_2 = -0.5$. Note that the influence of the viscosity ratio is neglected in the calculation of $e_{r,i}$, $e_{avr}$, and $\chi_i(t', z')$.

The surface area and the conductance downstream of the meniscus read,

$$ A'_{i,I}(z') = A'_{i}(z') \text{ and } g'_{i,I}(z') = g'_{i}(z'). \quad (32) $$

By analogy with the cylindrical case, the surface area and the conductance in the presence of wetting films upstream the meniscus we have,

$$ A'_{i,II}(t', z') = \chi_i(t', z') A_i(z') \text{ and } g'_{i,II}(z') = g_i(z') \left( 1 + c \frac{\lambda_i(t', z')}{R'_h(z'_i)} \right). \quad (33) $$

In Eq. (33), c is a numerical constant that weakly changes with the aspect ratio. For a rectangular shape (Ebert and Sparrow, 1965) computed $c \approx 6$. At first approach, we take the effective slip length as equal to the effective mean film thickness,

$$ \lambda_i(t', z') = e'_{avr,i}(t', z') M. \quad (34) $$

This is a significant assumption, but this allows us to find a value for $\phi$ in the rectangular case while keeping the calculations relatively simple. The mean curvature of the meniscus is approximated using (Magnini et al., 2022)

$$ \kappa'_{i}(t', z') = 2 \left( \frac{1}{w'(z') - 2e'_{r,i}(t', z')} + \frac{1}{D' - 2e'_{r,i}(t', z')} \right), \quad (35) $$

where $e_{r,i}(t', z)$ is the thickness of the films at the center of the walls as shown in Fig. 1.

Finally, after a first-order Taylor expansion in $e$ and some manipulation, we obtain (see the details in Appendix C),

$$ \phi(t', z') = \frac{M}{2\chi_t^2 a'(1+c\lambda)} \left( \ln \left( \frac{(w'(z_i) - \beta)(w'(z_i)(1+c\lambda) + c\lambda)}{(1-\beta)(1+2c\lambda)} \right) - \frac{c\lambda - \beta(1+c\lambda)}{A} \ln \left( \frac{(w'(z_i) - \beta)(1+2c\lambda)}{(w'(z_i)(1+c\lambda) + c\lambda)(1-\beta)} \right) \right) + \frac{1}{\alpha' \chi_t} \ln \left( \frac{w'(l) - \beta}{w'(z_i) - \beta} \right), \quad (36) $$

with $A = (c\lambda + \beta(1 + c\lambda))$.

## 2.4. Numerical implementation

To implement the equations derived in the previous section, we use Matlab with an ODE solver to compute the movement of the menisci in a pore doublet with and without the presence of wetting layers. We choose the ode45 function of Matlab with its default options, as it provides the best compromise of accuracy and computational speed for non-stiff problems (Dormand and Prince, 1980). This solver is based on an explicit fourth order Runge-Kutta method.

To solve these equations, we impose that the interfaces cannot move backward in the channel. This is crucial to preserve the stability of the numerical solutions (Al-Housseiny et al., 2014). To create the perturbation leading to the instability, the initial positions of the interfaces are set differently: $z_1(t=0) = 0$ but $z_2(t=0) = \tau' = 10^{-2}$. The only effect of this parameter is to accelerate or delay the onset of the instability. For a smaller $\tau$, if the system is unstable, the instability will appear

later during the invasion, whereas for a larger $\tau$, it will become visible earlier. However, changing the value of $\tau$ does not affect the overall stability of the invasion. To find the values of the wet area and the film thickness, we compute the capillary number of each channel defined as $Ca_i(t) = Ca \frac{U_i(t)A_k}{Q_p}$, where $U_i(t)$ and $Q_p$ are calculated from the time-averaged speed of the menisci in the previous time-steps. When one of the interfaces reaches the outlet (when $z'_i = 1$), the speed of the menisci is set to zero to signify the end of the invasion.

## 2.5. Stability criteria

The system stability analysis on the experiments and the simulations is done by looking at the difference in the menisci positions, $\Delta z'$, once one interface reaches the outlet. We consider three different regimes:

*   Stable: $\Delta z' < \tau'$, both interfaces arrive at the outlet almost simultaneously.

*   Strongly unstable: $\Delta z' > 0.99$, one of the two menisci remains stuck at the inlet.

*   Weakly unstable: $\tau' < \Delta z' < 0.99$, both channels are invaded, but only one of the interfaces reaches the outlet.

Alternatively, Al-Housseiny et al. (2014) derived a theoretical condition to predict whether the system is stable or unstable. In their derivation, they look at the sign of the growth rate $\frac{1}{\Delta z'} \frac{d\Delta z'}{dt'}$. On the one hand, if it is negative, any difference in the position decreases with time and the invasion is stable. On the other hand, if it is positive the perturbation grows in time until it diverges leading to an unstable invasion. For a pore doublet with rectangular channels, they show that unstable flow arises for:

$$ \frac{1-M}{1-\beta} + \frac{2\zeta\alpha \cos\theta}{Ca} - \ln \frac{w(l) - \beta}{1-\beta} \geq 0. \quad (37) $$

This relation, however, neither informs on the weakly unstable regime nor accounts for the wetting layers. This is further discussed in Section 4.

# 3. Experimental material and methods

In this section, we describe the experimental procedure to investigate the two-phase flow stability conditions in porous media. It includes the microfabrication of pore-doublet in microfluidic devices, the fluid properties, and the image analysis used to interpret the results.

## 3.1. Fabrication of microfluidic pore doublets

A photolithography method is used to create silicon wafer molds with the desired geometry. Then, PDMS (PolyDiMethylSiloxane) soft lithography is used to produce the microfluidic chips containing the pore doublet (Roman et al., 2012) (see Fig. 1). The pore-doublet length, $L$, from the bifurcation to the point where the child branches reconnect, the etching depth, $D_0$, the initial and final channel width, $w_0$ and $w_f$ (corresponding to $\alpha = (w_f - w_0)/L$) are found in Table 1. Since the produced PDMS chips have a heterogeneous hydrophobic behavior, we use either Polyvinyl alcohol (PVA) deposition or silanization techniques to get a uniformly hydrophilic or hydrophobic behavior of the channels. These two approaches are described in Trantidou et al. (2017) and Karadimitriou et al. (2013), respectively. After the PVA deposition, the micromodel becomes hydrophilic with $\theta \approx 40^\circ$ and stays hydrophilic for a long period (no change after 30 days) (Trantidou et al., 2017). The silanization process, on the other hand, creates a strongly and uniformly hydrophobic PDMS surface. Karadimitriou (2013) report a contact angle close to 180°. Using these two methods enables a wider range of viscosity ratios by choosing proper fluid pairs for which water is either the wetting or the non-wetting fluid.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Experimental images of the displacement of water from the pore doublet. The experiments are re-scaled with $t = 0$ at the beginning of the invasion and $t = 1$ at the end of the invasion. The injection occurs from the left to the right.</figcaption></figure>

<table><caption>**Table 1**
Dimensions of the micromodels used in the experiments.</caption><thead><tr><th>L (μm)</th><th>α</th><th>D<sub>0</sub> (μm)</th><th>w<sub>0</sub> (μm)</th><th>w<sub>f</sub> (μm)</th></tr></thead><tbody><tr><td>1800</td><td>0.037</td><td>45</td><td>45</td><td>112</td></tr><tr><td>3400</td><td>0.037</td><td>90</td><td>90</td><td>216</td></tr></tbody></table>

## 3.2. Fluid properties

Different fluids are picked to investigate the stability of the invasion: water with a viscosity of $1.005 \times 10^{-3}$ Pa s at room temperature, air with a viscosity of $1.85 \times 10^{-5}$ Pa s, and silicon oils with a viscosity of $5.0 \times 10^{-3}$ Pa s and $1.0 \times 10^{-2}$ Pa s. Water is used as a wetting fluid in PVA-treated micromodel and is displaced by air (viscosity ratio of $M = 1.84 \times 10^{-2}$) or by silicon oils ($M = 5$ or $M = 10$). Silicon oils are used as wetting fluids for hydrophobic micromodels and are displaced by water, leading to viscosity ratios of $M = 2.0 \times 10^{-1}$ or $M = 1.0 \times 10^{-1}$. The air-water interfacial tension is 73 mN m$^{-1}$ at 20 °C (Vargaftik et al., 1983). The interfacial tensions between the silicon oils and water were measured using a OneAttention Theta Flex tensiometer from Biolin Scientific and are both around 47 mN m$^{-1}$. During the experiments, dye (Royal Blue) is added to the water for an easier distinction of the fluids. We also performed experiments with untreated PDMS, where silicon oil is used as a wetting fluid and the non-wetting fluid is a mix of glycerol and DI water. In this case, the contact angle of the silicon oil with PDMS is close to zero. Glycerol and water can be mixed in different proportions to reach different values of viscosity, from 1.005 mPa s (pure water) to 1410 mPa s (pure glycerol) at 20 °C. The silicon oil chosen as the wetting fluid has a viscosity of 5 mPa s. The experiments in Section 4 were done with a mix containing 80 wt% glycerol and 20 wt% DI water. The mix has a viscosity of $60.1 \times 10^{-3}$ Pa s (Segur and Oberstar, 1951) and the measured interfacial tension between the mix and the silicon oil is 33 mN m$^{-1}$ at room temperature.

## 3.3. Experimental procedure and image analysis

The micromodel inlet is connected to a 4-way valve using a needle. One of the ports of the valve is closed. The two remaining ports are connected to syringes of non-wetting and wetting fluids, and A LEGATO 110 syringe pump from kdScientific is used to inject the fluids. The micromodel is placed under an inverted microscope (Eclipse Ti2-U, from Nikon), and a fast camera (HS7 high-speed camera, from Fastec) is used to visualize and record the flow.

The micromodel is first saturated with the wetting fluid, then the valve is switched and the non-wetting fluid is injected at a fixed flow rate (between $10^{-4}$ mL min$^{-1}$ and 1 mL min$^{-1}$). When the interface

approaches the pore doublet, the recording is started. After the invasion, the pore doublet is again saturated with the wetting fluid before starting the procedure again with another flow rate.

The image sequences are then processed to isolate the interfaces on each frame. Then, using an in-house Matlab code, we track the position of the three interfaces (one in the inlet channel and one in each channel) frame by frame to deduce their velocities. Using this information, we can then find the capillary number of the invasion.

# 4. Results and discussion

In this section, we investigate the role of wetting layers in pore doublets using the models introduced in Section 2 and the microfluidics experiments (Section 3). First, we present experimental evidence of stable and unstable flow regimes in rectangular pore doublets. Then, we highlight the role of wetting films in both rectangular and cylindrical channels using our numerical model coupled with experimental results in the rectangular channels case.

## 4.1. Stable, weakly unstable and unstable regimes

We conduct pore-doublet experiments that cover capillary number and viscosity ratio ranging from $Ca \in [1.2 \times 10^{-6}; 4 \times 10^{-2}]$ and $M \in [1.84 \times 10^{-2}; 1.2 \times 10^{1}]$. In all experiments, the non-wetting fluid drains the main channel and then splits into two menisci at the bifurcation to invade the two daughter channels. Typical invasions are shown in Fig. 2. We observe three different invasion patterns: (i) viscous stable invasions, if both daughter channels are invaded simultaneously at the same rate, (ii) weakly unstable drainage, if one channel is only partially invaded, and (iii) capillary unstable invasions, if only one channel is invaded by the non-wetting fluid, see Fig. 2.

Viscous stable invasion occurs if the viscosity of the injected fluid is larger than the defending fluid ($M > 1$) and if viscous forces dominate capillary forces ($Ca > 10^{-2}$). In these conditions, both channels are invaded at the same time and the menisci have the same velocity. Eventually, they re-join each other at the end of the pore doublet. All the wetting fluid is recovered. This regime, observed here in a pore doublet, corresponds to the viscous stable front propagation in porous media (Lenormand et al., 1988).

Weakly unstable invasion occurs at intermediate capillary numbers ($10^{-4} < Ca < 10^{-2}$) and low viscosity ratio ($M < 1$), for which both capillary and viscous forces are of similar intensity (Mansouri-Boroujeni et al., 2023). Viscous forces are strong enough for the injected fluid to enter both daughter branches. The displacement rate of the two menisci, however, is not equal. As a consequence, one meniscus reaches the pore-doublet outlet before the other one, stopping the progression
<figure><img src="image_5.png" /><figcaption>Fig. 3. Experimental results compared with the theoretical stability curve without films. We observe unstable experiments in the stable domain, suggesting that wetting layers impact the flow stability.</figcaption></figure>

of the delayed meniscus, and leaving residual saturation in the partially invaded channel (Dullien, 2012). This unstable regime is reminiscent of the viscous fingering introduced by Saffman and Taylor (1958) that leads to the trapping of fluid ganglia in an ensemble of adjacent pores (Lysyy et al., 2022).

Unstable invasion occurs at low capillary numbers ($Ca < 10^{-4}$) for which surface-tension forces dominate. In this case, only one of the child branches is invaded, the other remains completely saturated with the initial fluid. This instability appears because the path of least resistance is the channel with the lower capillary pressure, which is also the channel with the leading meniscus. Unlike the weakly unstable regime for which one of the branches is partially drained, the residual trapping in this regime is greater as one of the channels is still fully saturated with the wetting fluid when the meniscus reaches the outlet. This regime is the capillary instability highlighted by Lenormand and Zarcone (1985).

The experimental results are summarized in Fig. 3 using a $Ca - M$ diagram where a square represents an unstable invasion, a circle a stable one, and a triangle corresponds to a weakly unstable experiment. The theoretical frontier between stable and unstable regimes proposed by Al-Housseiny et al. (2014) is added for comparison (see Eq. (37)). Some experiments are not correctly predicted by the theoretical criteria, especially for $M \ge 10$ where we observe unstable experiments in the stable domain. We posit that this contraction between our experiments and the regime boundary is due to the impact of the wetting layers on the flow, which is not considered in the theoretical criteria. Indeed, in our experiments, we observe a wetting layer left behind the meniscus (Mansouri-Boroujeni et al., 2023). In the following sections, we study the impact of these films and how they alter the stability criteria.

## 4.2. Impact of wetting layers on two-phase flow in cylindrical channels

We look at the impact of the wetting layers in pore doublets with the cylindrical channels. The channel geometry is set to $\alpha = 0.05$, $L = 1.8 \times 10^{-3}$ m and $r_0 = 4.5 \times 10^{-5}$ m. We run two kinds of simulations: with and without the films. In Fig. 4, we compare the evolution of the meniscus positions for four ($Ca, M$) couples. In all cases, we see an impact of the wetting layer on the meniscus dynamics. This is expected as the wetting layers change both capillary and viscous forces. For small viscosity ratios ($M < 1$), the films tend to stabilize the invasion, as the difference in meniscus position predicted by the model with films is smaller than the prediction of the model without films (see Fig. 4a and b). For large viscosity ratios ($M > 1$), the opposite occurs. In Fig. 4d, the model without the film predicts a stable invasion, as the meniscus positions follow the same curve, while the model with films predicts a weakly unstable invasion. This shows the impact of the viscous coupling between the two fluids inside the channel. If the wetting fluid is less viscous than the invading fluid ($M > 1$), the invasion is more unstable, while if the wetting fluid is more viscous ($M < 1$), the invasion is stabilized by the viscous coupling.

To see the impact of the films on a wider range of parameters, we draw the $M-Ca$ phase diagram similar to Lenormand et al. (1988). We run $500 \times 500$ simulations to cover values logarithmically distributed over the ranges $M \in [10^{-3}; 10^3]$ and $Ca \in [10^{-6}; 10^{-1}]$. Each simulation is then identified as either stable, unstable or weakly unstable. The results are presented in Fig. 5. The colored regions correspond to the invasion regimes obtained using the model with films. In dashed and plain lines, we draw the boundaries between the domains without films. We observe an important difference for $M > 100$, where a horizontal asymptote bounds the stable regime including wetting layers at $Ca \approx 10^{-2}$. This is also the case for the strongly unstable domain that plateaus
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Simulations with and without wetting layers for cylindrical channels. The channel geometry is $a = 0.05$, $L = 1.8 \times 10^{-3}$ m and $r_0 = 4.5 \times 10^{-5}$ m. The two models do not always predict the same stability regime.</figcaption></figure>

at $Ca \approx 10^{-3}$. Our predictions with films for large viscosity ratio are in very good agreement with the prediction of Lenormand et al. (1988) which considered drainage in a network of pores.

For $M < 1$, we observe that the boundary of the stable regime with the film shifts slightly toward smaller $M$ and $Ca$ compared to the prediction of the model without films. This is due in part to the viscous coupling between the fluids, where the more viscous wetting layers increase the viscous dissipation, and also the reduction of the available space for the flow of the non-wetting fluid because of the presence of the wetting films, that occupy a significant portion of the channel at high $Ca$, increasing the viscous pressure drop. Both of these phenomena increase the relative strength of the viscous forces and thus stabilize the invasion.

Our results highlight the differences between the models with and without films, and in the next section, we will compare them to microfluidics experiments to see which one better describes the invasion.

### 4.3. Impact of wetting layers on two-phase flow in rectangular channels

In this section, we use the experimental data obtained in Section 4.1 to check the accuracy of the model developed in Section 2 for rectangular channels. For that, we use the Matlab code introduced in Section 3 to extract the menisci positions in time and compare them with the theoretical models, see Fig. 6. In blue and yellow are the experimental curves, and in green and red are the meniscus positions predicted by the model without the films (dashed lines) and with the films (full lines). On Fig. 6(a), the experimental curves show a weakly unstable behavior, with one meniscus lagging behind the other, as predicted by our model with films. For the same parameters, the model without the wetting layers predicts a stable invasion, which does not fit with the experiment. Fig. 6(b), the experimental curves show a strong instability, with one of the meniscus not being invaded at all,

which is coherent with the prediction of the model with the wetting layers, while the model without the films presents significant errors. This shows that our model developed in Section 2 is in good agreement with experiments and can be used to predict the invasion regimes.

We discretize the $Ca - M$ phase diagram into a grid of $500 \times 500$ points, with values logarithmically distributed over the ranges $M \in [10^{-3}; 10^3]$ and $Ca \in [10^{-6}; 10^{-1}]$ to observe the impact of the films for a wider range of parameters for channel with rectangular sections. We then run two simulations for each point: one excluding the films and one including them. Based on the position difference between the two menisci, each simulation is classified as stable (when $\Delta z_{final} < \Delta z_{initial}$), unstable (when $\Delta z > 0.99L$), or weakly unstable. We then draw the borders between each regime using the models with and without films. Finally, the experimental points are added for comparison.

As can be seen in Fig. 7, the model with films leads to an unstable regime dominating at low capillary numbers, a stable domain at high $Ca$ and $M > 1$, and finally a weakly unstable regime in between them. Similar to the classic phase diagram first introduced by Lenormand and Zarcone (1985), we recover a stable displacement for $M > 1$ and $Ca > 10^{-3}$, an unstable regime for $Ca < 10^{-4}$, associated with capillary fingering that arises when capillary forces are dominant, and the weakly unstable regime for $M < 1$ and $Ca > 10^{-4}$, dominated by the destabilizing effect of the unfavorable displacement.

Similarly to Fig. 5, the difference between the two models is higher when the viscosity ratio is far from 1. But in the rectangular case, the effects of the films are more pronounced, as more wetting fluid stays behind the meniscus, thus increasing the impact of the wetting layers, see Section 4.2. The impact of the wetting films is complex and nonlinear. On the one hand, these layers reduce the available surface area for the flow of the non-wetting fluid and change the flow condition via viscous coupling, as can be seen in Eq. (23). For $M > 1$, the viscous coupling leads to a slip-length larger than the flow surface reduction,
<figure><img src="image_7.png" /><figcaption>Fig. 5. Ca–M phase diagram of a pore doublet with cylindrical channels with $\alpha = 0.05$, $L = 1.8 \times 10^{-3}$ m and $r_0 = 4.5 \times 10^{-5}$ m. The three regimes, stable (green), weakly unstable (blue), and unstable (red) are determined using the model with the wetting layers. Without the wetting layers, the limit between stable and weakly unstable is shown in black lines, and the dotted line is the limit between weakly unstable and unstable.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Comparisons between experimental curves and the two models for, (a) silicon oil displacing water in a PVA treated pore doublet, and (b) a mix of 80 w% glycerol and 20 w% DI water displacing silicon oil. The blue and yellow curves correspond to the experimental menisci position. The dotted and continuous red and green curves correspond to the numerical prediction of the model without and with the wetting layers, respectively.</figcaption></figure>

thus decreasing the viscous drop. For $M < 1$, the slip-length is smaller than the surface reduction, increasing the viscous drop and stabilizing the invasion. On the other hand, the wetting layers also increase the curvature of the menisci, see Eq. (24), which increases the capillary force, leading to a more unstable invasion.

We see this effect for $M > 10$, where the viscous coupling and the increased menisci curvature destabilize the flow and shift the lower border of the stable regime to a higher capillary number compared to the model without films. This explains the weakly unstable experiments observed in that region, see Fig. 7. This lower border reaches a plateau at $Ca = 10^{-3}$, where increasing the viscosity ratio further does not increase the stability of the invasion, as reported by Gong et al. (2023).

For $0.3 < M < 1$ and $10^{-2} < Ca < 1$, the presence of the wetting layers also decreases the range of parameters for which the invasion is stable. Indeed, the model without the film allows stable invasion for

$M < 1$ if the capillary number is high enough, see the continuous black line in Fig. 7. For these values of $M$, the slip-length leads to a viscous coupling that compensates for the reduction in flow surface area, meaning that the viscous forces are roughly unchanged compared to the model without films, but the capillary forces are increased by the curvature increase, which destabilizes the invasion.

This complex effect of the films means that there should also be a regime for $M \ll 1$ and a high $Ca$ where the invasion would be stable. Indeed, the meniscus curvatures increases with the thickness of the film ($\Delta P_{l,cap} \propto \frac{1}{u'(z_i)-2e'_{l,r}}$) but not as quickly as the viscous drop caused by the decrease of the available surface for the flow of the non-wetting fluid ($\Delta P_{visc,II} \propto \frac{1}{(u'(z_i)-2e'_{l,av})^2}$). This means that for a capillary number high enough, we should have a stable regime despite a very unfavorable displacement. We do not observe this behavior in Fig. 7 because it
<figure><img src="image_9.png" /><figcaption>Fig. 7. Ca–M phase diagram for a pore doublet with rectangular channels with the different regimes determined using the model with wetting layer. The stability limit determined using the model without the wetting layer model is shown, as well as the experimental points. Significant differences between the two models can be seen, especially when M > 10.</figcaption></figure>

appears for $Ca > 0.1$, which is a range of capillary numbers where the hypotheses taken for the derivation of our model are no longer valid. This could partially explain the stable experiments that we observed for $M < 1$, even though some other phenomena are needed to explain these experimental results.

In Fig. 7, we observe that, while the model with films demonstrates good agreement at high viscosity ratios, some discrepancies remain in predicting certain experimental points. This could be attributed to the approximations made during the derivation of the slip length in rectangular channels (see Section 2.3) or to the assumption of steady flow, which becomes less accurate at the relatively high flow rates considered.

The different regimes predicted by the model with films exhibit significant similarities with the phase diagram obtained for complex porous structure by Lenormand et al. (1988). In their work, a stable regime is observed only when both Ca and M exceed certain thresholds ($Ca > 0.1$ and $M > 10$). Similarly, our model with films predicts a stable regime only when a, Ca, and M are above specific thresholds. In contrast, the model without the wetting layers predicts a stable regime at a low capillary number when the viscosity ratio is sufficiently high, which contradicts the conclusions of Lenormand et al. (1988).

## 5. Conclusion and perspectives

Immiscible two-phase flow in porous media remains one of the most challenging research topics in fluid mechanics. Complex flow patterns and fluid distributions arise from the intricate interplay between surface tension and viscous forces. During drainage, residual wetting layers are often deposited after the passage of the fluid-fluid interface and lie in corners and crevices. These films are ubiquitous in porous media. In this paper, we investigate their impact on flow stability using a pore-scale analysis.

Our study relies on the concept of pore doublets. This theoretical tool represents the elementary feature of a porous medium, i.e. a channel separating into two parallel daughter channels that sandwich

a solid grain and rejoin each other further on. The difference in the meniscus position in both channels informs the flow stability at larger scales. State-of-the-art pore-doublet models, however, are not able to reproduce the transition between the stable, viscous, and capillary fingering regimes observed in complex pore networks. We proposed a novel kind of dynamic pore-doublet model that accounts for the presence of the wetting film using an effective slip condition that characterizes the mean film thickness and the mutual shear stress at the interface between both fluids. The slip length is a non-linear function of the capillary number and the viscosity ratio. We proposed models for pore doublet of identical channels with circular or rectangular cross-sections.

Our model, supported by microfluidic experiments, highlights the role of the films in the flow stability, especially if the viscosity ratio is larger than one. The film acts as a lubricant for the invading non-wetting fluid. On the one hand, it decreases the viscous pressure drop if the wetting fluid has a lower viscosity. On the other hand, the viscous drop increases if the wetting fluid is more viscous than the invading fluid. This effect is amplified in rectangular geometries because the films form thicker layers in the corners inducing stronger viscous coupling between the two fluids. We obtained phase diagrams that are similar to Lenormand diagrams built using complex pore networks. In particular, we recover a similar regime transition for high viscosity ratios, which is not the case for pore-doublet models that neglect the effect of the films.

We demonstrated that despite its simplicity, pore-doublet models reconcile the complex non-linear two-phase flow behavior in heterogeneous porous media if the wetting films are considered. Because film thickness strongly depends on the pore geometry, the impact observed here in straight and smooth channels is expected to amplify in geometries more representative of porous systems.
### CRediT authorship contribution statement

**Nathan Bernard:** Writing – review & editing, Writing – original draft, Software, Methodology, Investigation, Formal analysis, Conceptualization. **Cyprien Soulaine:** Writing – review & editing, Writing – original draft, Supervision, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. **Sophie Roman:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

Experiments have been performed at Nanolab, a nano-microfluidic laboratory in Orléans (France) that is part of the ECCSEL consortium and of the project STORIES (STorage Research Infrastructure Ecosystem) - European Research Infrastructures capacities and services to address European Green Deal challenges. The research leading to these results has received support from the French Agency for Research (Agence Nationale de la Recherche, ANR) through the grant GeoMIME ANR-21-CE50-0038, the Labex Voltaire ANR-10-LABX-100-01. Cyprien Soulaine is funded by the European Research Council through the grant COCONUT (101043288) from the European Union. This work was also supported by LAAS-CNRS micro and nanotechnologies platform, a member of the Renatech french national network. This project has received financial support from the CNRS through the MITI interdisciplinary programs and Cellule Energie of CNRS through the project EsCapaDeS.

## Appendix A. P(M), Q(M) and ψ(z) functions for calculating the value of the films thickness

According to Balestra et al. (2018), for the cylindrical geometry,

$$P(M) = \frac{0.643}{2} (1 + 2^{2/3} + (2^{2/3} - 1) \tanh(1.2 \log_{10} M + 0.1657)) \quad (A.1)$$

and

$$Q(M) = \frac{1.37 M^3 + 33.84 M^2 + 111.25 M + 2.21}{M^3 + 54.50 M^2 + 44.86 M + 0.89} \quad (A.2)$$

According to Wong et al. (1995), the wet fraction in the case $Ca \to 0$ is given by,

$$\psi(z) = 1 - \frac{1 - \frac{\pi}{4}}{\epsilon(z)} \left( \frac{2\epsilon(z)}{\epsilon(z) + 1 + \sqrt{(\epsilon(z) - 1)^2 + \pi\epsilon(z)}} \right)^2 \quad (A.3)$$

## Appendix B. Solution of $\phi(z_i)$ in the cylindrical case

To solve Eq. (25), we need to integrate the first term of the right-hand side of the equation, which is non-trivial. To do this we first need to do the variable change $X = r(z) = 1 + \alpha z$, which leads us to the integration of the inverse of a fourth-order polynomial:

$$M \int_{0}^{z_{i}(t)} \frac{dz}{r(z)^{4} + 4r(z)^{3}e(M-1)} = \frac{M}{\alpha} \int_{1}^{1+\alpha z_{i}(t)} \frac{dX}{X^{4} + 4X^{3}e(M-1)} \quad (B.1)$$

Then, we must write that function as a sum of simpler terms. As we know, a polynomial P of order m whose expression is

$$P(X) = \alpha X^m + \beta X^{m-1} + \gamma X^{m-2} + \dots + \omega, \quad (B.2)$$

with $\alpha, \beta, \gamma$ and $\omega$ constants, can always be written under the form

$$P(X) = \alpha \sum_{k=1}^{n} (X - \delta_{k})^{a_{k}}, \quad (B.3)$$

where $\delta_k$ are the roots of P, n the number of root of P, and $\lambda_k$ is the order of the root $\delta_k$, which is a non null positive integer. Following this equation, we also have

$$m = \sum_{k=1}^{n} \alpha_{k}, \quad (B.4)$$

where m is the order of P. We use these properties to solve our integral. For that, we first need to re-write our polynomial, which we call Q:

$$Q(X) = X^3 (X + 4e(M - 1)). \quad (B.5)$$

### B.1. Case $e \neq 0$ and $M \neq 1$

One root of Q is zero, which is of order 3. And the other is $K = -4e(M - 1)$, so for $M \neq 1$, we can re-write Q as a sum of simpler fractions:

$$\frac{AX^2 + BX + C}{X^3} + \frac{D}{X - K} = \frac{1}{X^3 (X + 4e(M - 1))} \quad (B.6)$$

with

$$A = \frac{1}{K^3}, \quad B = \frac{-1}{K^2}, \quad C = \frac{-1}{K} \quad \text{and} \quad D = \frac{-1}{K^3} = -A.$$

And we can now integrate the term that we are working on, which gives us

$$\int_{1}^{1+\alpha z_{i}(t)} \left( \frac{AX^2 + BX + C}{X^3} + \frac{D}{X + 2\lambda} \right) dX = A \ln r(z_i) + B \left( 1 - \frac{1}{r(z_i)} \right) + \frac{C}{2} \left( 1 - \frac{1}{r(z_i)^2} \right) + D \ln \left( \left| \frac{r(z_i) - K_1}{1 - K_1} \right| \right). \quad (B.7)$$

### B.2. Case $e = 0$ or $M = 1$

In this case, $Q(X)$ has only the root of order 4, 0. In this case, the problem degenerates to the problem introduced by Al-Housseiny et al. (2014) and the integration is straightforward. So we are finally able to re-inject the previous equation in the Eq. (16), which gives us the following expression of $\phi(z)$,

$$\phi(z_i) = \frac{r'(z_i')^2}{(r'(z_i') - e_i')^2} \times \frac{1 - \frac{r'(z_i')^3}{r'(1)^3}}{3\alpha' r'(z_i')^3} - \frac{M}{8\alpha' \lambda_i^3} \left( -\ln r'(z_i') + 2\lambda_i \left( 1 - \frac{1}{r'(z_i')} \right) + 2\lambda_i^2 \left( 1 - \frac{1}{r'(z_i')} \right)^2 \right) + \ln \left( \frac{r'(z_i') + 2\lambda_i}{1 + 2\lambda_i} \right). \quad (B.8)$$

## Appendix C. Solution of $\phi(z_i)$ in the rectangular case

The Eq. (16) for a rectangular channel can be written as:

$$\phi(z_i) = \alpha \left( M \int_{0}^{z_i(t)} \frac{dz}{\chi_i^2 D^2(w'(z) - \beta)(1 + c\lambda/R_h'(z))} + \frac{1}{\chi_i} \int_{z_i(t)}^{t} \frac{dz}{D^2(w'(z) - \beta)} \right), \quad (C.1)$$

with $\beta = 6(2/\pi)^5(D_0/w_0)$. For the rest of the derivation, we set $D' = 1$ as $w_0 = D_0$ in our pore doublet. We focus on the first term of the right-hand side of the equation since it is the hardest one to integrate. We can re-write this term as:

$$\frac{1}{(w'(z) - \beta)(1 + c\lambda/R_h'(z))} = \frac{w'(z)}{(w'(z) - \beta)(w'(z)(1 + c\lambda) + c\lambda)}. \quad (C.2)$$
We now need to do a decomposition into simple elements to integrate this term,

$$ \frac{w'(z)}{(w'(z) - \beta)(w'(z)(1 + c\lambda) + c\lambda)} = \frac{1}{2\alpha(1 + c\lambda)} \frac{2(1 + c\lambda)\alpha w'(z) + c\lambda\alpha - \beta(1 + c\lambda)\alpha}{w'(z)^2(1 + c\lambda) + w'(z)(c\lambda - \beta(1 + c\lambda)) - c\beta\lambda} + \frac{\beta(1 + c\lambda) - c\lambda}{2\alpha(1 + c\lambda)} \left( \frac{\alpha}{(w'(z) - \beta)(w'(z)(1 + c\lambda) + c\lambda)} \right). \quad (C.3) $$

The second term on the right hand side can be further decomposed :

$$ \frac{1}{(w'(z) - \beta)(w'(z)(1 + c\lambda) + c\lambda)} = \frac{1}{A(w'(z) - \beta)} + \frac{-(1 + c\lambda)}{A(w'(z)(1 + c\lambda) + c\lambda)}. \quad (C.4) $$

with $A = (c\lambda + \beta(1 + c\lambda))$. We can now integrate Eq. (16), which gives us:

$$ \phi(z') = \frac{M}{2\chi_i^2\alpha'(1 + c\lambda)} \left( \ln\left(\frac{(w'(z_i) - \beta)(w'(z_i)(1 + c\lambda) + c\lambda)}{(1 - \beta)(1 + 2c\lambda)}\right) - \frac{c\lambda - \beta(1 + c\lambda)}{A} \ln\left(\frac{(w'(z_i) - \beta)(1 + 2c\lambda)}{(w'(z_i)(1 + c\lambda) + c\lambda)(1 - \beta)}\right) \right) + \frac{1}{\alpha'\chi_i} \ln\left(\frac{w'(l) - \beta}{w'(z_i) - \beta}\right). \quad (C.5) $$

## Data availability

Data will be made available on request.

## References

Aaboud, B., Saouab, A., Park, C.H., 2019. Modeling of air bubble dynamics during resin transfer molding by pore doublet model. Int. J. Adv. Manuf. Technol. 105, 2343-2355. http://dx.doi.org/10.1007/s00170-019-04435-1.

Al-Houseiny, T.T., Hernandez, J., Stone, H.A., 2014. Preferential flow penetration in a network of identical channels. Phys. Fluids 26, 042110. http://dx.doi.org/10.1063/1.4872275.

Al-Houseiny, T.T., Stone, H.A., 2013. Controlling viscous fingering in tapered hele-shaw cells. Phys. Fluids 25, Publisher: AIP Publishing.

Aussillous, G., Zhu, L., Gallaire, F., 2018. Viscous Taylor droplets in axisymmetric and planar tubes: from Bretherton's theory to empirical models. Microfluidics 22, 67. http://dx.doi.org/10.1007/s10404-018-2084-y.

Bacri, J.-C., Chaouche, M., Salin, D., 1990. Modèle simple de perméabilités relatives croisées. (modelling cross relative permeabilities). C. R. l'Acad. Sci. Sér. II 311.

Balestra, G., Zhu, L., Gallaire, F., 2018. Viscous Taylor droplets in axisymmetric and planar tubes: from Bretherton's theory to empirical models. Microfluidics 22, 67. http://dx.doi.org/10.1007/s10404-018-2084-y.

Blunt, M.J., 2017. Multiphase Flow in Permeable Media: A Pore-Scale Perspective. Cambridge University Press, http://dx.doi.org/10.1017/9781316145098.

Bretherton, F.P., 1961. The motion of long bubbles in tubes. J. Fluid Mech. 10, 166-188. http://dx.doi.org/10.1017/S0022112061000160, publisher: Cambridge University Press.

Chatzis, I., Dullien, F.A.L., 1983. Dynamic immiscible displacement mechanisms in pore doublets: Theory versus experiment. J. Colloid Interface Sci. 91, 199-222. http://dx.doi.org/10.1016/0021-9797(83)90326-0.

Chen, K., Liu, P., Wang, W., Chen, Y., Bate, B., 2023. Effects of capillary and viscous forces on two-phase fluid displacement in the microfluidic model. Energy Fuels 37, 17263-17276. http://dx.doi.org/10.1021/acs.energyfuels.3c03170, publisher: American Chemical Society.

Chen, Y.-F., Wu, D.-S., Fang, S., Hu, R., 2018. Experimental study on two-phase flow in rough fracture: Phase diagram and localized flow channel. Int. J. Heat Mass Transfer 122, 1298-1307. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.02.031.

Cottin, C., Bodiguel, H., Colin, A., 2010. Drainage in two-dimensional porous media: From capillary fingering to viscous flow. Phys. Rev. E 82, 046315. http://dx.doi.org/10.1103/PhysRevE.82.046315, publisher: American Physical Society.

Delshad, M., Pope, G.A., Sepehrnoori, K., 1996. A compositional simulator for modeling surfactant enhanced aquifer remediation, 1 formulation. J. Contam. Hydrol. 23, 303-327. http://dx.doi.org/10.1016/0169-7722(95)00106-9.

Dormand, J.R., Prince, P.J., 1980. A family of embedded Runge-Kutta formulae. J. Comput. Appl. Math. 6, 19-26. http://dx.doi.org/10.1016/0771-050X(80)90013-3.

Dullien, F.A.L., 2012. Porous Media: Fluid Transport and Pore Structure. Academic Press, Google-Books-ID: pPJCOkfZGPkcC.

Ebert, W.A., Sparrow, E.M., 1965. Slip flow in rectangular and annular ducts. J. Basic Eng. 87, 1018-1024. http://dx.doi.org/10.1115/1.3650793.

Ehrlich, R., 1993. Viscous coupling in two-phase flow in porous media and its effect on relative permeabilities. Transp. Porous Media 11, 201-218. http://dx.doi.org/10.1007/BF00614812.

Ferrari, A., Lunati, I., 2014. Inertial effects during irreversible meniscus reconfiguration in angular pores. Adv. Water Resour. 74, 1-13. http://dx.doi.org/10.1016/j.adwatwres.2014.07.009.

Gong, W., Liu, Y., Lei, W., Ju, Y., Wang, M., 2023. Viscous coupling effect on hydraulic conductance in a square capillary tube. Adv. Water Resour. 104568. http://dx.doi.org/10.1016/j.adwatwres.2023.104568.

Gredicak, V., Douat, C., Slodczyk, A., Dozias, S., Roman, S., 2025. Wettability alteration of closed glass microfluidic devices by in situ plasma. Microfluid. Nanofluidics 29, 22.

Karadimitriou, N.K., 2013. Two-Phase Flow Experimental Studies in Micro-Models. Utrecht Studies in Earth Sciences, ISBN: 9789062663262, Accepted: 2013-04-25T09:05:33Z. Publisher: UU Department of Earth Sciences.

Karadimitriou, N.K., Musterd, M., Kleingeld, P.J., Kreutzer, M.T., Hassanizadeh, S.M., Joekar-Niasar, V., 2013. On the fabrication of PDMS micromodels by rapid prototyping, and their use in two-phase flow studies. Water Resour. Res. 49, 2056-2067. http://dx.doi.org/10.1002/wrcr.20196.

Kumar, S., Foroozesh, J., Edlmann, K., Rezk, M.G., Lim, C.Y., 2020. A comprehensive review of value-added CO2 sequestration in subsurface saline aquifers. J. Nat. Gas Sci. Eng. 81, 103437. http://dx.doi.org/10.1016/j.jngse.2020.103437.

Lenormand, R., 1990. Liquids in porous media. J. Phys.: Condens. Matter. 2, SA79. http://dx.doi.org/10.1088/0953-8984/2/S/008.

Lenormand, R., Touboul, E., Zarcone, C., 1988. Numerical models and experiments on immiscible displacements in porous media. J. Fluid Mech. 189, 165-187. http://dx.doi.org/10.1017/S0022112088000953, publisher: Cambridge University Press.

Lenormand, R., Zarcone, C., 1985. Invasion percolation in an etched network: Measurement of a fractal dimension. Phys. Rev. Lett. 54, 2226-2229. http://dx.doi.org/10.1103/PhysRevLett.54.2226, publisher: American Physical Society.

Li, H., Pan, C., Miller, C.T., 2005. Pore-scale investigation of viscous coupling effects for two-phase flow in porous media. Phys. Rev. E 72, 026705. http://dx.doi.org/10.1103/PhysRevE.72.026705.

Lundström, T.S., Gustavsson, L.H., Jekabsons, N., Jakovics, A., 2008. Wetting dynamics in multiscale porous media. porous pore-doublet model, experiment and theory. AIChE J. 54, 372-380. http://dx.doi.org/10.1002/aic.11387, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/aic.11387.

Lysy, M., Ersland, G., Fernø, M., 2022. Pore-scale dynamics for underground porous media hydrogen storage. Adv. Water Resour. 163, 104167. http://dx.doi.org/10.1016/j.adwatres.2022.104167.

Magnini, M., Munichi, F., El Mellas, I., Icardi, M., 2022. Liquid film distribution around long gas bubbles propagating in rectangular capillaries. Int. J. Multiph. Flow 148, 103939. http://dx.doi.org/10.1016/j.jmultiphaseflow.2021.103939.

Mansouri-Boroujeni, M., Soulaine, C., Azaroual, M., Roman, S., 2023. How interfacial dynamics controls drainage pore-invasion patterns in porous media. Adv. Water Resour. 171, 104353. http://dx.doi.org/10.1016/j.adwatwres.2022.104353.

Moebius, F., Or, D., 2012. Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries. J. Colloid Interface Sci. 377, 406-415. http://dx.doi.org/10.1016/j.jcis.2012.03.070.

Moebius, F., Or, D., 2014. Inertial forces affect fluid front displacement dynamics in a pore-throat network model. Phys. Rev. E 90, 023019. http://dx.doi.org/10.1103/PhysRevE.90.023019, publisher: American Physical Society.

Moore, T.F., Slobod, R.L., 1955. Displacement of Oil by Water-Effect of Wettability, Rate, and Viscosity on Recovery. OnePetro, http://dx.doi.org/10.2118/502-G.

Norouzisadeh, M., Leroy, P., Soulaine, C., 2024. A lubrication model with slope-dependent disjoining pressure for modeling wettability alteration. Comput. Phys. Comm. 298, 109114. http://dx.doi.org/10.1016/j.cpc.2024.109114.

Pasquier, S., Quintard, M., Davit, Y., 2017. Modeling flow in porous media with rough surfaces: Effective slip boundary conditions and application to structured packings. Chem. Eng. Sci. 165, 131-146. http://dx.doi.org/10.1016/j.ces.2017.01.063.

Picchi, D., Battiato, I., 2018. The impact of pore-scale flow regimes on upscaling of immiscible two-phase flow in porous media. Water Resour. Res. 54, 6683-6707. http://dx.doi.org/10.1029/2018WR023172, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2018WR023172.

Protière, S., Bazant, M.Z., Weitz, D.A., Stone, H.A., 2010. Droplet breakup in flow past an obstacle: A capillary instability due to permeability variations. Europhys. Lett. 92, 54002. http://dx.doi.org/10.1209/0295-5075/92/54002.

Roman, S., Abu-Al-Saud, M.O., Tokunaga, T., Wan, J., Kovscek, A.R., Tchelepi, H.A., 2017. Measurements and simulation of liquid films during drainage displacements and snap-off in constricted capillary tubes. J. Colloid Interface Sci. 507, 279-289. http://dx.doi.org/10.1016/j.jcis.2017.07.092.

Roman, S., Lorthos, S., Duru, P., Risso, F., 2012. Velocimetry of red blood cells in microvessels by the dual-slit method: Effect of velocity gradients. Microvasc. Res. 84, 249-261, Publisher: Elsevier.

Roman, S., Soulaine, C., Kovscek, A.R., 2020. Pore-scale visualization and characterization of viscous dissipation in porous media. J. Colloid Interface Sci. 558, 269-279. http://dx.doi.org/10.1016/j.jcis.2019.09.072.
Rose, W., 2000. Myths about later-day extensions of Darcy's law. J. Pet. Sci. Eng. 26, 187-198. http://dx.doi.org/10.1016/S0920-4105(00)00033-4.

Rose, W., Witherspoon, P.A., 1956. Studies of Waterflood Performance II. Trapping Oil in a Pore Doublet. Circular no. 244.

Saffman, P., Taylor, G.I., 1958. The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous liquid. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 245, 312-329. http://dx.doi.org/10.1098/rspa.1958.0085.

Segur, J.B., Oberstar, H.E., 1951. Viscosity of glycerol and its aqueous solutions. Ind. Eng. Chem. 43, 2117-2120.

Shams, M., Raeini, A.Q., Blunt, M.J., Bijeljic, B., 2018. A study to investigate viscous coupling effects on the hydraulic conductance of fluid layers in two-phase flow at the pore level. J. Colloid Interface Sci. 522, 299-310. http://dx.doi.org/10.1016/j.jcis.2018.03.028.

Sorbie, K.S., Wu, Y.Z., McDougall, S.R., 1995. The extended Washburn equation and its application to the oil/water pore doublet problem. J. Colloid Interface Sci. 174, 289-301. http://dx.doi.org/10.1006/jcis.1995.1394.

Soulaine, C., Horgue, P., Franc, J., Quintard, M., 2014. Gas-liquid flow modeling in columns equipped with structured packing. AIChE J. 60, 3665-3674. http://dx.doi.org/10.1002/aic.14550.

Soulaine, C., Maes, J., Roman, S., 2021. Computational microfluidics for geosciences. Front. Water 3, http://dx.doi.org/10.3389/frwa.2021.643714, publisher: Frontiers.

Tabeling, P., Zocchi, G., Libchaber, A., 1987. An experimental study of the Saffman-Taylor instability. J. Fluid Mech. 177, 67-82. http://dx.doi.org/10.1017/S0022112087000867, publisher: Cambridge University Press.

Taylor, G.I., 1961. Deposition of a viscous fluid on the wall of a tube. J. Fluid Mech. 10, 161-165. http://dx.doi.org/10.1017/S0022112061000159, publisher: Cambridge University Press.

Trantidou, T., Elani, Y., Parsons, E., Ces, O., 2017. Hydrophilic surface modification of PDMS for droplet microfluidics using a simple, quick, and robust method via PVA deposition. Microsyst. Nanoeng. 3, 1-9. http://dx.doi.org/10.1038/micronano.2016.91, number: 1 Publisher: Nature Publishing Group.

Tsuji, T., Jiang, F., Christensen, K.T., 2016. Characterization of immiscible fluid displacement processes with various capillary numbers and viscosity ratios in 3d natural sandstone. Adv. Water Resour. 95, 3-15. http://dx.doi.org/10.1016/j.advwatres.2016.03.005.

Tuller, M., Or, D., 2001. Hydraulic conductivity of variably saturated porous media: Film and corner flow in angular pore space. Water Resour. Res. 37, 1257-1276. http://dx.doi.org/10.1029/2000WR900328.

Vargaftik, N.B., Volkov, B.N., Voljak, L.D., 1983. International tables of the surface tension of water. J. Phys. Chem. Ref. Data 12, 817-820. http://dx.doi.org/10.1063/1.555688.

Wong, H., Radke, C.J., Morris, S., 1995. The motion of long bubbles in polygonal capillaries. part 2. drag, fluid pressure and fluid flow. J. Fluid Mech. 292, 95-110. http://dx.doi.org/10.1017/S0022112095001455.

Yiotis, A.G., Psihogios, J., Kainourgiakis, M.E., Papaioannou, A., Stubos, A.K., 2007. A lattice boltzmann study of viscous coupling effects in immiscible two-phase flow in porous media. Colloids Surfaces A: Physicochem. Eng. Asp. 300, 35-49. http://dx.doi.org/10.1016/j.colsurfa.2006.12.045.

Zhang, C., Oostrom, M., Wietsma, T.W., Grate, J.W., Warner, M.G., 2011. Influence of viscous and capillary forces on immiscible fluid displacement: Pore-scale experimental study in a water-wet micromodel demonstrating viscous and capillary fingering. Energy Fuels 25, 3493-3505. http://dx.doi.org/10.1021/ef101732k, publisher: American Chemical Society.

Zhao, B., Alizadeh Pahlavan, A., Cueto-Felgueroso, L., Juanes, R., 2018. Forced wetting transition and bubble pinch-off in a capillary tube. Phys. Rev. Lett. 120, 084501. http://dx.doi.org/10.1103/PhysRevLett.120.084501, publisher: American Physical Society.

Zhao, B., MacMinn, C.W., Juanes, R., 2016. Wettability control on multiphase flow in patterned microfluidics. Proc. Natl. Acad. Sci. 113, 10251-10256. http://dx.doi.org/10.1073/pnas.1603387113, publisher: Proceedings of the National Academy of Sciences.