Mahdi Mansouri-Boroujeni

Cyprien Soulaine

Mohamed Azaroual

Sophie Roman

How interfacial dynamics controls drainage pore-invasion patterns in porous media

Keywords: two-phase flow, porous media, pore invasion, pore doublet, front instability, microfluidics, pore-scale investigation

How interfacial dynamics controls drainage pore

Introduction

Immiscible two-phase flow in porous media is ubiquitous and of utmost importance in many environmental and industrial phenomena, e.g. CO 2 se-questration in geological reservoirs [START_REF] Riazi | Visualisation of mechanisms involved in co2 injection and storage in hydrocarbon reservoirsand water-bearing aquifers[END_REF][START_REF] Bakhshian | Pore-scale characteristics of multiphase flow in heterogeneous porous media using the lattice boltzmann method[END_REF], soil remediation [START_REF] Philippe | Experimental study of the temperature effect on twophase flow properties in highly permeable porous media: Application to the remediation of dense non-aqueous phase liquids (dnapls) in polluted soil[END_REF], as well as enhanced oil recovery [START_REF] Amirian | Pore scale visualization of low salinity water flooding as an enhanced oil recovery method[END_REF]. In industrial applications, the two-phase Darcy model that considers a direct correlation between mobility of each phase and its saturation, is frequently used to predict the displacement process. Nevertheless, various structural and interfacial porous media properties (affinity of porous surface to one fluid, i.e., wettability, pore and throat size, pore connectivity, etc.), fluid properties (viscosity, density, interfacial tension, etc), and flow conditions result in different displacement behavior [START_REF] Peters | The Onset of Instability During Two-Phase Immiscible Displacement in Porous Media[END_REF]. In particular, when the displacing phase has a lower viscosity than the displaced phase, a ramified displacement front tends to be developed, which is a central factor in displacement efficiency. This is the case in the context of CO 2 sequestration in deep saline aquifers, where the injected supercritical CO 2 into a deep geological reservoir of porous rock has a lower viscosity than the displaced native fluid in the pore space (usually brine). Today, the different mechanisms leading to the immobilization of CO 2 in subsurface formations are well-identified. However, the effective CO 2 storage capacity worldwide is uncertain because the influence of complex coupled processes occurring during and after injection is often neglected or unknown. Therefore, the success of a secure and permanent storage of CO 2 in subsurface formations depends largely on our understanding of two-phase fluid displacements leading to residual CO 2 saturation. In particular, the development of fluid front instabilities questions the validity of the multiphase extension of Darcy's equation to predict immiscible displacement using relative permeabilities [START_REF] Gupta | An experimental study of immiscible displacement with an unfavorable mobility ratio in porous media[END_REF][START_REF] Peters | The Onset of Instability During Two-Phase Immiscible Displacement in Porous Media[END_REF].

Although the two-phase Darcy model is largely used, a number of studies have emphasized its limit [START_REF] Rose | Myths about later-day extensions of darcy's law[END_REF][START_REF] Hassanizadeh | Toward an improved description of the physics of two-phase flow[END_REF]. A key assumption in the derivation of the model by averaging Stokes equations consider that the fluid-fluid interfaces are stable [START_REF] Whitaker | Flow in porous media i: A theoretical derivation of darcy's law[END_REF]. However,recent observations report intermittent pore-invasion patterns that question this hypothesis [START_REF] Moebius | Pore scale dynamics underlying the motion of drainage fronts in porous media[END_REF][START_REF] Spurin | Intermittent fluid connectivity during two-phase flow in a heterogeneous carbonate rock[END_REF]. Moreover, experimental and numerical investigations have shown that inertial forces are not always negligible [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF][START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF]. Besides, novel experimental techniques such as micro-PIV (Particle Image Velocimetry) revealed viscous dissipation mechanisms that are not included in Darcy two-phase model [START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF][START_REF] Roman | Pore-scale visualization and characterization of viscous dissipation in porous media[END_REF].

In one of the pioneer studies of front instabilities in unsaturated porous media [START_REF] Lenormand | Numerical models and experiments on immiscible displacements in porous media[END_REF], fluid invasions were categorized as stable, viscous, and capillary flow regimes. These regimes are based on two dimensionless numbers describing active forces, the viscosity ratio M = µ (1) µ [START_REF] Bakhshian | Pore-scale characteristics of multiphase flow in heterogeneous porous media using the lattice boltzmann method[END_REF] , and the capillary number, Ca = µ (2) ϑ σ cos θ , where µ, ϑ, σ and θ represent viscosity (of the displacing (1) and displaced (2) phase), velocity, interfacial tension and contact angle. For M > 1 and low Ca the displacement of a wetting fluid by a non-wetting fluid is governed by capillary forces. The non-wetting fluid is able to enter a pore throat of radius r when the capillary pressure exceeds the entry value P e = 2σ/r. For this capillary flow regime, the pore size distribution controls the displacement, the fluid/fluid interface front advances by invading the largest pore available that have the lowest capillary resistance. This results in the onset of capillary fingering. For greater Ca, viscous force become important, the pressure differences may force the non-wetting fluid into smaller throats, and the displacement front is stable. When M < 1 and low Ca the capillary flow regime is observed, by increasing Ca viscous fingering dominates because of the lower viscosity of the invading non-wetting fluid [START_REF] Zhang | Influence of viscous and capillary forces on immiscible fluid displacement: Pore-scale experimental study in a water-wet micromodel demonstrating viscous and capillary fingering[END_REF]. Front instabilities were observed and studied thoroughly in later studies both in 3D core flood studies [START_REF] Aryana | Experiments and analysis of drainage displacement processes relevant to carbon dioxide injection[END_REF] and microfluidics experiments [START_REF] Zhang | Influence of viscous and capillary forces on immiscible fluid displacement: Pore-scale experimental study in a water-wet micromodel demonstrating viscous and capillary fingering[END_REF]. However, the large number of factors that influences two-phase immiscible flow in porous systems makes the prediction of fluid behavior difficult. Thus, pore-scale studies are needed to identify the processes at stakes and their consequences on the macroscopic front behavior. Despite recent advancements in pore-scale simulation tools, they still lack in thoroughly predicting the displacement process. In pore network models (PNM) geometry of pores and throats are replaced by regular shapes, which reduces the complexity of pore space. In addition, in both quasi-static PNM (considering single network element invasion each time) and dynamic PNM (allowing gradual and simultaneous displacements in several network elements), capillarity is considered a local phenomenon and specific to each meniscus. However, the experiments of Armstrong and Berg [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF] showed non-locality of the event and that acting forces are distributed in an interconnected pore network. Although computational barrier for computational fluid dynamics (CFD) has been almost overcome and performing high-resolution simulations has become cost-efficient, CFD methods fail to resolve physics of interfacial phenomena and describe invasions for flows dominated by capillary effects [START_REF] Soulaine | Computational microfluidics for geosciences[END_REF][START_REF] Meakin | Modeling and simulation of pore-scale multiphase fluid flow and reactive transport in fractured and porous media[END_REF]. So far, almost none of the existing models can resolve and predict accurately displacement processes at the pore-scale in porous media [START_REF] Soulaine | Computational microfluidics for geosciences[END_REF]. To predict interface behavior and have a model that incorporates interfacial phenomena, we need an in-depth understanding, accurate observation, and relevant measurement at the scale of single pores.

