Bubble removal in microfluidic channels surrounded by gas-permeable media: experiments and a predictive model

Abstract

Controlling the removal of bubbles from channels is crucial in microfluidics, either to eliminate air pockets if they are unwanted, or in pumpless microfluidic applications where receding bubbles is a way to induce liquid flows. To provide a better physical understanding of air removal in microchannels, we study the dynamics of invasion of wetting liquids in dead-end microchannels surrounded by an air-permeable medium. Using polydimethylsiloxane (PDMS)-based devices, we demonstrate that gas permeation through the channel walls drives an exponential decay in trapped air length with time (in marked contrast with the so-called Lucas–Washburn law of imbibition in porous media), providing a straightforward route to bubble elimination. Systematic experiments varying channel width, height, and PDMS thickness reveal how geometric and material factors modulate the refilling timescale. A simple analytical model, coupling capillarity and gas diffusion, captures these results quantitatively. For this purpose, we introduce an explicit expression for the interfacial curvature in microchannels with heterogeneous wettability (e.g., PDMS-on-glass). This framework offers practical guidelines for microfluidic engineers aiming to prevent or remove trapped bubbles without relying on active pumping.

---

1 Introduction

Spontaneous or forced imbibition in porous media plays a central role in many natural and industrial processes, such as the infiltration of water in soils or the recovery of hydrocarbons in oil reservoirs. In idealised microfluidic channels open to an external air reservoir, the imbibition dynamics can often be described by extensions of the Lucas–Washburn law.^1–5 Under these conditions, the principal limitation typically arises from viscous dissipation within the invading liquid phase,^1–3 leading to a characteristic growth of the liquid front proportional to the square root of time. Imbibition can also occur through mechanisms driven by thermodynamics at imposed pressure, for instance via condensation in ink-bottle pores. These nanoporous constrictions lead to liquids existing under negative pressures, as highlighted by Vincent et al.⁶ in a context relevant for xylem embolism refilling in plant physiology.⁷ Understanding such phenomena is essential for explaining how cavitated conduits in plants might refill despite internal pressure deficits.

Microfluidic systems present additional motivations for studying liquid imbibition, particularly when undesired air bubbles enter a chip. Such bubbles can result from incomplete priming, dissolved gas emerging under pressure/temperature shifts, or imperfect sealing of the microchannel network. Once formed, they can severely disrupt flow patterns, obstruct the transport of important reagents or cells, and even cause oxidative or chemical reactions if retained in oxygen-sensitive processes. Consequently, there is strong incentive to eliminate these trapped air pockets, which connects directly to the idea of liquid imbibition removing gas from otherwise inaccessible parts of the device.

Kang et al.⁸ studied the elimination of spurious bubbles in microfluidic channels, by permeation through PDMS walls under pressurisation. Their experiments revealed an exponential decay in bubble volume, which is faster when the overpressure in the channels increases. Another strategy to eliminate bubbles in microchannels has been proposed by Guo et al.,⁹ using three parallel channels separated by constrictions. Interestingly, they mention that their device is bioinspired by the repair of embolisms in angiosperms,^7,10–12 although this biological question is still hotly debated and that the possible repair mechanisms in plants remain unclear. They recover an exponential decay of the bubbles, and show experimentally that the decay time increases with the channel width and the PDMS thickness.

Xu et al.¹³ provided a broader review of air transport through PDMS for applications, mentioning that air can be eliminated from the channels either by first degassing the PDMS itself, or by applying a low-pressure environment in separate, adjacent chambers, thus drawing the gas out from the channels of interest. Not only can this be used to eliminate air bubbles, but it opens the way to vacuum-driven power-free pumpless microfluidics, with applications e.g. for blood separation or cell cultures. This possibility to operate microfluidic devices for sample preparation without complicated and bulky pumping has developed recently a lot of interest.^14–18 A comprehensive overview of pervaporation-based strategies and other mass-transfer-driven microfluidics has also been provided by Bacchin et al.¹⁹

However, the dynamics of liquid refilling (or bubble removal) in PDMS microchannels has hitherto not been rationalised by a comprehensive model, whilst it is necessary to control the timescale associated to e.g. pumpless microfluidics. In particular, the dynamics is often quantified in terms of permeability through a thin membrane of thickness δ.^8,9,13 However, in most practical situations in microfluidics, the PDMS surrounding the channels cannot be identified to a thin membrane, because of side effects (fluxes from the side walls of the channel) and because the top PDMS part of the channel must generally not be too compliant. To account for these complications (as compared to the simple membrane case), empirical correction factors have been proposed,⁸ but they have not been related to the device geometry by a predictive model. This calls for a proper study of the influence of all geometrical parameters: not only PDMS thickness, but also channel width and height, on the flux.

In this paper, we revisit the problem of bubble removal and liquid imbibition in PDMS microchannels, focusing on dead-end geometries in which air cannot be simply swept away by a primary flow. We present experiments under well-controlled capillary overpressure conditions, systematically varying channel width, height, and PDMS thickness. We then propose a theoretical framework that accounts for both capillary and imposed pressures, combined with mass transport through the entire PDMS bulk. This allows us to quantify precisely the timescale of bubble shrinkage and channel refilling, moving beyond the previously reported exponential fitting. Our objective is to provide microfluidic users with a simple, yet robust predictive tool for controlling bubble dynamics in permeable environments, thereby improving the reliability of pumpless microfluidic technologies and other air-permeable-medium-based applications.

2 Materials and methods

2.1 Microchannel production

A silicon wafer serves as the base for photolithography using an SU-8 photoresist (Gersteltec). After depositing a uniform layer of resist by spin coating (to achieve a target thickness h), the wafer is selectively exposed through a photomask (Selba). The unexposed areas are subsequently removed in propylene glycol methyl ether acetate (PGMEA, Sigma-Aldrich), revealing a topography corresponding to the final channel geometry. A post-exposure bake at 200 °C completes the mold fabrication.

The polydimethylsiloxane, a silicone elastomer hereafter referred to as PDMS, (Sylgard™ 184, Dow) is prepared by thoroughly mixing a 10 : 1 mass ratio of base to curing agent. The mixture undergoes vacuum degassing for approximately one hour to eliminate trapped air. A controlled volume (e.g., 2.5 mL) of this degassed solution is then spread over the structured wafer by spin coating for 10 s at 500 rpm, followed by 40 s at speeds between 700 and 2000 rpm, depending on the desired PDMS thickness H. The coated wafer is then placed for an additional hour to ensure uniform levelling of the film, before the curing step.

The spin-coated PDMS is then cured in an oven at 65 °C for 24 h. Once polymerised, the flexible PDMS layer is gently released from the silicon wafer, exposing the channel features in negative topography. A brief atmospheric plasma treatment on a glass slide provides a surface amenable to bonding. Finally, the PDMS is carefully placed onto the plasma-treated slide.