Although X-ray computational tomography (CT) has improved in terms of resolution and real-time characterization in recent years, it still lacks cap-turing the dynamics of small-scale fast mechanisms of invasions. Coupled with a high-resolution monitoring system, microfluidics has shown a great potential to capture key mechanisms involved in displacement processes. Especially when the displacement concerns very fast dynamics [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF][START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF]. Microfluidics devices, also called micromodels, allow direct visualization of the flow processes in rock microstructure replica confined between two parallel plates separated by a small distance [START_REF] Morais | Studying key processes related to co 2 underground storage at the pore scale using high pressure micromodels[END_REF]. Using microfluidics devices provide the ability of designing the geometry and having precise control of pore and throat shape and ratio. In addition, it brings new insight in pore-scale dynamics by using techniques such as micro-PIV.

Pore-scale experimental studies of invasion dynamics for slow displacements show that interface moves through the porous media by pinning at constrictions followed by an abrupt movement of the interface associated with pressure bursts, named after the work of Haines [START_REF] Haines | Studies in the physical properties of soil. v. the hysteresis effect in capillary properties, and the modes of moisture distribution associated therewith[END_REF] as Haines jumps [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF]. These jumps are associated with inertial effects and are suspected to rule the invasion mechanisms and change invasion patterns [START_REF] Ferrari | Inertial effects during irreversible meniscus reconfiguration in angular pores[END_REF]. Velocimetry of these jumps showed that the interface velocity can reach up to 50 times the mean interface velocity [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF][START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF]. Sudden pressure and velocity fluctuation and associated inertial movement of the interface are accompanied by relaxation and retraction of the interface inside pores, showing non-locality of capillary effects [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF]. These small-scale phenomena add to the complexity of displacement prediction and rise the question of whether it affects the global fluid displacement. Singh et al. [START_REF] Singh | The role of local instabilities in fluid invasion into permeable media[END_REF] reported a direct relation between ramification of the front and frequency of Haines jumps. Moebius and Or [START_REF] Moebius | Inertial forces affect fluid front displacement dynamics in a pore-throat network model[END_REF] reported the alteration of displacement pattern and the change in invaded pore size distributions as a result of the inertial effect of Haines jumps in their simulations. Roman et al. [START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF] and Blois et al. [START_REF] Blois | A microscopic particle image velocimetry method for studying the dynamics of immiscible liquidliquid interactions in a porous micromodel[END_REF] reported downstream perturbations and flow fluctuations during two-phase displacement process, showing local pressure gradients instabilities. Despite significant contributions in studying displacement behavior, the majority of the experimental works have focused on the overall macroscopic behavior of the interface and there is no conclusive understanding of the pore-scale invasion for precisely defined pore geometry. Besides, most of the previous works were performed for highly connected pore network systems, making it almost impossible to associate observed front disturbance to a specific invasion. Using microfluidics, we are able to design simplified porous media to study invasions at a single pore-scale and neighboring pore effects in details.

A pore doublet is a simplified consideration of a porous media, consisting of two connected channels. It was originally designed to explain the entrapment of oil in water flooding, fluid distribution, and preferential flow paths during immiscible displacement in porous networks [START_REF] Benner | Nature and importance of surface forces in production of petroleum[END_REF][START_REF] Moore | Displacement of oil by water-effect of wettability, rate, and viscosity on recovery[END_REF][START_REF] Rose | Studies of waterflood performance ii. trapping oil in a pore doublet[END_REF]. Despite its simplified and minimal replication of the porous media, it still captures the essence of the transport and interface displacement in the porous media and results can be extended to larger scales [START_REF] Moore | Displacement of oil by water-effect of wettability, rate, and viscosity on recovery[END_REF]. Theoretical consideration of the dynamics of front behavior in pore-doublet geometries has been investigated accounting for viscous and capillary forces [START_REF] Moore | Displacement of oil by water-effect of wettability, rate, and viscosity on recovery[END_REF][START_REF] Rose | Studies of waterflood performance ii. trapping oil in a pore doublet[END_REF][START_REF] Chatzis | Dynamic immiscible displacement mechanisms in pore doublets: theory versus experiment[END_REF][START_REF] Washburn | The dynamics of capillary flow[END_REF]. Moebius and Or [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF] considered effects of inertial forces on interface behavior in their model and Al-Housseiny et al. [START_REF] Al-Housseiny | Preferential flow penetration in a network of identical channels[END_REF] modeled the preferential flow paths in smoothly varying rectangular cross-section channels. Millimeter-scale experimental investigations of meniscus dynamics were carried out by Chatzis and Dullien [START_REF] Chatzis | Dynamic immiscible displacement mechanisms in pore doublets: theory versus experiment[END_REF] using pore doublets etched in glass and by Moebius and Or [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF] using glass bead micromodels. Nevertheless, direct comparison and the link between a systematic and dynamic theoretical pore-scale interface behavior and micro-scale experiments is missing. In this work, we designed a poredoublet geometry made of sinusoidal channels to improve our understanding of pore-invasion mechanisms. This work aims to investigate experimentally and numerically pore-scale dynamics, key pore-invasion mechanisms, and their contribution in defining global front behavior for different flow regimes. In this regard, we conduct experiments under drainage conditions (when the non-wetting phase displaces the wetting phase) for various fluid pairs, for unfavorable viscosity ratio at meticulously controlled flow conditions in a well-defined pore-doublet geometry. We develop a theoretical model based on Navier-Stokes volume averaging at Representative Elementary Volume (REV) that takes inertial, viscous and capillary forces into account and we evaluate invasion processes predicted for a pore-doublet model. Details on theoretical method and experimental setup are elaborated in Section 3 and Section 4 respectively. Then, experimental and numerical modeling results are presented and discussed in Section 5. We close with conclusions and perspectives.

2. Pore-doublet models: concept and geometry

The concept of pore-doublet models

Pore-doublet models are theoretical tools to investigate the fundamental mechanisms of two-phase flow in porous media [START_REF] Benner | Nature and importance of surface forces in production of petroleum[END_REF][START_REF] Moore | Displacement of oil by water-effect of wettability, rate, and viscosity on recovery[END_REF][START_REF] Rose | Studies of waterflood performance ii. trapping oil in a pore doublet[END_REF]. They conceptualize a porous medium as a pore that splits into two parallel pores -a pore doublet -rejoining each other after a distance, as illustrated in Fig. 1. The two channels are labeled i and ii, respectively. The volume in between the pore doublet represents a solid obstacle, e.g. a mineral grain. In theory, each pore can have different lengths and different cross-sections. In our study, however, we consider that the pore doublet is made of two pores of identical length L. The channel cross-section, A, is not necessarily constant, therefore it varies with the location, and we note A(x), where x is the distance along the channel centre line from the channel inlet.

The position of the interface between displacing and displaced fluids in channels i and ii are denoted h i and h ii , respectively. Therefore, the velocity of the meniscus in each channel is ḣi = dh i dt and ḣii = dh ii dt . In the pore doublet, the flow rate Q injected at the inlet is conserved throughout the system and the overall pressure drops, ∆p overall , along channel i and ii are equal. These two principles are used in Section 3 to derive the governing equations of dynamic pore-doublet models, i.e. equations ruling the evolution of h i , ḣi , h ii , and ḣii . This kind of model was originally proposed to investigate the fluid invasion behavior at the bifurcation between two pathways, and to explain entrapment of oil in water flooding in porous media [START_REF] Chatzis | Dynamic immiscible displacement mechanisms in pore doublets: theory versus experiment[END_REF][START_REF] Rose | Studies of waterflood performance ii. trapping oil in a pore doublet[END_REF]. Despite its simplified and minimal replication of the pore structure, it still captures the essence of the transport and interface displacement mechanisms in porous media and the insights obtained with pore-doublet model are used to predict field-scale results [START_REF] Moore | Displacement of oil by water-effect of wettability, rate, and viscosity on recovery[END_REF].

Pore-doublet model for microfluidic devices

To characterize fluid front displacements and improve our understanding of the interface dynamics at the pore-scale, we used pore-doublet microfluidic devices made of a succession of pores and throats (see Fig. 2) defined by a sinusoidal function that varies with the distance from the inlet (h),

r(h) = r p + r th 2 - r p -r th 2 cos 2πh H , (1)
where r p = D p /2 and r th = D th /2 are the pore and throat width, and H is the pore length. The pore profile is rectangular because of the etching process. Therefore, if D 0 is the micromodel depth, then the cross-section area is,

A(h) = 2r(h)D 0 .
(

) 2
The curvature of the meniscus at each position is

κ(h) = 2 D 0 + 1 r(h) cos(θ), (3)
where θ is the contact angle, and the permeability of the system is approximated by,

K(h) = D 2 0 12 1 -3 2 π 5 D 0 r(h) . (4)
Note that this expression corresponds to the permeability of a rectangular channel and is only valid for D 0 /(2r(h)) < 1. It is obtained from the first term of the Fourier series solution for the velocity profile of a viscous flow in a rectangular channel [START_REF] Al-Housseiny | Preferential flow penetration in a network of identical channels[END_REF].

Figure 2: Schematic of experimental setup: a pressure controller is used to control the flow in a PDMS microchip and the displacement process is monitored using a high-resolution camera mounted on a microscope and synced to a monitoring system. The schematic representation of porous media geometry of the microchip as a pore doublet with sinusoidal channels and rectangular cross-section is shown. H, D 0 , D p , D th and L respectively represent pore length, pore depth, pore width, pore throat width and channel total length.

J o u r n a l P r e -p r o o f

Journal Pre-proof 3. Theoretical framework: dynamic pore-doublet models

In this section, we develop a dynamic pore-doublet model to investigate two-phase flow mechanisms in porous media. The model is used in the following sections and confronted with microfluidic experiments.

Mass and momentum balance equations

Assuming that both phases are incompressible, the sum of the mass flow rates Q a (with a = i, ii) in each channel is equal to the mass flow rate, Q, at the pore doublet inlet -also called the total mass flow rate. We have,

Q i + Q ii = Q. ( 5
)
In each pore, at any location x along the pore channel, the flow rate is Q a (x) = va (x)A(x) (with a = i, ii) where va (x) is the section-averaged fluid velocity. Within a pore, the fluid velocity is equal to the velocity of the interface, va = dha dt = ḣa . Therefore the mass balance writes,

A(h i ) ḣi + A(h ii ) ḣii = Q, (6)
which is the first of the two equations that made the dynamic pore-doublet model.

The second equation is obtained by determining the overall pressure drop in both channels. Because of Laplace law, we have ∆p overall = ∆p (1) a (h a ) + ∆p (2) a (h a )σκ(h a , θ) with a = i, ii,

where ∆p

a (h a ) is the pressure drop in the invading fluid from the inlet of channel a = i, ii to the meniscus, and ∆p

(2) a (h a ) is the pressure drop in the displaced fluid from the meniscus to the channel outlet. The last term in the equation defines the capillary pressure at the interface position, where κ(h a ), σ, and θ represent respectively the interface curvature in each position, the interfacial tension, and the contact angle formed between the solid wall and the fluid-fluid interface. The interface curvature depends on the cross-section geometry. Therefore, if the cross-section varies along the channel length, so does the interface curvature.

The pressure drops ∆p

a (h a ) and ∆p

a (h a ) are determined by averaging the Navier-Stokes equation over a pore cross-section to form a sectionaveraged momentum equation that govern the evolution of the mean flow along the channel center line. For a Newtonian fluid with no-slip conditions at the fluid-solid interface and assuming that the flow is laminar and unidirectional along the pore center line, the integration of Navier-Stokes momentum equations following the work of Gray [START_REF] Gray | A derivation of the equations for multi-phase transport[END_REF], Whitaker [START_REF] Whitaker | Flow in porous media i: A theoretical derivation of darcy's law[END_REF] and Civan [START_REF] Civan | Porous media transport phenomena[END_REF] over a pore cross-section reads as,

ρ ∂v a (x) ∂t = - ∂ pa ∂x -µ va (x) K(x) with a = i, ii, (8)
where ρ and µ are the fluid density and viscosity, va (x) is the section-averaged flow velocity, pa is the section-averaged fluid pressure, and K(x) is the channel permeability that described the mutual friction between the wall and the fluid.

The pressure drop, ∆p

a , within the invading fluid between the channel inlet and the meniscus is therefore,

∆p (1) a = ha 0 ∂ p (1)
a ∂x dx = -ρ (1) ∂ ha 0 va (x) dx ∂t -µ (1) ha 0 va (x) K(x) dx. (9)
Notes that the flow rate, Q a , is constant within the invading fluid (from 0 to h a ). Therefore, as va (x) = Qa A(x) and Q a = va (h a )A(h a ) = ḣa A(h a ), the pressure drop within the displacing fluid becomes:

∆p (1) a = -ρ (1) ∂ ḣa A(h a ) ha 0 1 A(x) dx ∂t -µ (1) ḣa A(h a ) ha 0 1 A(x)K(x) dx. ( 10
)
Similarly, one can derive the formula for the pressure drop, ∆p

a , in the displaced phase,

∆p (2) a = -ρ (2) ∂ ḣa A(h a ) L ha 1 A(x) dx ∂t -µ (2) ḣa A(h a ) L ha 1 A(x)K(x) dx. (11)
Combining Eq. ( 7) with Eqs. ( 10)-( 11), the overall pressure drop is,

∆p overall = -ρ (2) d dt A(h a ) ḣa Θ(h a ) -µ (2) A(h a ) ḣa Φ(h a ) -σκ(h a , θ), (12) where
Θ(h) = Γ h 0 1 A(x) dx + L h 1 A(x) dx, (13)
and

Φ(h) = M h 0 1 A(x)K(x) dx + L h 1 A(x)K(x) dx. ( 14
) in which M = µ (1) µ (2) , Γ = ρ (1) ρ (2)
are respectively the viscosity and density ratios between the displacing and displaced phases.

Finally, because the pressure drops over the two connected channels in a pore doublet are equal, we have,