2.2 Microchannel design

We investigate two channel geometries: first channels of constant width, and then channels of variable width. The first geometry consists of straight channels of length L₀ with a constant width w, examined in sec. 3 and modelled in sec. 4. The height of the channels is h, and the thickness of the total PDMS membrane is H, such that a layer of PDMS of thickness δ = H − h separates the top wall of the channels from the surrounding air (Fig. 1). Five such channels, with similar height and PDMS thickness but different widths, are actually placed as an array on the same chip (Fig. 2a). We checked that they are far enough apart that they can be considered as independent as far as the flux through PDMS is concerned.^20,21 The second geometry builds on insights from the first configuration and is developed in sec. 6: it consists in channels with a width w(x) that varies linearly along their longitudinal direction x (Fig. 2b). This extended design enables further tests of our theoretical approach under spatially non-uniform geometries, and promises reliable applications on networks with arbitrary geometries.

Fig. 1 Sketches of the channel, not to scale. (a) Sketch of the cross-section, with the relevant dimensions and the coordinates. (b) Longitudinal side-view sketch along the channel. The meniscus, sketched as an arc of circle inside the channel, separates a gas pocket on its right side and a liquid-filled region on its left side. The green dashed arrows depict the air flux leaving the channel by permeation through PDMS, owing to the overpressure in the gas pocket. (c) Top-view sketch of the channel. In (b) and (c), the drop deposited at the entrance of the channel is sketched as a portion of a circle.

Fig. 2 Snapshots of the channels used in this study. Scale bars correspond to 1 mm. (a) Single channels of constant width w: from bottom to top, w = 60, 80, 100, 120 and 140 μm. The inlet of the channels is at the left of the picture, where menisci can be observed (red dashed arrows). The dead ends of the channels are the right of the picture, indicated by blue solid arrows. (b) Two channels of linearly variable width, decreasing at left and increasing at right. They are connected at the beginning of the experiment, but once the meniscus (at the bottom of the snapshot, indicated by a red dashed line) reaches the bifurcation, it splits in two menisci, one in each channel, which then evolve independently towards the two dead ends (in the top, indicated by blue solid arrows).

2.3 Experimental protocol

The experiments proceed as follows. First, the channels are unsealed at one end by incising the PDMS layer with a sharp scalpel. The incised channels are placed under a binocular. We first checked by visual inspection that the incision is sharp and clean; indeed, uncontrolled roughness is to be avoided, because it is known to affect flow resistance²² and transfer rates, e.g. in the case of evaporation.²³ A drop of liquid of volume of order 1 mL is deposited at the open end of the channels. Imbibition ensues, and is recorded with a CCD camera. In the images, the advancing menisci separating dark bubbles from bright liquid-filled portions of the channels are readily detectable: they are highlighted by the red dashed arrows in Fig. 2. It is worth noting that the liquid-filled parts of the channels are almost invisible, owing to the very similar values of the refractive index of PDMS (n = 1.43) and of pentanol (n = 1.41), the liquid used in most of our experiments (see sec. 2.4). Once imbibition is completed and all bubbles have disappeared, the movies are analysed using home-made Matlab scripts to extract the time evolution of the bubble length L(t). All experiments are run at 25 ± 1 °C.

2.4 On the choice of the liquid

The choice of the liquid is important in this study. We first checked that it is impossible to fill the microchannels with pure water. By placing a water droplet at the channel entrance, we then observed that the meniscus remains stuck at the entrance. Indeed, water on PDMS has a very high contact angle (θ = 108°),²⁴ hence PDMS is too hydrophobic for spontaneous imbibition with pure water to take place. However, it sufficed to add a minute amount of soap in the water droplet to initiate filling: soap decreases the surface tension of water, hence its contact angle with PDMS, enough for spontaneous imbibition to start.

Since we wanted to avoid the potential complications associated with solutions, this early observation led us to search for pure liquids that present good wetting with PDMS, hence liquids of low surface tension. The wetting liquid should be not too soluble in PDMS to minimise swelling.²⁵ Lastly, it must be not too volatile, to keep a little reservoir of liquid at the entrance of the channel throughout the experiments. The alcohol family presents some good candidates, and we chose pentan-1-ol (henceforth called pentanol for simplicity) to produce all of the experimental data and results presented in this work, except in sec. 3.4 where a comparison with octanol is made. We note for future reference that the contact angle of pentanol on PDMS is around θ = 30°.²⁶

3 Experimental results

3.1 An exponentially decaying dynamics

We show in Fig. 3a the time evolution of the length of the air pocket L(t) for five channels of fixed total length L₀, height h and upper wall thickness δ, but varying width w. All curves display a decrease of the air length with time, with a rate of decrease |L̇| = −dL/dt which itself decreases with time. The dynamics is faster for narrower channels.

Fig. 3 (a) Time evolution of the length of the air-filled part of five channels, with dimensions L₀ = 18 mm, h = 60 μm, H = 66 μm, and w = 60 (○), 80 (□), 100 (◊), 120 (△) and 140 μm (▽). Curves are fits by an exponential law L(t) = L₀e^−t/τ. Inset: plot of the velocity of the air/liquid interface as a function of the length of the air-filled part of the five channels. Lines are fits by a linear law |L̇| = L/τ_L. (b) Five independent time evolutions (with different colours) of the length of the air-filled part of a channel with dimensions L₀ = 7 mm, h = 11.3 μm, H = 68.6 μm and w = 100 μm. Curves are fits by an exponential law L(t) = L₀e^−t/τ.

The experimental dynamics is extremely well fitted by an exponential decay:

L(t) = L(t = 0)e^−t/τ, (1)

and it turns out to be the case for all other experiments analysed in this study. This exponential decay confirms previous studies on bubble removal in microfluidic channels.^8,9 We shall henceforth focus mostly on the value of the parameter τ, which we call the imbibition timescale.

To confirm the validity of the exponential decay, we also extract the rate of decrease |L̇| by finite differences on the experimental data of L(t), and plot it as a function of L, in the inset of Fig. 3a. Although the data are more noisy, as always when estimating a time derivative of discrete data, they agree very well with a linear law of the form |L̇| = L/τ_L. Moreover, the values of τ_L measured as the best fitting slope in the inset of Fig. 3a agree with the value of τ measured as the best fitting parameter in Fig. 3a within 3%.

To assess the robustness and reproducibility of our results, we recorded five independent time evolutions L(t) in a channel with a given geometry in Fig. 3b. The corresponding five curves all collapse together, and they are all extremely well fitted by the exponential law (1). From the statistics of the five fitted timescales τ, we get τ = 47.0 ± 0.7 s: the very small value of the standard deviation compared to the mean confirms the reproducibility of our measurements.

3.2 Comparison of imbibition in closed-end and open-end channels

To clearly demonstrate the difference in imbibition dynamics with closed or open ends, we focused on one channel. We first recorded its imbibition dynamics with a dead end (open symbols in Fig. 4). Once the imbibition was completed, we let the pentanol pervaporate out of the channel, then incised the dead end of the channel to open it. We then redid the imbibition experiment (close symbols in Fig. 4). The imbibition dynamics is strikingly different: not only is the process about 500 times faster, but L(t) follows closely the square root behaviour typical of the Lucas–Washburn law in tubes of constant section,^1–5 in marked contrast with the exponential dynamics for the closed-end channel. This simple comparison convincingly demonstrates that whilst the driving force for the imbibition is the same in both cases, stemming from the capillary pressure, the resistive force is completely different: viscous dissipation for the open-end channels vs. gas transfer through PDMS for the close-end channels.

Fig. 4 Time evolution of the length of the air-filled part of a channel with dimensions w = 100 μm, h = 40 μm and H = 125 μm, either with a closed end (○) or with an open end (●). The time T equals 2500 s for the closed-end experiment and 5 s for the open-end experiment. The curves are fit by the exponential law (1) for the closed-end channel, and by a square-root law: for the open-end channel.

3.3 Dependence of the imbibition timescale on the geometrical parameters

We now study the influence of the various geometrical parameters of the channel: its total length, width, height and upper wall thickness, on the imbibition timescale τ. We first plot the imbibition timescale as a function of the total channel length, for different widths, at fixed height and thickness, in Fig. 5a. It shows that the imbibition timescale is almost independent of the total length in our studied range of parameters (here L₀ is varied between 4.5 and 18 mm): its relative variation, defined as the ratio between the standard deviation and the mean of τ for each series (see caption of Fig. 5a for the values) does not exceed 2.5%. Fig. 5a also confirms that τ increases at increasing width.

Fig. 5 (a) Plot of the imbibition timescale as a function of the channel length, for channels of dimensions h = 60 μm, H = 116.5 μm, and w = 60 (○), 80 (□), 100 (◊), 120 (△) and 140 μm (▽). Lines are average of the data; the statistics over each series of four data points gives: T = 525 ± 4, 714 ± 4, 861 ± 7, 1000 ± 15 and 1105 ± 28 s respectively for w = 60, 80, 100, 120 and 140 μm. (b) Plot of the imbibition timescale as a function of the channel width, for channels of dimensions h = 60 μm, L₀ = 18 mm, and δ = 57 (◊) and 127 μm (△). The linear size of the symbols is proportional to the channel width. The curves are fits by (11), where a = (DS^∞_airγ)⁻¹ is taken as a free fitting parameter, with value a = 0.0258 and 0.0231 μm⁻³ s respectively for δ = 57 and 127 μm. (c) Plot of the imbibition timescale as a function of the channel height, for fixed channel length L₀ = 18 mm and width w = 100 μm, and fixed upper wall thickness δ = 57.9 ± 1.8 μm. The curve is a fit by (11), where a = (DS^∞_airγ)⁻¹ is taken as a free fitting parameter, with value a = 0.0259 μm⁻³ s. For (a), (b) and (c), the chosen values of the parameters (L₀, w, h and δ) were selected to enable the independent variation of each geometric factor, ensuring that their individual influence on the imbibition timescale could be clearly assessed.

Next, we plot the imbibition timescale as a function of the channel width, for two different upper wall thicknesses, at fixed channel length and height, in Fig. 5b. It evidences the increase of the imbibition timescale at increasing width, with a slight concave trend.

Finally, we plot the imbibition timescale as a function of the channel height, at fixed channel width and length and fixed upper wall thickness, in Fig. 5c. It is much more difficult to obtain these data points, because each channel height requires making a new mold, and because aiming at a constant upper wall thickness requires to adjust the spin coating velocity of PDMS with some trial and error (whence the error bar on δ mentioned in the caption of Fig. 5c, which comes from the statistical dispersion of the three measured values of δ). For this reason, we could obtain only three data points in Fig. 5c. However, they are noteworthy because they show that the imbibition timescale increases strongly, and in a superlinear fashion, with the channel height.

3.4 Effect of the liquid

To ensure that the dynamics discussed in this work does not only apply to a specific liquid, we studied imbibition in four channels of same height and thickness and different width with two liquids: pentanol and octanol (Fig. 6). We recover the exponential decay of L(t) for octanol, and the data points are almost undistinguishable from those for pentanol. The ratio of the imbibition timescales for pentanol τ_pent and for octanol τ_oct is plotted in the inset of Fig. 6, showing a constant value of 1.03 independent of the channel width.

Fig. 6 Time evolution of the length of the air pocket, for four channels with dimensions L₀ = 18 mm, h = 60 μm, H = 180 μm, and w = 80 (blue), 100 (green), 120 (violet) and 140 μm (red). The channels were filled with pentanol (○) and octanol (▷). The lines (simple for pentanol, dotted for octanol) are fits by eqn (1). Inset: ratio of the imbibition timescales of pentanol τ_pent and octanol τ_oct, as a function of channel width.

4 Theory

The goal of this Section is to predict the imbibition timescale as a function of the geometrical parameters of the channels and of the material parameters of the PDMS.

Since pentanol is wetting the channel walls, Laplace pressure implies that the air pocket inside the channel is in overpressure relative to the liquid in contact with the meniscus, the pressure difference being equal to γ[capital script C] with γ the surface tension and [capital script C] the curvature of the meniscus. Assuming that pressure variations are negligible in the liquid, an assumption that we will shortly justify, the air trapped in the channel is therefore in overpressure relative to the ambient air around the channel. Hence, the air bubble tends to permeate outwards, which justifies the decrease of L.

4.1 Air flux through PDMS

We assume that the channels do not deform under the action of the overpressure, and we neglect the transport of pentanol through PDMS with respect to the transport of air. Hence, the instantaneous volume of the air pocket is hwL, and the conservation of air volume imposes that:

hwL̇ = −Q, (2)

where Q is the volumetric flux of air leaving the channel. In the absence of significant pressure variation throughout the experiment, we assume that the air density remains constant, whence the validity of volume (instead of mass) conservation. Moreover, the air pocket is slender, in the sense that its length is much greater than all cross-section dimensions h, w and H, except at the very end of the filling process. Hence, neglecting end effects, the permeation rate per unit length q is constant all along the air pocket, which is at constant pressure, and the outwards air flux is simply proportional to L: Q = qL. Together with (2), this justifies the experimentally observed exponential decay (1), and predicts that the imbibition timescale is:

The next step consists in predicting q. This permeation rate per unit length arises from the outwards diffusion of air within the PDMS resulting from the pressure difference between the air pocket and the atmosphere; using Fick's law, we have:

q = −∫D∇c·nds, (4)

where D = 3.4 × 10⁻⁹ m² s⁻¹ is the diffusion constant of dissolved air within PDMS²⁷ and c its concentration, and where the integral is taken along the channel/PDMS boundary (in the cross-section), of curvilinear abscissa s and unit outwards vector n. Notice that since q is a volumetric flux per unit length, it has the same dimensions as D, hence c here is dimensionless: its meaning is the dissolved volume of gas per unit volume of PDMS. To proceed, we must determine the concentration field, which obeys the diffusion equation. Since the typical order of magnitude of the transverse dimensions (e.g. the channel width or the PDMS thickness) is [small script l] ≈ 10⁻⁴ m, the typical diffusion time in the cross-section is [small script l]²/D ≈ 3 s, which is two orders of magnitude smaller than the measured values of τ in almost all experiments (see sec. 3.3), except for the smallest channel heights. Hence, we may consider that the diffusion is quasisteady, and the concentration field obeys the two-dimensional Laplace equation in the cross-section.

To get the boundary conditions, it is required to know the amount of dissolved air in the PDMS in the range of pressures of the experiments. For this, we use the solubility coefficient S̄,²⁷ defined such that c = S̄p at the PDMS/gas boundaries, where p is the partial pressure of air. This linear relationship is an example of Henry's law, and it is valid because the air concentration in PDMS is small. Hence, we take the following boundary conditions: at the PDMS/atmosphere interface, c = S̄p_atm, with p_atm = 1 atm = 10⁵ Pa, and at the PDMS/channel interface, c = S̄p_in, with p_in the pressure inside the gas pocket, neglecting the influence of the pentanol vapour in the gas pocket. Finally, there is no flux at the glass/PDMS interface.

Let us define the following dimensionless function related to the concentration field:

then f obeys the two-dimensional Laplace equation (see Fig. 1a for the definition of the coordinates):
with boundary conditions:

f = 1 at y = h and |x| ≤ w/2, and at |x| = w/2 and 0 ≤ y ≤ h.

The second of these three conditions expresses that glass is impermeable to air. This problem has been solved analytically in Appendix A of ref. 21. In particular, as long as δ/w ≤ 1.5 or h/w ≤ 0.1, a condition met in our experiments, the following approximation was found to hold within 1%:

The two terms in the right-hand side can be interpreted as follows. The term w/δ represents the “direct” transfer through the PDMS wall of width w and thickness δ separating the channel top from the atmosphere. The second term represents the transfer from the sides of the channel to the atmosphere through the thicker body of PDMS. In our experimental conditions, this side-transfer term is far from negligible.

Coming back to the definition (5) of f and plugging (7) into (4), we thus get the flux:

4.2 Capillary pressure

To predict the pressure difference p_in − p_atm remaining in (8), we first split it in two parts: p_in − p_atm = (p_in − p_{[small script l]}) + (p_{[small script l]} − p_atm), where p_{[small script l]} in the pressure in the liquid in contact with the meniscus. Hence, p_in − p_{[small script l]} is the Laplace pressure across the meniscus. A difficulty in our experiments is that the meniscus is in contact with two different materials: PDMS at the side walls and at the top wall, and glass at the bottom wall. Since glass is more hydrophilic than PDMS, the contact angle of pentanol on glass θ′ is lower than the contact angle θ = 30° of pentanol on PDMS. To account for this presence of two different materials, which is an ubiquitous situation in microfluidics owing to the frequent gluing of PDMS channels on glass slides, the following expression of the capillary pressure has been claimed to apply without justification²⁸ and quoted in subsequent review articles:²⁹p_in − p_{[small script l]} = γ[(cos θ + cos θ′)/h + 2 cos θ/w]. However, this unjustified formula is erroneous, and we derive in the Appendix the correct formula for arbitrary, extending a method introduced by Mason and Morrow.³⁰ Taking for simplicity θ′ = 0, hence assuming that pentanol is perfectly wetting on glass, this formula reads:In our experiments, γ = 25.7 mN m⁻¹ for pentanol and h is at most 80 μm, hence the order of magnitude of p_in − p_{[small script l]} is γ/h ≈ 3 × 10² Pa.

We now argue that |p_{[small script l]} – p_atm| ≪ p_in − p_{[small script l]}. Three terms may contribute to p_{[small script l]} – p_atm: the Laplace pressure at the surface of the drop, a gravitational term over the height of the drop, and a viscous pressure drop in the liquid set into motion by the bubble shrinkage. In our experiments, the drop is flat and Laplace pressure is certainly orders of magnitude lower than at the meniscus between the bubble and the liquid (Fig. 1 is not at scale). The height of the drop is less than 1 mm, hence hydrostatic pressure remains lower than 10 Pa. To estimate the viscous pressure drop, we assume that a length L₀ − L of the liquid in the channel is entrained at the velocity of the meniscus |L̇| (hence we neglect pervaporation effects), and we use Poiseuille law in a rectangular channel: |L̇| = [scr S, script letter S]∇p/η, with ∇p the pressure gradient, η = 4 × 10⁻³ Pa s the liquid viscosity and an effective section given by (see e.g.³¹):

which remains close to h²/12 in our range of experimental parameters. Hence, a good estimation of the viscous pressure drop is 12η(L₀ − L)|L̇|/h². With the exponential law L(t) = L₀e^−t/τ, the quantity (L₀ − L)|L̇| has as upper bound L²₀/(4τ). Hence, the viscous pressure drop is lower than 3ηL²₀/(h²τ). This quantity is equal to 2 Pa with L₀ = 18 mm, h = 60 μm and τ = 5 × 10² s (an order of magnitude given by Fig. 5b). This remains more than two orders of magnitude lower than Laplace pressure (9). Altogether, the pressure variations within the liquid remain negligible with Laplace pressure (9), which we can therefore identify to the term p_in − p_atm in (8). Therefore, we get the final prediction

Inserting this expression in (3), we predict the imbibition time:

which constitutes the main theoretical result of our study. In particular, although we shall compare this prediction to data on pentanol and octanol in hydrophobic PDMS, eqn (11) can be applied to any liquid and any surface treatment, provided the contact angle θ is measured. Some studies also reported that surface treatment can also affect the diffusion coefficient, which decreases for freshly oxidised PDMS:³² this possible effect is also accounted for in our model, since eqn (11) shows that the imbibition time is inversely proportional to the diffusion coefficient.

5 Comparison between experiments and theory