-ρ (2) d dt A(h i ) ḣi Θ(h i ) -µ (2) A(h i ) ḣi Φ(h i ) -σκ(h i , θ) = -ρ (2) d dt A(h ii ) ḣii Θ(h ii ) -µ (2) A(h ii ) ḣii Φ(h ii ) -σκ(h ii , θ). ( 15
)
Equations ( 6) and (15) form a system of Ordinary Differential Equations (ODE) whose unknown variables are the meniscus positions, h i and h ii , in each pore, and their displacement velocity, ḣi and ḣii . They constitute the dynamic pore-doublet model used in this paper. The parameters K(x), A(x), and κ(x) depends on the pore doublet geometry.

Dimensionless equations

To generalize the problem and investigate numerical results based on a set of common dimensionless parameters (Ca, Re, M , Γ), the mass and momentum balance equations, Eq. ( 6) and Eq. ( 15), are nondimensionalized using dimensionless variables denoted with superscript * and characteristic values denoted by subscript c: t

* = t tc , p * = p Πc , K * = K Kc , h * = h hc , κ * = κ κc , A * = A
Ac and r * = r rc . For a microfluidic device with rectangular crosssections as defined in Section 2.2, the characteristic permeability is

K c = D 2 0 12 1 -3 2 π 5 D 0 r th
, the characteristic length h c as L, the total length of the channels, and for the curvature we consider κ c = 2 D 0 + 1 r th . A c is the average of the cross-section, A c = D 0 × (r p + r th ), and r c = r th . The characteristic time-scale is t c = Achc Q (from Eq. ( 6)), and the characteristic pressure is Π c = µ (2) hcQ AcKc (from Eq. ( 12)). Then, the dimensionless mass balance equation writes

A * (h i ) ḣ * i + A * (h ii ) ḣ * ii = 1, (16)
J o u r n a l P r e -p r o o f

Journal Pre-proof and the dimensionless pressure drop is,

∆p * (overall) = -Re d dt * A * (h * a ) ḣ * a Θ * (h * a ) -A * (h * a ) ḣ * a Φ * (h * a ) -Ca -1 κ * (h * a , θ), (17)
where a = i, ii and

Θ * (h * a ) = Γ h * a 0 1 A * (x) dx + 1 h * a 1 A * (x) dx, (18)
and

Φ * (h * a ) = M h * a 0 1 A * (x)K * (x) dx + 1 h * a 1 A * (x)K * (x) dx. ( 19
)
The momentum equation, Eq. ( 17), involves four dimensionless numbers, namely the viscosity ratio M = µ (1) µ (2) , the density ratio Γ = ρ (1) ρ (2) , the Reynolds number Re = ρ (2) QKc µ (2) hcAc and the capillary number Ca = µ (2) hcQ σκcKcAc . Note that the two first dimensionless numbers depend on the fluid properties only, and that the Reynolds and capillary numbers depends on the pore geometry and flow conditions.

Finally, the dimensionless momentum the equation for a pore-doublet system writes

Re d dt * A * (h * i ) ḣ * i Θ * (h * i ) -A * (h * ii ) ḣ * ii Θ * (h * ii ) = -A * (h * i ) ḣ * i Φ * (h * i ) + A * (h * ii ) ḣ * ii Φ * (h * ii ) -Ca -1 (κ * (h * i , θ) -κ * (h * ii , θ)) . ( 20
)
The analytical solution of Θ(h * ), Φ * (h * ), and d dt * A * (h * ) ḣ * Θ * (h * ) can be found in supplementary material (8.1).

Numerical solution

The dimensionless equations Eq. ( 16) and Eq. (20) form a system of second-order nonlinear ordinary differential equations (ODE). The formula ) associated with the pore geometries used in this paper are presented in Section 8. The pore-doublet equations are solved numerically using an in-house MATLAB ® code in which the ODE system is transformed in to a first-order ODE set with mass matrix whose unknown variables are the meniscus positions and their displacement velocity. The developed model is applied to a pore doublet with sinusoidal channels and rectangular cross-section, mimicking our experimental porous system and flow conditions. Porous media dimensions are given in Table 1. A small surface area difference is imposed between channels to ensure that one of the menisci breaks through first (channel 2 is 1% wider than the other) [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF].

Experimental methodology

In this section, we present our experimental method, setup, and measurement techniques along with investigated fluid pairs and implemented model properties.

Microchip geometry and fabrication

A modified pore-doublet geometry composed of two sinusoidal channels representing successions of pores and throats in subsurface formations is considered as a simplified representation of geological porous media. A schematic representation of the porous media geometry is presented in Fig. 2, details of the designed geometry are given in Table 1. Each pore has a length of H, pore width of D p , throat width of D th , and depth of D 0 . A plasma cryogenic etching method [START_REF] Dussart | Plasma cryogenic etching of silicon: From the early days to today's advanced technologies[END_REF] is adopted to prepare silicon molds, and PDMS (Polydimethylsiloxane) soft lithography is used for micro-chip fabrication [START_REF] Roman | Velocimetry of red blood cells in microvessels by the dual-slit method: Effect of velocity gradients[END_REF]. As PDMS has a heterogeneous hydrophobic nature [START_REF] Karadimitriou | On the fabrication of pdms micromodels by rapid prototyping, and their use in two-phase flow studies[END_REF], plasma treatment, PVA (polyvinyl alcohol) [START_REF] Trantidou | Hydrophilic surface modification of pdms for droplet microfluidics using a simple, quick, and robust method via pva deposition[END_REF], and silane deposition [START_REF] Karadimitriou | On the fabrication of pdms micromodels by rapid prototyping, and their use in two-phase flow studies[END_REF] are used to obtain a homogeneous hydrophilic or hydrophobic wettability. The static contact angle and interfacial tension of the fluid pairs with PDMS have been measured using a KRÜSS ADVANCE drop shape analyzer, see Table 2. The contact angle reported are between the wetting phase and PDMS in presence of the non-wetting phase.

Fluids

Different fluid pairs are considered to obtain a range of viscosity ratios, -2 < log M < -1, representative of unfavorable displacement condition

J o u r n a l P r e -p r o o f

Journal Pre-proof for CO 2 sequestration in subsurface environments [START_REF] Nordbotten | Injection and storage of CO2 in deep saline aquifers: Analytical solution for CO2 plume evolution during injection[END_REF], see Table 2. To distinguish between different phases, the aqueous phase is dyed (PME royal blue dye). The aqueous phase is seeded with micro-particles (Polybead Carboxylate Microsphere 1 µm diameter, Polysciences) at approximately 0.06% volume concentration for micro-PIV measurements [START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF].

Fluidic setup

The experimental setup is presented in Fig. 2. A high-resolution camera (Andor Neo, sCMOS 5.50 Megapixel, up to 500 fps) is coupled with an upright Nikon Eclipse Ni microscope to capture image sequences of the displacement process. The flow is regulated using a pressure controller (OB1 MK3+, Elveflow). For each experiment, first, the wettability of the micromodel is uniformly modified to the desired contact angle, then the micro-chip is saturated with the wetting phase. Finally, the wetting phase is withdrawn from the micromodel. By adjusting the pressure gradient over the micro-chip, ranging from 30 mbar to 110 mbar, the mean flow rate varies from 2.8×10 -14 to 1.3 × 10 -11 m 3 .s -1 . All experiments are conducted at room pressure and temperature.