We now compare this prediction to the experimental measurements, first focusing on the effect of geometry. First, we remark that the imbibition timescale is predicted not to depend on the channel length, which agrees with Fig. 5a. Next, we fit the data from Fig. 5b and c by the formula (11) with a = (DS̄γ)⁻¹ taken as a fitting parameter. These fits show that the model captures qualitatively all the trends, namely the increase of the imbibition timescale at increasing width, increasing upper wall thickness and increasing channel height. In particular, it reproduces excellently the dependence on the channel height (Fig. 5c), whilst it predicts slightly stronger variations with the channel width than the measurements (Fig. 5b). We also notice that the obtained values of the fitting parameter a, taken independently for each series of data, are consistent within 7%; averaging them yields a = 0.0240 ± 0.0016 μm⁻³ s. This suggests that our model includes the essential ingredients to understand the role of the geometrical parameters on the imbibition dynamics.

We can also test the predictions for different liquids, here pentanol and octanol (Fig. 6). According to (11), at given channel geometry, the ratio of the imbibition times of the two liquids should be constant and equal to the inverse ratio of their surface tensions: τ_pent/τ_oct = γ_oct/γ_pent at given contact angle θ. The imbibition timescale ratio was found to be almost constant in experiments (see inset in Fig. 6). Using the values γ_pent = 25.4 mN m⁻¹, γ_oct = 27.1 mN m⁻¹, given for T = 25 °C in,³³ and comparing with the values obtained from curve fitting, we find that the results are predicted well by the theory: we get τ_pent/τ_oct = 1.03 whilst γ_oct/γ_pent = 1.07. The remaining small discrepancy likely comes from a difference of contact angle between pentanol and octanol on PDMS.

We now discuss the value of the fitting parameter a. To have a numerical value for the solubility, we use the data of Merkel et al.,²⁷ who measured the solubility S of various gases in PDMS. Merkel et al. showed that the solubility is well represented by an affine law S = S^∞ (1 + np), and provide the following values (see Tab. II in²⁷) for N₂: and n_N2 = 3 × 10⁻³ atm⁻¹, and for O₂: and n_O2 = 5 × 10⁻³ atm⁻¹. Since our pressures remain of order 1 atm = 10⁵ Pa, these numerical data show that we may take S = S^∞ as a good approximation. In the value of the solubility, 1 cm³_STP is the number of moles that would occupy one cubic centimeter at standard temperature and pressure, as calculated via the ideal gas law. In our model, since c is dimensionless, we use the following unit conversion: S̄ = 1 cm³_STP cm⁻³ atm⁻¹ corresponds to S̄ = 1 atm⁻¹ = 10⁻⁵ Pa⁻¹. Hence, if we consider a simplified atmosphere composed of 79% of O₂ and of 21% of N₂, we take the linear combination of the contributions of these two gases (which is justified by the additivity of the partial pressures), whence the following value of the air solubility: . Hence, we shall use S̄ = 1.1 × 10⁻⁶ Pa⁻¹. With the values already mentioned: D = 3.4 × 10⁻⁹ m² s⁻¹ and γ = 25 mN m⁻¹, this gives (DS̄γ)⁻¹ = 0.0105 μm⁻³ s. This value agrees in order of magnitude with our fitting value a = 0.024 μm⁻³ s, but it is two times lower. We shall discuss later some reasons for this discrepancy, in relation with some simplifying assumptions of the model. However, a probable explanation for this factor-of-two difference is the specific PDMS recipe used here, which differs from the preparation reported by Merkel et al.²⁷ In particular, Lamberti et al.³⁴ demonstrated that the gas permeability of PDMS can vary by a factor of four when the prepolymer-to-curing-agent ratio is changed. Physically, this variation arises from alterations to the free-volume “holes” in the polymer matrix–a notion supported by earlier observations of Carrillo et al.³⁵ and Stafie et al.,³⁶ which highlight how increasing the crosslinker content reduces the size and number of these voids, thereby curbing gas transport. Consequently, even moderate shifts in PDMS preparation could readily explain the gap between our parameter a and the value inferred from Merkel's data. It is also worth noting that Merkel et al. purchased their PDMS membranes and do not provide quantitative details about the curing agent, hence it is probably illusory anyway to expect a 1 : 1 comparison between their reported solubilities and the ones inferred from our study.

6 Channels of varying width

We now extend our study to the case of single channels of varying width. We start by the extension of the theory to such a case which, as we shall see, is relatively straightforward. An analogous treatment was realised in the context of drying channels by pervaporation.³⁷ We use a streamwise coordinate x along the channel, choose its origin x = 0 at the close end of the channel, and orient x towards its open end; with this convention, the meniscus is located at x = L, and its velocity L̇ is negative. We may then extend the conservation of air (2) simply by evaluating the width at the meniscus location:

hw(x = L)L̇ = −Q, (12)

with Q the volumetric flux of air leaving the channel. The latter is easy to evaluate if the cross-section dimensions are small with respect to the longitudinal scale L_w along which w varies significantly. The cross-section dimensions (w, δ and h) determine the order of magnitude of the length over which the concentration field varies across the channel, whilst L_w determines the order of magnitude of the length over which the concentration field varies along the channel. Hence, if L_w is much larger than w, δ and h, we can still use the expression (8) of the permeation rate per unit length q derived in sec. 4.1, provided that in the term between braces in (8), we use the local value of the channel width. The volumetric flux Q is then obtained by integrating q over the length of the air pocket. On the other hand, the Laplace pressure term p_in − p_atm, given by (9), depends on the width at the meniscus location. Hence,

Once the profile w(x) is known, it is a simple matter to compute the remaining integral, and to insert the expression of Q into the equation of conservation of air (12) to obtain a (nonlinear) differential equation for the dynamics of the meniscus.

To test this idea experimentally, we designed two channels with linearly variable width (Fig. 2b). They are both of length L_c = 5.5 mm, and their width varies linearly between 80 and 225 μm, one with increasing width, and the other with decreasing width; hence, they contain the same volume. We plot in Fig. 7a the time evolution of the length of the gas pocket for this two channels. Interestingly, whilst this length decreases initially faster in the channel with increasing width, there is a crossover at t = 4 × 10³ s in Fig. 7a, and eventually the gas pocket disappears first in the channel with decreasing width. Moreover, the shape of the curves is no longer exponential as it was for channels of uniform width (Fig. 3). Indeed, plotting |L̇| as a function of L shows that whilst we get straight lines for channels of uniform width, we get a convex curve for the channel of increasing width, and a concave curve for the channel of decreasing width (Fig. 7b).

Fig. 7 (a) Time evolution of the length of the air-filled part of two channels, one of increasing width (◁) and one of decreasing width (▷). The recording for this experiment has been stopped at 7000 s. (b) Corresponding plot of the velocity of the air/liquid interface as a function of the length of the air-filled part; curves are fits by eqn (14). The inset zooms in the region of low values of L and |L̇|.

To capture the experimental results, we use the following expression for the width varying linearly between w₁ at the close end (x = 0) and w₀ at the open end (x = L_c): w(x) = w₁ + (w₀ – w₁)x/L_c. Inserting this expression in (13) and using (12), we get the following differential equation:

with dimensionless parameters L̄ = L/L_c, w̄ = w₀/w₁ and h̄ = h/w₁. This is a separable equation which can be solved analytically to yield an implicit solution of the form t(L̄), but this is a tedious task, and it is much more direct to fit the data from Fig. 7b with the right-hand side of (14), with a as a free fitting parameter. The fitting curve is in excellent agreement with the data. This validates the extension of our analysis of imbibition for channels of varying width.

7 Discussion

7.1 Approximations in the model

Owing to the wetting properties of pentanol, four liquid gutters are expected to form in the corners of the channel surrounding the bubble. It is known that when they are present, such gutters influence strongly confined drops and bubbles: they play a significant role in the motion of drops and bubbles in rectangular capillaries,³⁸ as well as in droplet break-up processes involving T-junctions^39–41 or gradients of confinement.⁴² They are also essential in our derivation of the capillary pressure (see Appendix). In our experiments, they are difficult to detect owing to their small cross-section, but they show by slightly thicker dark boundaries along the sides of the channel ahead of the meniscus. However, their extension is limited to about one millimetre only. Such a limitation is not observed in glass capillaries. It is likely due to the quick pervaporation of the pentanol within the gutters through the PDMS. It is because of their small length that we neglected the presence of gutters to compute the air flux through PDMS in our model.

Moreover, the channel may deform and bulge out because the gas pocket is capillary overpressure p_in − p_{[small script l]} with respect to the atmosphere. To get an estimate of this deformation, we proceed as in,⁴³ and assume that the PDMS layer between the channel top and the atmosphere behaves as a thin elastic plate. Even though this layer is not slender in our experimental conditions, because δ is not much smaller than w, this should give the correct order of magnitude of its deflection. The latter is maximal at the centre of the channel, where it equals:⁴⁴ζ_max = (p_in − p_{[small script l]})w⁴/ (384B), with a bending modulus B = Eδ³/[12(1 – ν²)] with E = 2 MPa the Young modulus of PDMS, and v = 0.5 its Poisson ratio. For the typical height h = 60 μm of our channels and their largest width (for which the deflection is maximal) w = 140 μm, and the expression (9) of the capillary pressure, we get the numerical evaluation: ζ_max = 24 nm, which is negligible in comparison with the channel height. Even for the channel of smallest height h = 11 μm (the left point in Fig. 5c), for which w = 100 μm, we get ζ_ax = 26 nm, which remains negligible. Therefore, the channel deformation can be safely neglected in our study.

7.2 Predicting bubble removal timescale for microfluidics engineering

Having established a formula that describes spontaneous imbibition in our specific configuration (a wetting liquid invading a dead-end channel under capillary action), we now discuss how to extend this approach to bubble removal in broader microfluidic contexts. In many practical devices, especially those featuring complex networks or bypass geometries, trapped air pockets cannot be evacuated by a main flow alone. These bubbles typically inhabit dead-end branches or become pinned in regions unfavourable for advection-driven clearance. When the contact angle is not sufficiently low, the capillary pressure may be too weak to eliminate such gas pockets on its own. A commonly employed solution in microfluidics engineering is to impose an inlet pressure p_{[small script l]} exceeding ambient pressure. In such a setup, operators normally stop or block other flows: the outlet, if present, is either sealed or set to the same pressure as the inlet. The gas pocket therefore remains stationary, whilst its internal pressure p_in becomes:

p_in = p_{[small script l]} + p_cap, (15)

where p_cap is the net capillary contribution (which can be positive for favourable wetting or negative otherwise), and is given by (9) in its range of validity. Even a modest overpressure can significantly speed up the gas diffusion through the surrounding PDMS, thereby removing the bubble more quickly. This generalises our earlier analysis of spontaneous imbibition: rather than relying solely on the meniscus curvature to drive diffusion, the total bubble pressure now reflects both p_cap and the imposed inlet pressure p_{[small script l]}. Retaining the same treatment of the net diffusion flux through the channel walls, we obtain a generic expression for the characteristic timescale of bubble elimination:where h and w denote the channel height and width, D is the diffusion coefficient of air in PDMS, and S̄ the air solubility. The term p_env is the external (often atmospheric) pressure, which can be made down to very low values by operators using vacuum to boost bubble removal. The term Γ(δ, w, h) encapsulates geometric factors dictating the overall permeation flux, as derived in a previous article related to pervaporation:²¹with δ the thickness of the upper wall and H the total thickness of the chip, such that H = h + δ.

From an experimental standpoint, the parameter not always tabulated is the product DS̄ which can differ significantly based on fabrication factors – such as curing-agent ratio or post-baking steps – especially for elastomers like PDMS (see our discussion at the end of sec. 5). Once a single calibration experiment is performed for a given material (and preparation recipe), operators can extract this product, then apply eqn (16) to predict bubble-removal timescales under new geometric or pressure conditions.

Although eqn (16) was derived assuming uniform-width channels, an analogous treatment can be extended to more complex geometries like those in sec. 6, or even intricate network topologies. Moreover, whilst our derivation focused on PDMS due to its wide usage and high gas permeability, the same conceptual approach can be extended to other chip materials so long as their air-diffusion and solubility properties are known or can be calibrated. This versatility makes the overall framework broadly applicable to a range of bubble-management challenges in microfluidic devices.

Conclusions

Our work provides a comprehensive experimental and theoretical framework for predicting the time required to remove air pockets from dead-end microchannels via gas diffusion in PDMS. We showed that spontaneous imbibition of a wetting liquid, combined with gas permeation through the channel walls, follows a simple exponential law whose timescale is set by geometry, capillarity, and PDMS material properties. Crucially, we quantified how variations in channel width, height, and material thickness can shift the dynamics of bubble shrinkage by orders of magnitude. By developing a physically grounded model and validating it against experimental measurements, we offer a robust tool for microfluidic operators confronted with trapped air pockets. This approach not only simplifies device operation and maintenance but also lays the groundwork for future advances in pumpless and passive microfluidics. Additionally, in deriving the interface curvature for channels composed of different wall materials, a configuration commonly seen in PDMS-on-glass microfluidics, we correct some existing erroneous predictions and provide an expression with broad relevance, applicable to diverse droplet-based or two-phase flow configurations well beyond this study.