Image processing

Image processing techniques are applied using a self-developed MATLAB ® code to separate phases, and monitor the fluid-fluid interface behavior during the displacement process. Interface velocities are measured by tracking the interface. We also benefit from the micro-PIV technique to measure the velocity field of the wetting fluid [START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF] using the PIVlab toolkit [START_REF] Thielicke | Particle image velocimetry for matlab: Accuracy and enhanced algorithms in pivlab[END_REF].

Results and discussion

In this section, we present modeling and experimental results of drainage in a pore-doublet system. First, we identify and discuss different flow regimes. Then, we discuss the pore-invasion mechanisms governing the front instability and interface behavior during the two-phase displacement.

Flow regimes

Based on the front behavior and interface movement in our model porous medium, we identified three different flow regimes that depend on Ca and M (Fig. 3), namely viscous, capillary, and cross-over flow regime. They are observed both experimentally and numerically. To illustrate the different regimes, a spatio-temporal graph representing pore invasion where time is correlated with a color bar is plotted in Fig. 4. In the following, we provide details on the properties of these regimes.

The viscous flow regime corresponds to greater flow rates (log Ca > -0.5), the invading phase displaces the wetting phase in both channels simultaneously, see movie in supplementary materials 8.2. The interface tends to maintain the same speed in both channels with no pinning behavior of the interface at throats, as shown by the parallel superimposed displacement lines in Fig. 3b. This is in agreement with previous observations of viscousdominated flows in glass beads packs [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF]. For this flow condition, although the absolute capillary forces increase at constrictions, they are still negligible compared to viscous forces [START_REF] Lenormand | Numerical models and experiments on immiscible displacements in porous media[END_REF]. The acceleration and deceleration of the interface when it reaches and passes pore throats is observed clearly in simulations but is less pronounced in experiments Fig. 3b. In fact, the model assumes that while the interface moves through the channel, it drains all the fluid behind it. In experiments, however, at high flow rates, the interface does not displace all the fluid inside a pore and a portion of the wetting fluid remains at pore curvatures (around the region, where r = r p ), corners, and walls as wetting films, as shown in Fig. 5. As a result, in experiments, the meniscus experiences lower extent of curvature changes, hence lower capillary pressure variations, resulting in less pronounced acceleration and deceleration of the meniscus. The viscous flow regime in Fig. 4a is characterized by smooth color transitions due to no abrupt change of velocity during pore invasion. This regime is captured experimentally by Moebius and Or [START_REF] Moebius | Pore scale dynamics underlying the motion of drainage fronts in porous media[END_REF] as simultaneous invasion events of the menisci with partial pore invasion due to the rapidly moving drainage front that limits the time for pore evacuation. The residual wetting phase might increase the viscous dissipation at the fluid-fluid interface. The shear stress along the trapped wetting phase interface while the non-wetting phase is flowing in the center of the channel induces a momentum transfer at the fluid-fluid interface leading to particles recirculation in the trapped phase. In the work of Roman et al. [START_REF] Roman | Pore-scale visualization and characterization of viscous dissipation in porous media[END_REF], the momentum transfer is experimentally highlighted at the steady state condition, using micro-PIV after the non-wetting phase percolated the system. Wetting films and layers seem to play a determining role in existence of these internal circulations. In addition, we observed particles escaping the internal recirculation and moving from one pore to another through corner films. These observations are in contrast with multiphase Darcy assumptions that do not consider energy dissipation at the interface between two fluids.