Future explorations would profit from extending our approach to complex network architectures, such as branched and interconnected microfluidic layouts in which individual branches can compete or couple, making bubble removal more intricate. Another promising direction involves incorporating the deformability of soft channels, where elastocapillary effects become significant and channel walls can distend in response to local capillary pressures, potentially leading to unexpected dynamics. Furthermore, introducing strong geometric constrictions, as opposed to the incremental width variations explored here (sec. 6), may produce sharply nonlinear imbibition scenarios, much like what could be observed in the reverse configuration of air penetration in liquid-filled channels with constrictions.⁴³ Indeed, such constrictions parallel the functionality of bordered pits in plant xylem,^9,45,46 so investigating their influence in PDMS microchannels could provide fresh insights relevant to both biomimetic design and microfluidic applications.

A. Derivation of the capillary pressure

In this Appendix, we derive the capillary pressure p_in − p_{[small script l]} in the channel bounded by PDMS at the top, left and right edges (contact angle θ) and by glass at the bottom edge (contact angle θ′). Crucially, we assume that the contact angles are low enough that wedge-like gutters form ahead of the meniscus,⁴⁷ as shown in the sketch of the cross section of the channel in Fig. 8. The success of our method relies on this assumption: as we shall see, it enables to bypass the direct computation of the meniscus curvature which can be done only numerically. For simplicity, we shall reason on a half-width of the cross section, without loss of generality.

Fig. 8 Representation of the half-width cross section of the channel in the case w/h = 1.5, θ′ = 25° and θ′ = 10°. The gutters are the circular triangles, filled in light blue, at the top right and bottom right of the channel. The perimeter of the air-filled regions comprises three parts: a part separating air and PDMS, of length [scr P, script letter P]_P (green lines FA and BA′); a part separating air and glass, of length [scr P, script letter P]_g (red line F′B′); and a part separating air and liquid at the gutter boundaries, of length [scr P, script letter P]_L (blue curves AB and A′B′).

On Fig. 8, we denote [scr A, script letter A] the area filled with air, which spans the cross section excluding the gutters. The outer perimeter [scr P, script letter P] of this air-filled region is split in three parts: [scr P, script letter P] = [scr P, script letter P]_L + [scr P, script letter P]_P + [scr P, script letter P]_g, where [scr P, script letter P]_L is the perimeter between air and liquid, [scr P, script letter P]_P the perimeter between the air and PDMS, and [scr P, script letter P]_g the perimeter between the air and glass. Considering a virtual displacement dx of the meniscus in the longitudinal direction of the channel yields the following virtual work balance: (p_in − p_{[small script l]})[scr A, script letter A]dx = γ([scr P, script letter P]_L + [scr P, script letter P]_P cos θ + [scr P, script letter P]_g cos θ′)dx. This equation is a straightforward extension of the approach of Mason and Morrow,³⁰ who considered a single contact angle, to the case of two contact lines on two solid materials with different contact angles. The capillary pressure is related to the mean radius of curvature r of the meniscus by p_in − p_{[small script l]} = γ/r; and r is also the radius of curvature of the gutters far ahead the meniscus, as their cross section becomes independent of x. Hence, r obeys the relation:

We now consider the geometry of the gutters, to transform (18) into a equation of the second degree relating r to h, w, θ and θ′. First, [scr P, script letter P]_L is the curvilinear length of the two arcs of circle AB and A′B′:

[scr P, script letter P]_L = r(π − 3θ −θ′). (19)

Then, from elementary geometry, O′D′ = r cos θ, D′A′ = r sin θ, O′E′ = r cos θ′ and E′B′ = r sin θ′. Hence, B′C′ = E′C′ − E′B′ = O′D′ − E′B′ = r(cos θ − sin θ′). Since B′C′ must be positive for a gutter to form, a necessary condition for the existence of a gutter at the corner between PDMS and glass is therefore cos θ > sin θ′. We thus obtain:

We also have A′C′ = O′E′ − A′D′ = r(cos θ − sin θ′). Similarly, we obtain OD = OE = r cos θ, and DA = EB = r sin θ, hence AC = BC = r(cos θ − sin θ), which shows that a necessary condition for the existence of a gutter at the corner between two PDMS sides is θ < π/4, as already known.⁴⁷ Hence, . Moreover, A′B′ = h – BC − A′C′ = h – r(cos θ + cos θ′ − 2 sin θ). We can thus express the perimeter between air and PDMS:

Finally, . Now, from elementary geometry, . Similarly, . Hence,

Inserting (19)–(22) in (18), and simplifying after some algebra, yields the equation of second degree for the gutter curvature κ = 1/r:

which can be solved to obtain κ. We distinguish in the following three different cases, depending on the value of θ′. First expressing the curvature for two particular cases (θ = θ′ and θ′ = θ), and then developing the general expression of κ for arbitrary values of θ and θ′.

A.1 case 1: θ = θ′

In the simple case, corresponding to the four walls made with the same materials, one has θ = θ′. eqn (23) is readily solved, and a result obtained by Wong et al.⁴⁸ is recovered:

A.2 case 2: θ′ = θ

Most relevant for our study is the case θ′ = θ, which corresponds to the case of perfect wetting on glass. The equation of second degree on κ then becomes:which has a discriminant , which is obviously positive in the relevant range 0 ≤ θ < π/4 for the contact angle. Hence, there are two real solutions κ_±, and a quick study of the signs of the coefficients of (24) show that the two solutions are positive, hence more work is required to determine which one gives the true curvature. However, for a positive solution to be physically admissible, two adjacent gutters must not merge, which requires A′B > 0 on Fig. 8. Hence, from the previous geometrical calculations, this sets a lower bound for the curvature: κ ≥ κ_c = (1 + cos θ – 2 sin θ)/h. A quick calculation yields , which is obviously negative. But f(κ) defined in (24) is a second-order polynomial whose coefficient of order 2 is positive, hence κ_c lies between its two roots: κ₋ < κ_c < κ₊, which proves that the only admissible curvature is:
whence the expression of the capillary pressure p_in − p_{[small script l]} = γκ₊ used in the main text in section 4.

A.3 case 3: arbitrary values for θ and θ′

In the general case, for any value of θ and θ′, one can follow the same procedure, which yields the discriminant from eqn (23). Keeping only the largest root as the correct physical curvature, one gets:

This expression for the curvature can be taken directly to deduce the value of ΔP_cap = p_in − p_{[small script l]} = γκ₊ in the final discussion (sec. 7) when discussing the generic expression to be used by microfluidic operators.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article, including data files and scripts to generate the figures, are available at Zenodo at https://doi.org/10.5281/zenodo.16356173.

Acknowledgements

This work has been supported by the French government, through the UCA Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01 and through the grant ANR-19-CE30-0010 (PHYSAP). The authors thank technical help from D. Centanni and M. Van Melle-Gateau.