The capillary flow regime occurs for the lower values of Ca (log Ca < -2.5), see movie in supplementary materials 8.3. In this regime, the invading phase only enters one of the channels. The pore-invasion mechanism in this flow condition is by sudden interfacial instability and burst to the drained pore by an abrupt decrease in local capillary pressure, i.e., Haines jumps, see Fig. 3c. It is usually assumed that a displacement of the interface in the direction opposite of the imposed flow can never arise in pore doublets [START_REF] Chatzis | Dynamic immiscible displacement mechanisms in pore doublets: theory versus experiment[END_REF][START_REF] Al-Housseiny | Preferential flow penetration in a network of identical channels[END_REF]. Nevertheless, backward displacements have been observed during drainage experiments in micromodels [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF][START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF], see Section 5.3. In our experiments, in the inlet channel the interface moves backward until the capillary pressure is balanced with the leading meniscus. This retraction of the interface is not captured in the graphs (Fig. 3c) because the retraction happens out of the camera field of view. However, the retraction of the interface in the channel that is not invaded appears well in the simulations (Fig. 3 c, retractions in channel 1). In this flow regime, the interface is pinned at each pore throat before it reaches the entry pressure of the pore. Capillary forces dominate the invasion patterns and interface invades pores that have lower entry pressure threshold. These observations are in agreement with the drainage experiments of Moebius and Or [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF] where they observed a transition from avalanche-like invasions (viscous regime) into individual pore invasions (Haines jumps) by decreasing the flow rate. For this flow regime, most of the wetting phase is displaced from each pore that is invaded, as seen in Fig. 5. The pinning-jumping behavior of the interface and the preferential flow path is shown by the sharp color transitions between pores in Fig. 4b. The cross-over flow regime is observed for intermediate capillary numbers (-2.5 < log Ca < -0.5), see movie in supplementary materials 8.4. While in previous studies, cross-over refers to the transition from viscous to capillary fingering [START_REF] Zhang | Influence of viscous and capillary forces on immiscible fluid displacement: Pore-scale experimental study in a water-wet micromodel demonstrating viscous and capillary fingering[END_REF][START_REF] Wang | Experimental study of crossover from capillary to viscous fingering for supercritical co2-water displacement in a homogeneous pore network[END_REF], in this work, it characterizes a new flow regime in terms of pore-invasion mechanisms. For this regime, the interface is able to enter both channels, but not at the same time. Each channel is drained at a different rate, see Fig. 3d. We observed a mixed behavior between pinning-jumping (i.e. Haines's jumps) and invasion (i.e. viscous behavior with no pinning at throats) in both channels, experimentally and numerically (Fig. 3 d). In experiments, however, both channels can participate in the pinning-jumping invasion whereas in simulations the invasion starts by a viscous behavior (both channels are invaded), then only one channel (the larger) is invaded by a pinning-jumping behavior (as seen on Fig. 3 d). We believe that for the microfluidic device there are small differences in capillary entry pressure from one pore to the other due to the microfabrication process, whereas for 2), although the displacing phase has invaded the pore, a portion of the wetting phase remains at curvatures of the pore. (b) Invaded pore for the cross-over flow regime, log Ca = -2.3 with fluid pair 1 (Table 2). Most of the wetting phase has been displaced, (c) Invaded pore at capillary flow regime, log Ca = -3 (related to the experiment with fluid pair 1, see Table 2). Most of the wetting phase has been displaced.

the simulation one of the channels is slightly larger on purpose. Following each jump, we observed a retraction of the interface in the channel that is not invaded at this time, as shown by small peaks on the interface position curve in Fig. 3 d,g. Further discussion on the influence of Haines jumps on the flow dynamics is provided in Section 5.3. The cross-over flow regime is characterized by sharp color transitions between pores and the invasion of both channels in Fig. 4a.

Experimental results and simulation predictions of the front behavior are presented in a log Calog M plot known as Lenormand's phase diagram [START_REF] Lenormand | Numerical models and experiments on immiscible displacements in porous media[END_REF], see Fig. 6. The transition from different flow regimes are in line with previous studies [START_REF] Lenormand | Numerical models and experiments on immiscible displacements in porous media[END_REF][START_REF] Zhang | Influence of viscous and capillary forces on immiscible fluid displacement: Pore-scale experimental study in a water-wet micromodel demonstrating viscous and capillary fingering[END_REF], however, there is a difference in transition onsets of the different flow regimes between experiments and simulations. Overall, numerical predictions of the transitions are higher in terms of Ca. The difference might be explained by the model assumptions that the meniscus always fills the whole pore cross-section, in contrary to experimental observations where a layer of wetting phase remains at walls and corners (Fig. 5), and that consider that momentum transfer at the fluid-fluid interface is negligible. Recently, it has been shown that viscous dissipation may develop at such a fluid-fluid interface [START_REF] Roman | Pore-scale visualization and characterization of viscous dissipation in porous media[END_REF], resulting in an increase of the mobility of the wetting phase during drainage. Roman et al. [START_REF] Roman | Pore-scale visualization and characterization of viscous dissipation in porous media[END_REF] suggested that viscous coupling is more significant when many small globules of the wetting fluid are immobilized in a porous medium and will impact on larger scales. Thus, this work provide new evidence that the Darcy's extended law for two-phase flows should not neglect the mutual interactions between the fluids in some cases (typically low Ca and M close to 1).

Invasion dynamics at pore throats

In this section, our aim is to identify the pore-invasion dynamics characteristic of each flow regime. For that purpose, we measured the invasion time of one pore τ (i.e. time from one pore throat to the next one) for various flow conditions, see Fig. 7. This invasion time τ is the sum of the pinning time at the pore throat t pin , and the time for the interface to move from the pore throat to the next pore throat t jump , with τ = t pin + t jump . The pinning time, during which capillary pressure increases, is indicative of the dominance of capillary forces. Whereas, a quick succession of jumps with t pin ≈ 0 is indicative of the dominance of viscous forces. The measured invasion times are normalized by the average invasion time of one pore based on the mean flow rate in the invaded channel, τ * = τ τc , with τ c = H Q i /Ac , with Q i being the average flow rate in one channel. Thus, τ * = 1 when each pore is invaded at the same rate. We plotted τ * for successive pore invasions, we present also the histograms for the values of τ * for the three flow regimes in Fig. 7. We obtain a mean value for τ * of 1.00 and a standard deviation, δ, of 0.16 for the viscous flow regime. The mean value of τ * is equal to 0.982 with δ = 0.18 for the capillary flow regime, whereas the mean value for τ * is 0.95 and the δ = 0.52 for the cross-over flow regime1 .

For the viscous flow regime there is no pinning at pore throats, thus τ = t jump , and the invasion times are very regular from one pore to the other. In this case the viscous forces are dominant and the jump velocities are dependant on bulk flow rate, as previously described by [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF]. For capillary and cross-over flow regimes, the interface is pinned at each pore throat before it invades the pore by Haines's jump. For the capillary flow regime, the duration of each jump is negligible compared with the pinning time, i.e. τ ≈ t pin . Moreover, the standard deviation of τ * from the mean invasion time is low, this is characteristic here of regular pinning times for all pores. These characteristics of invasion times were reported by Moebius and Or [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF], indeed, they observed clear waiting times between two interfacial jumps for low withdrawal rates. For the capillary flow regime, the invasion time depends on the timescale of capillary forces, i.e. on the time during which the capillary pressure rise at a pore throat. Here the invasion times are regular because the pore throat size is homogeneous through the porous media. In this work, we newly characterized the intermediate crossover regime for which invasion times are not constant through the pore invasion process. This is reflected by a wider range of τ * and higher deviation from the mean value, see Fig. 7. For this flow regime, τ = t pin + t jump , where t pin or t jump can be dominant or equally important. Consequently, the invasion time for some pores is very short (τ * < 1) , while a longer invasion time (τ * > 1) is observed for other pores. This is attributed to the interaction between the different active forces in the invasion process. For this flow regime, the velocity of the interface reached during a pore invasion may initiate a quick invasion of the next pore (τ * < 1), as observed for the viscous flow regime. On the contrary, as the two menisci advance, one event (i.e. one jump) in one channel dominates and can inhibit invasion in the other channel, as observed for the capillary flow regime and by Armstrong and Berg [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF]. Thus, this flow regime is characterized by cooperative pore filling events that will be further discussed in the next section using micro-PIV measurements.

The pinning behavior, i.e. when the interface stops at a pore throat, was observed experimentally. However, in our simulations the interface is always moving with no actual pinning time. We can see in Fig. 3c that for the capillary flow regime, at pore throats the interface creeps slowly toward the tightest part of the channels. We believe that this is due to the imposed flow rate in the model and the effect of dynamics of corner and film flow during experiments (not permitted in the model) for this flow regime. Indeed, when the interface is pinned at a pore throat, we observed a flow of wetting phase through corners of the channels in the direction of the imposed flow rate. The dynamics of corner flows and their consequences on pore-invasion mechanisms will be part of a future work.

In this section, we demonstrate the role of interfacial dynamics on defining pore-invasion mechanisms. These interfacial effects are further discussed in the next section.

Influence of Haines jumps on flow dynamics in neighboring pores

In this section our aim is to determine what are the consequences of Haines jumps on the macroscopic transport properties of multiphase flow in porous media (e.g. fluid entrapment, energy dissipation). For that, we focus on the interface and fluid dynamics for slow drainage conditions (i.e., capillary and cross-over flow regime) when Haines jumps is the dominant pore-invasion mechanism. First, we discuss the dynamics of the fluid-fluid interfaces during invasion, then we present results on the dynamics of the fluids using Particle Image Velocimetry measurements.

Fluid-fluid interface dynamics during Haines jumps

The dynamics of the interface for the cross-over flow regime is shown on Fig. 8. In this case, channel 1 is invaded by Haines jumps while the interface is pinned at the pore throat in channel 2. On Fig. 8b we show measurements of the position of the interface as a function of time in channel 2. In the channel that is not invaded at that time (channel 2), the interface shows a back and forth movement that is synced with Haines jumps happening in the adjacent channel (channel 1). The sequence of events in channel 1 and 2 is as follows, (1) the fluid-fluid interfaces are pinned at a pore throat in both channels while the capillary pressure increases, (2) the pore entry pressure is reached for channel 1 leading to the invasion of the pore by a Haines jump, the capillary pressure drops and the interface in channel 2 retreats toward the widest part of the constriction, (3) the interface advance until the narrowest part of the pore throat in both channels, and the capillary pressure increases again until a new invasion either in channel 1 or 2.

The oscillations of the interface after abrupt invasions have been observed previously in microfluidic and glass beads experiments [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF][START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF][START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF]. In this work, we are able to link each fast event (Haines jump) to its associated backward movement in the adjacent channel due to pressure release. In the micromodel experiments of Armstrong and Berg [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF] the interface is retracting to the nearest pore body (i.e. nearest larger cross-section) after an invasion by Haines jump. In our observations, the interface is also retracting to larger sections of the channel, but the retraction is limited to the pore throat vicinity, see Fig. 8. Armstrong and Berg [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF] observed meniscus retractions over a distance of up to 4 pores from the pore drainage events. In this work, we observe that the oscillations of the interface can arise very far from the invaded pore (up to the full channel length that contains multiple pores). Thus, the capillary pressure difference during drainage in porous media has probably effects over many pores around the invasion events. This should be kept in mind when defining a two-phase representative elementary volume [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF] for upscaling purpose. It is also worth mentioning that the interfacial oscillations show a size-dependant property; the larger the pore sizes are, the more dominant the inertial oscillations will be. That is why in experiments of Moebius and Or [START_REF] Moebius | Interfacial jumps and pressure bursts during fluid displacement in interacting irregular capillaries[END_REF] (mm scale) the repetitive inertial oscillation is clearly distinguishable. However, in our geometry and work of Armstrong and Berg [START_REF] Armstrong | Interfacial velocities and capillary pressure gradients during haines jumps[END_REF] (µm scale) the oscillations are rapidly damped. This size dependency of interfacial dynamics is also observed in the model. When larger systems show the full oscillation wave, while in small systems the interface rapidly finds its equilibrium state.

Several studies have also observed that for capillary-dominated regimes, two-phase flow occurs along pathways that periodically disconnect and reconnect in a series of snap-off and reconnection events called intermittent pathway flows [START_REF] Spurin | Intermittent fluid connectivity during two-phase flow in a heterogeneous carbonate rock[END_REF]. We believe that oscillations of the interface due to Haines jumps are one mechanisms responsible for intermittent pathway flows.

As mentioned previously, a backward movement of the interface is usually not permitted in pore-doublet models [START_REF] Chatzis | Dynamic immiscible displacement mechanisms in pore doublets: theory versus experiment[END_REF][START_REF] Al-Housseiny | Preferential flow penetration in a network of identical channels[END_REF]. Our observations clearly indicate the need to allow such displacements in the models.

Fluid dynamics at slow drainage conditions

In this section, our aim is to identify the consequences of interface dynamics (i.e. Haines jumps) on fluid flow , through the porous media. We performed micro-PIV measurements in the aqueous phase to relate interface dynamics to fluid flow in our system. In Fig. 9 we report measurements of fluid velocities using micro-PIV at different locations in the channels, as well as the velocity of the interfaces, for the capillary and crossover flow regime respectively2 . In Fig. 9 we show that front velocity during Haines jumps reaches velocities much higher than the average front velocity, up to 34 times greater for the capillary flow regime, and 15 times greater for the crossover flow regime here. , the fluid velocity is disturbed and shows an oscillating movement linked with each jump of the interface (see zooms in Fig. 9 b andc). This back and forth movement is observed close to the invading interface, as well as far from it (more than 15 pores away). This back and forth movement of the wetting fluid has been previously observed during drainage in micromodels by Roman et al. [START_REF] Roman | Particle velocimetry analysis of immiscible two-phase flow in micromodels[END_REF] and Blois et al. [START_REF] Blois | A microscopic particle image velocimetry method for studying the dynamics of immiscible liquidliquid interactions in a porous micromodel[END_REF]. It has been mentioned that it is probably a consequences of Haines jumps. However, due to the complex pore structure used in these studies it was not possible to relate interface dynamics events to fluid flow behavior. Here we provide the confirmation that the back and forth movements of the fluid are a result of pressure instabilities caused by Haines jumps. We also show that these pressure instabilities can propagate through the system more than ten pores apart from the drainage events.

We identified a major difference between the capillary and crossover flow regime. For the capillary flow regime, the flow is disturbed by each Haines jumps with no effect on the interface dynamics, indeed, the jump dynamics is dominated by capillary forces. This results in jumps regularly spaced in time and invading only one capillary. For the crossover flow regime, the flow is more disturbed after each jump with consequences on the interface dynamics. A jump of the interface comes with a backward movement of the fluid. This backward flow has the same intensity for the capillary and the crossover flow regime, around -1 times the mean velocity, see Fig. 9 b,c. The backward movement is followed by an acceleration of the forward flow, that reaches velocity around 1 times the mean velocity for the capillary flow regime, and up to 4 times greater than the mean velocity for the crossover flow regime. This sudden acceleration of the flow is sometimes linked with avalanche-like pore invasion for the crossover flow regime. Thus, by increasing the flow rate and moving from the capillary to the crossover flow regime, it causes interference between the dynamics of pore invasions in both channels, and cooperative pore filling is observed. Indeed, invasion in one channel might provoke or delay invasion in the adjacent channel and cause different poreinvasion dynamics from pore to pore.

Thus, it seems that flow perturbations throughout the whole porous medium are driven by local pressure gradients that can change in intensity and direction. In this work we provide a clear link between each interface dynamics and flow instabilities. Flow perturbations have consequences on the interface dynamics as we increase the flow rate and move from the capillary to the crossover flow regime. These phenomena emphasize the importance of taking into account pore-scale mechanisms related to the drainage pressure J o u r n a l P r e -p r o o f Journal Pre-proof [START_REF] Lenormand | Mechanisms of the displacement of one fluid by another in a network of capillary ducts[END_REF] for the crossover flow regime.

Conclusion

In this work, the multiphase displacement front instability in porous media for unfavorable viscosity ratio (M < 1) at drainage conditions in a pore doublet was revisited. A microfluidic system was used to study fluids and meniscus displacement in PDMS microchips with controlled pore and throat shape and size. A theoretical model based on the volume-averaged Navier-Stokes equation was proposed to study the interface behavior during the displacement process. A micron-scale pore-doublet porous system composed of a succession of pores and throats in form of sinusoidal channels was used to investigate the front behavior for various fluid pairs and flow conditions. With this setting, having a meticulous observation and measurement of displacement properties and pore invasions was done and compared directly with the model. We have captured and modelled fundamental pore-invasion mechanisms associated with each flow regime and we have highlighted the role of viscous and capillary forces on the dynamics of interface and pore invasion.

Major differences were observed between viscous and capillary-dominated flow regimes in terms of front behavior, interface dynamics, and pore-invasion mechanisms. For the viscous flow regime, viscous forces are dominant and control the front and displacement. Displacement occurs in both channels with a stable front with continuous interface displacement and it does not stop at throats due to its velocity and momentum forces. Meanwhile, the interface does not find enough time to expand over the whole pore surface area and a body of residual wetting fluid remains at curvatures of pore space, normal to the invasion direction. The retention of wetting layer in corners and curvatures of the geometry for viscous flow regime results in a greater residual saturation [START_REF] Hoogland | Drainage dynamics controlled by corner flow: Application of the foam drainage equation[END_REF][START_REF] Naghavi | An experimental study of the factors influencing apparent wetting phase residual saturation in dynamically drained porous media[END_REF]. Haines jumps are the dominant pore-invasion mechanism at the capillary flow regime. For the capillary flow regime, capillary forces are dominant and the pore in which the interface reaches the invasion threshold pressure will be invaded. This leads to suppression of the invasion in one of the channels. For this pore-invasion mechanism, the interface expands laterally and displaces all the wetting phase, except for corner wetting fluid. The crossover flow regime newly identified presents a mixed behavior of capillary and viscous flow regime. Both channels are invaded; however, they will not J o u r n a l P r e -p r o o f Journal Pre-proof participate similarly in the displacement. Pore-by-pore invasion mechanism is the main invasion process. By focusing on pore-scale dynamics, we showed that the flow is highly disturbed by each Haines jump for the cross-over flow regime. It results in pore invasions that can be provoked or prevented, and variable duration between each jumps.

The presented model is able to capture main pore-invasion mechanisms and show similar displacement behaviors for all flow regimes. Results of the model are used to predict flow regime transitions in terms of log Calog M phase diagram. There is still room for improvement of the model by taking into account viscous slip between displacing and displaced phase and including corner flows, that will be a part of the following work.

Acknowledgement

Experiments have been performed at Nanoµlab, a nano-microfluidic laboratory in Orléans (France) that is part of the ECCSEL consortium and of the project StoRIES (Storage Research Insfrastructure Eco-System) -European Research Infrastructures capacities and services to address European Green Deal challenges. MMB acknowledges the support from French Geological Survey (Bureau de Recherches Géologiques et Minières, BRGM) for PhD grant. The research leading to these results has received support from the French Agency for Research (Agence Nationale de la Recherche, ANR) through the grant Planex ANR-11-EQPX-36, the Labex Voltaire ANR-10-LABX-100-01, and through the FraMatI project under contract ANR-19-CE05-0002. Authors show their sincere gratitude toward Dr. Arnaud Stolz, Dr. Rémi Dussart and Jack Nos for giving us the opportunity to prepare microfluidics master molds in the cleanroom of GREMI (UMR7344). 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». MMB thanks Dr. Pierre-Yves Passaggia for lending a phantom v341 fast camera for image sequence acquisitions.

Figure 1 :

 1 Figure 1: A schematic representation of a general pore-doublet geometry. A constant flow rate of Q is injected to the system. L, h and A(h) respectively represent total length of the pore doublet, invaded length in each channel and the cross-section area of the channel as a function of the position.

  phase log M θ ( • ) σ(mN.m -1 )

Figure 3 :

 3 Figure 3: (a) A snapshot of the porous media and fluid displacement for a cross-over flow regime. Dynamics of interface position in the pore doublet for different viscosity ratios (log M = -1.7, -1.3, -0.56) and capillary numbers (log Ca = -4 to -0.8). Dimensionless interface positions at each channel with respect to dimensionless time for (b) viscous flow regime (log (M, Ca) = (-1.7, -0.78), (-1.3, -1.4) and (-0.56, -1.2)) experimentally and numerically (-1.7, -0.27), (c) capillary flow regime (log (M, Ca) = (-1.7, -3.97), (-1.3, -2.82) and (-0.56, -3.47)) experimentally and numerically (-1.7, -2.87) and (d) crossover flow regime (log (M, Ca) = (-1.7, -2.3), (-1.3, -2.77) and (-0.56, -2.3)) experimentally and numerically (-1.7, -0.56). The experimental snapshots on each figure shows the position of the interfaces for M = -1.7 at the time that intersects with the curves.

Figure 4 :

 4 Figure 4: Spatiotemporal graph of the interface position in pore-doublet geometry at (a) viscous flow regime, (b) capillary flow regime, and (c) crossover flow regime. In this graph the interface position at each time is correlated with a color as indicated by the color bar and the image show progression of the interface with respect to time.

Figure 5 :

 5 Figure 5: Images of an invaded pore for the different flow regimes. (a) Viscous flow regime, log Ca = -1.27 and fluid pair 3 (Table2), although the displacing phase has invaded the pore, a portion of the wetting phase remains at curvatures of the pore. (b) Invaded pore for the cross-over flow regime, log Ca = -2.3 with fluid pair 1 (Table2). Most of the wetting phase has been displaced, (c) Invaded pore at capillary flow regime, log Ca = -3 (related to the experiment with fluid pair 1, see Table2). Most of the wetting phase has been displaced.

Figure 6 :

 6 Figure 6: The phase diagram in log(Ca)log(M ) coordinate based on simulation results along with our experimental results and results of Zhang et al. [17].

Figure 7 :

 7 Figure 7: Normalized invasion times for various flow conditions as a function of the number of invaded pores. For capillary flow regime and viscous flow regime, invasion times are almost constant, as shown with the sharp histogram peak. However, for the crossover flow regime, invasion time varies from pore to pore and has a wider range of values.

Figure 8 :

 8 Figure 8: a) Snapshots of the back and forth movement of the interface at a pore throat in channel 2 while the interface is advancing in channel 1. b) Measurement of the dimensionless position of the interface in channel 2 as a function of the dimensionless time. Each oscillation of the interface at the pore throat in channel 2 is synced with a Haines jumps in channel 1. log Ca = -2.08 and log M = -1.7, corresponding to cross-over flow regime.