Notes and references

1. J. M. Bell and F. K. Cameron, J. Phys. Chem., 1906, 10, 658–674.
2. R. Lucas, Kolloid-Z., 1918, 23, 15–22.
3. E. W. Washburn, Phys. Rev., 1921, 17, 273.
4. P. G. de Gennes, F. Brochard-Wyart and D. Quéré, Gouttes, bulles, perles et ondes, Belin, 2002.
5. J. B. Gorce, I. J. Hewitt and D. Vella, Langmuir, 2016, 32, 1560–1567.
6. O. Vincent, J. Zhang, E. Choi, S. Zhu and A. D. Stroock, Langmuir, 2019, 35, 2934–2947.
7. M. A. Zwieniecki and N. M. Holbrook, Trends Plant Sci., 2009, 14, 530–534.
8. J. H. Kang, Y. C. Kim and J. K. Park, Lab Chip, 2008, 8, 176–178.
9. L. H. Guo, J. Shan, R. H. Ran, S. Q. Yin, C. Liu and J. M. Li, Langmuir, 2022, 38, 12373–12381.
10. C. R. Brodersen, A. J. McElrone, B. Choat, M. A. Matthews and K. A. Shackel, Plant Physiol., 2010, 154, 1088–1095.
11. S. Delzon and H. Cochard, New Phytol., 2014, 203, 355–358.
12. K. H. Jensen, K. Berg-Sørensen, H. Bruus, N. M. Holbrook, J. Liesche, A. Schulz, M. A. Zwieniecki and T. Bohr, Rev. Mod. Phys., 2016, 88, 035007.
13. L. Xu, H. Lee, D. Jetta and K. W. Oh, Lab Chip, 2015, 15, 3962–3979.
14. M. Boyd-Moss, S. Baratchi, M. D. Venere and K. Khoshmanesh, Lab Chip, 2016, 16, 3177–3192.
15. A. Wang, D. Koh, P. Schneider, E. Breloff and K. W. Oh, Micromachines, 2019, 10, 543.
16. J. Park, D. H. Han and J. K. Park, Lab Chip, 2020, 20, 1191–1203.
17. L. F. Xu, A. Y. Wang, X. P. Li and K. W. Oh, Biomicrofluidics, 2020, 14, 031503.
18. V. Narayanamurthy, Z. E. Jeroish, K. S. Bhuvaneshwari, P. Bayat, R. Premkumar, F. Samsuri and M. M. Yusoff, RSC Adv., 2020, 10, 11652–11680.
19. P. Bacchin, J. Leng and J. B. Salmon, Chem. Rev., 2022, 122, 6938–6985.
20. X. Noblin, L. Mahadevan, I. A. Coomaraswamy, D. A. Weitz, N. M. Holbrook and M. A. Zwieniecki, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 9140–9144.
21. B. Dollet, J. F. Louf, M. Alonzo, K. H. Jensen and P. Marmottant, J. R. Soc., Interface, 2019, 16, 20180690.
22. Y. Liu, J. Li and A. J. Smits, J. Fluid Mech., 2019, 876, 1129–1145.
23. P. Kolliopoulos, K. S. Jochem, R. K. Lade Jr., L. F. Francis and S. Kumar, Langmuir, 2019, 35, 8131–8143.
24. A. Gökaltun, Y.-B. Kang and M. L. Yarmush, et al. , Sci. Rep., 2019, 9, 7377.
25. J. N. Lee, C. Park and G. M. Whitesides, Anal. Chem., 2003, 75, 6544–6554.
26. S. Basu, D. C. K. Rao, A. Chattopadhyay and J. Chakraborty, Phys. Rev. E, 2021, 103, 013101.
27. T. C. Merkel, V. I. Bondar, K. Nagai, B. D. Freeman and I. Pinnau, J. Polym. Sci., Part B: Polym. Phys., 2000, 38, 415–434.
28. E. Delamarche, A. Bernard, H. Schmid, A. Bietsch, B. Michel and H. Biebuyck, J. Am. Chem. Soc., 1998, 120, 500–508.
29. A. Olanrewaju, M. Beaugrand, M. Yafia and D. Juncker, Lab Chip, 2018, 18, 2323–2347.
30. G. Mason and N. R. Morrow, J. Chem. Soc., Faraday Trans., 1984, 80, 2375–2393.
31. H. Bruus, Theoretical Microfluidics, Oxford University Press, 2008.
32. D. A. Markov, E. M. Lillie, S. P. Garbett and L. J. McCawley, Biomed. Microdevices, 2014, 16, 91–96.
33. A. H. Johnstone, J. Chem. Technol. Biotechnol., 2007, 50, 1275–1276.
34. A. Lamberti, S. L. Marasso and M. Cocuzza, RSC Adv., 2014, 4, 61415–61419.
35. F. Carrillo, S. Gupta, M. Balooch, S. J. Marshall, G. W. Marshall, L. Pruitt and C. M. Puttlitz, J. Mater. Res., 2005, 20, 2820–2830.
36. N. Stafie, D. Stamatialis and M. Wessling, Sep. Purif. Technol., 2005, 45, 220–231.
37. K. N. Chagua Encarnación, P. Marmottant and B. Dollet, Microfluid. Nanofluid., 2021, 25, 71.
38. H. Wong, C. J. Radke and S. Morris, J. Fluid Mech., 1995, 292, 95–110.
39. V. van Steijn, C. R. Kleijn and M. T. Kreutzer, Phys. Rev. Lett., 2009, 103, 214501.
40. P. M. Korczyk, V. van Steijn, S. Blonski, D. Zaremba, D. A. Beattie and P. Garstecki, Nat. Commun., 2019, 10, 2528.
41. J. Zhou, Y.-M. Ducimetière, D. Migliozzi, L. Keiser, A. Bertsch, F. Gallaire and P. Renaud, Phys. Rev. Fluids, 2023, 8, 054201.
42. R. Dangla, S. C. Kayi and C. N. Baroud, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 853–858.
43. L. Keiser, P. Marmottant and B. Dollet, J. Fluid Mech., 2022, 948, A52.
44. L. D. Landau and E. M. Lifshitz, Theory of Elasticity, Elsevier, 3rd edn, 1986.
45. O. Vincent, D. Sessoms, E. Huber, J. Guioth and A. Stroock, Phys. Rev. Lett., 2014, 113, 134501.
46. L. Keiser, B. Dollet and P. Marmottant, J. R. Soc. Interface, 2024, 21, 20240103.
47. P. Concus and R. Finn, Proc. Natl. Acad. Sci. U. S. A., 1969, 63, 292–299.
48. H. Wong, S. Morris and C. J. Radke, J. Colloid Interface Sci., 1992, 148, 317–336.

---