Figure 9 :

 9 Figure 9: a) Snapshot of a drainage displacement showing the positions where micro-PIV measurements are reported (color dots, the same colors are used for the curves in b and c)

Q

  Nomenclature pa section-averaged fluid pressure (Pa) vi (x) section-averaged fluid velocity (m.s -1 ) ∆p (1) i (h i ) pressure drop in the invading fluid from the inlet of channel i to the meniscus (Pa) ∆p (2) i (h i ) pressure drop in the displaced fluid from the meniscus to the outlet of channel i (Pa) ∆p overall overall pressure drop in both channels (Pa) ḣi velocity of the meniscus in channel i (m.s -1 ) Γ density ratio κ(h) curvature of the meniscus at each position (m -1 ) κ * dimensionless curvature of the meniscus µ (1) viscosity of the invading phase (Pa.s) µ (2) viscosity of the displaced phase (Pa.s) ρ fluid density (kg.m 3 ) σ interfacial tension (N.m -1 ) τ invasion time of one pore (s) τ * dimensionless invasion time of one pore θ static contact angle ( • ) ϑ velocity of the fluid-fluid front (m.s -1 ) A(h) channel cross-section at position h (m 2 ) flow rate (m 3 .s -1 ) r(h) radius of the channel at position h (m) r * dimensionless radius of the channel r p half pore body width (m)

Table 1 :

 1 Dimensions of the designed porous geometry.

	Pore size (D p ) Throat size (D th ) Depth (D 0 ) Pore length (H) 160 µm 40 µm 40 µm 160 µm

for Θ * (h * ), Φ * (h * ) and d dt * A * (h * ) ḣ * Θ * (h *

Table 2 :

 2 Fluid properties.

J o u r n a l P r e -p r o o f

The first pore invasions (before the vertical line) are under the inlet effect, as seen by first pinning times that are lower than the mean values, these first pores are not taken into account for the averaged values

Note that PIV measurements are reported with an Eulerian representation (i.e. velocity of the fluid viewed from a fixed point in space), whereas interface velocity is reported in a Lagrangian way (i.e. by following the interface)

Declaratio if ioterettt

The authors declare that they have no known competng fnancial interests or personal relatonships ☒ that could have appeared to infuence the work reported in this paper.

The authors declare the following fnancial interests/personal relatonships which may be considered ☐ as potental competng interests: