Wettability-gradient-driven capillary filling dynamics in architected tapered microchannels

Abstract

Capillary-driven transport is central to soft and biological matter, from plant-xylem water ascent to autonomous flows in microfluidic networks. Here, we systematically investigate autonomous capillary filling dynamics in microchannels combining geometric tapering and spatially variable wettability. Using high-resolution computational fluid dynamics (Navier–Stokes equations and the level-set method), we quantify the impact of stepwise, linear, and quadratic contact-angle profiles on the Laplace pressure, interface morphology, and flow velocity. For uniform channels and contact angles, the simulations reproduce the classical Lucas–Washburn regime, characterized by a viscous slowdown. In contrast, geometric tapering amplifies the capillary pressure gradient, sustaining or accelerating interface advancement. Tailored wettability gradients enable further control: decreasing the contact angle maintains flow, while increasing the angle toward 90° robustly halts motion, enabling on-demand interface arrest. These results reveal how geometric and interfacial patterning can be coupled for precision fluid manipulation, offering broadly applicable design principles for advanced passive microfluidic systems and programmable soft-matter transport.

---

1 Introduction

Autonomous fluid motion is observed across a wide range of biological and physical systems, enabling fluid transport without the use of external energy sources. From the directed movement of water droplets on superhydrophobic self-cleaning surfaces¹ to the self-organization of bacterial colonies² and the spontaneous transport of nutrients in plants,³ autonomous fluid motion plays a central role in many natural processes. In living organisms, motile bacteria harness chemical gradients to propel themselves through fluids,⁴ while certain insects, such as water striders, exploit surface tension to glide effortlessly across water surfaces.⁵ Similarly, phalarope birds use geometric and wetting effects to move prey in droplet form within their mouth.^6,7 Fluid transport in capillary networks in plants also occurs spontaneously through intrinsic forces like surface tension and pressure gradients.⁸ Drawing inspiration from these natural mechanisms, autonomous fluid motion in microfluidic systems has emerged as a promising strategy for achieving precise flow control without the need for external energy input. In microfluidics, this principle forms the basis of passive pumping which is an essential mechanism in devices that operate without external power sources.

Furthermore, passive microfluidic pumps utilize natural forces such as capillary action,^9–12 gravity-driven flow,^13–15 vacuum suction,^16–18 and osmosis^19–21 to drive fluid movement. Inspired by energy-efficient natural systems such as plant vasculature, these methods eliminate the need for external actuators or complex control systems, enabling autonomous fluid transport. Passive pumping technologies are inherently portable, low-cost, and easily scalable, making them particularly suited for point-of-care diagnostics and applications in resource-limited settings, where simplicity and reliability are paramount.²² Consequently, passive microfluidic pumps have become integral to a wide array of applications, including biomedical diagnostics,²³ drug delivery,¹¹ chemical sensing,⁹ and integrated lab-on-a-chip platforms.²²

Among the various autonomous fluid-transport mechanisms, capillary-driven motion stands out for its simplicity, reliability, and ease of integration into microfluidic systems.²⁴ This mechanism relies on surface tension as the primary driving force, where adhesive forces between the liquid and the channel surface dominate over cohesive forces within the liquid, enabling spontaneous fluid movement. The classical Lucas–Washburn law provides the foundational description for capillary filling in uniform channels, predicting that the interface position grows as and that the velocity decays as due to increasing viscous resistance with filling length.^25,26 This regime, validated in numerous experimental and theoretical studies, sets the benchmark for understanding capillary transport in confined geometries.

Before delving deeper into the specifics of autonomous fluid motion driven by surface tension, it is essential to appreciate the rich body of foundational and contemporary research that has shaped this area of study. Pioneering studies on surface tension-driven motion include the work of Brochard,²⁷ which elucidated the role of surface-tension gradients in liquid spreading, and Chaudhury and Whitesides,²⁸ who demonstrated spontaneous droplet motion on surfaces with wettability gradients. Sandre et al.²⁹ further advanced this understanding by analyzing contact line dynamics in wetting-driven flows. Recent work has further highlighted how channel geometry, such as tapering, convergent/divergent sections, and non-circular cross-sections, can generate Laplace pressure gradients that enhance or modulate capillary action.^30–34 These studies collectively emphasize that both channel shape and surface properties are critical handles for optimizing capillary-driven filling, particularly in microfluidic applications.

Autonomous fluid motion driven by surface tension typically occurs in two forms: wicking through porous substrates and filling within solid microchannels. In porous materials, wicking persists as long as the substrate remains hydrophilic to the liquid,^21,35–40 whereas in solid-based systems, fluid motion is confined to narrow channels or capillaries, driven by surface tension gradients.^41–43 This latter process, referred to as capillary filling, forms the basis of the functioning of capillary-driven passive micropumps discussed earlier. Extensive experimental and numerical research has explored capillary filling dynamics,^{12,24,44–49} with studies examining diverse channel geometries, such as square, triangular, convergent/divergent, and circular capillaries, to enhance flow performance.³⁰ Notably, while forces like gravity may contribute, research demonstrates that channel geometry alone, through curvature-induced Laplace pressure differences, is often sufficient to drive fluid motion via surface tension.^31–33 These investigations highlight the critical role of channel shape and surface properties in optimizing capillary-driven filling, particularly in microfluidic applications.

Capillary filling motions encounter notable challenges that impact efficiency and scalability despite their numerous advantages. One primary challenge lies in achieving consistent flow velocities and precise control over fluid dynamics, primarily due to the inherent variability of surface-tension forces.^45,50 This is coupled to the emergence of multiple regimes during capillary filling, each characterized by distinct interfacial and flow behaviors.⁵¹ Consequently, both these effects lead to inaccuracies, particularly in applications requiring precise fluid handling and control. Researchers have investigated various strategies to regulate capillary-driven fluid filling to address these limitations. For instance, various studies have investigated the impact of different factors on capillary-driven filling, such as channel dimensions and configurations,^49,52,53 patterned wall structures,⁵⁴ contact-angle gradients,^34,50 contact-line dynamics,^45,55 and the application of pressure head at the inlet.⁵⁶ Significantly, these studies highlight how these parameters impact the surface-tension force driving the fluid motion in the microchannels. Additionally, some works have also attempted to appropriately regulate the capillary-driven filling by refilling the inlet droplet⁵⁷ or intermittently stopping and restarting the inlet flow.⁵³

Despite ongoing efforts to achieve consistent flow velocities in capillary-driven microfluidic systems, research aimed at regulating capillary filling dynamics through a combination of geometric modifications and axially varying wall wettability, particularly using high-fidelity numerical models, remains relatively limited. The key limitations in capillary-driven microfluidics persist, including an emphasis on uniform wettability and simple geometries while neglecting the combined effects of tapering and spatially varying wettability for enhanced flow control. Such regulation involves not only maintaining steady flow velocities, but also enabling precise control over the flow, including the ability to increase, decrease, or halt it at specific locations within the microchannel. This capability facilitates on-demand capillary filling. Furthermore, the physics underlying the transition from classical Lucas–Washburn behavior to regimes governed by geometry- or wettability-induced pressure gradients remains incompletely understood.^33,49,55,58

The present study introduces a unique approach that combines a tapered capillary design with axially varied wall contact angles to enhance and precisely control the filling of a fluid in a microchannel to address this gap. The proposed design features a microfluidic channel with a gradually narrowing width lined by solid walls, as illustrated in Fig. 1. Initially, the microchannel is filled with one fluid displaced by another less dense fluid. This geometry induces an additional pressure drop that drives the flow. At the same time, the axial variation in wall contact angles modifies the surface tension force, providing a robust mechanism to regulate capillary filling dynamics. One can establish a framework for achieving precise fluid control during the capillary filling of fluids inside capillaries by systematically tuning these parameters.

Fig. 1 Schematic representation of the tapered microchannel, depicting the capillary filling of fluid 1 into the microchannel filled with fluid 2. The wetted walls are represented using respective labels. The radii of the inlet and outlet, and the length of the microchannel are represented using the individual symbols R₁, R₃, and L, respectively. The position of the interface relative to the inlet is denoted as l. The contact angle of fluid 1 at the walls is denoted using θ_w, while the tapering angle is denoted using β. Fluids 1 and 2 are colored red and blue, respectively. The viscous and surface-tension forces are represented using oppositely directed arrows.

In this regard, a computational fluid dynamics (CFD) model is developed by coupling the continuity and Navier–Stokes equations with the level-set method. Our motivations for adopting this approach are to resolve the interface position (global dynamics) between the displacing and displaced fluids at various time intervals and the local physics, such as interface shape, viscous dissipation, pressure profiles, and transient forces.^59,60 Moreover, CFD helps to overcome limitations of analytical models, such as the inability to capture transient interface deformation, pressure and force resolution, and the static contact-angle assumption. Initially, a uniform contact angle is applied along the microchannel walls to quantitatively demonstrate the benefits of the tapered geometry using axial velocity profiles. Subsequently, the impact of non-uniform contact-angle variations along the microchannel walls on capillary filling dynamics is investigated. Three types of contact-angle variations, namely stepwise, linear, and quadratic, are analyzed, focusing on scenarios where the contact angle progressively decreases. Finally, a gradually increasing contact angle is applied to halt fluid motion at a specific point within the capillary. This ability to stop fluid motion on demand at desired locations provides a powerful mechanism for achieving precise control over capillary filling flow. The findings of this work provide valuable insights into tailoring capillary filling dynamics in tapered capillaries. These insights pave the way for precise and on-demand fluid motion control in microfluidic systems. Such control enhances the performance of passive microfluidic pumps, enabling more efficient and reliable fluid transport. This work bridges classical capillarity with modern microfluidic engineering, offering a fresh perspective for scalable, energy-efficient fluid manipulation by unifying geometric and surface-chemical strategies.

2 Problem formulation

In view of the background earlier discussed, a tapered microchannel is considered to model the capillary filling, as shown in Fig. 1. The length of the microchannel is labeled L, with the inlet and outlet radii represented by R₁ and R₃, respectively. The system is modeled in a 2-D axisymmetric cylindrical coordinate framework, with the longitudinal central axis serving as the axis of symmetry in order to simplify the analysis while preserving the essential physics. The position of the interface relative to the inlet of the microchannel is shown as l. The contact angle of fluid 1 at the microchannel walls is denoted by θ_w, and the tapering angle is denoted by β. The system is analyzed under the assumption of wetting conditions, where θ_w < 90° results in a concave meniscus at the interface. In contrast, under non-wetting conditions (θ_w > 90°), the interface becomes convex and resists forward motion; thus, only wetting scenarios are considered in this study to ensure unidirectional flow suitable for passive pumping applications. The motion of fluid 1 into fluid 2 is driven by the combined action of surface tension forces and pressure gradients induced by the tapered geometry. A capillary pressure difference arises due to the curvature of the meniscus, governed by the Young–Laplace equation, at the fluid–fluid interface. Moreover, in the tapered geometry, the curvature and thus the capillary pressure vary along the axial direction, generating a pressure gradient. Combined with adhesive forces at the channel walls, this gradient induces autonomous forward motion of fluid 1 within the microchannel. In the following sections, two cases are examined in detail: the uniform wettability case, where the value of θ_w remains constant along the entire axial length L, and the non-uniform wettability case, where θ_w varies along the length. Different variation patterns, including linear, quadratic, and stepwise changes, are investigated for the non-uniform wettability case. The densities and viscosities of fluids 1 and 2 are denoted by ρ₁, μ₁, ρ₂, and μ₂, respectively. Fluid behavior is assumed to be Newtonian, incompressible, and isothermal. The r and z directional velocities are indicated by u_r and u_z, respectively, while the net velocity vector is symbolized by u. The fluid pressure within the microchannel is denoted as p. The interfacial surface tension is given as γ. In the subsequent section, a relationship between the above-mentioned parameters is derived by considering a constant-volume assumption of the two immiscible fluids within the tapered microchannel.

3 Solution methodology

3.1 Reduced analytical model

In this section, a reduced analytical model is derived to estimate the rate of change of the fluid interface position, i.e., its velocity. Several assumptions are introduced to simplify the analysis, ensuring the integrity of the system is preserved and remains close to the real one without sacrificing the underlying physics of interest: (i) the contact-angle hysteresis, drag, and film resistance, which are often present in surface-tension-driven microchannel flows, are neglected. These forces are considered insignificant compared to the dominant frictional drag, as corroborated by prior research.⁴⁹ Consequently, contact-line dissipation, which may occur during capillary filling, is also considered negligible. This assumption is justified by the relatively high capillary numbers encountered in our simulations (Ca ∼ 10⁰–10¹). At these capillary numbers, prior work⁶¹ and recent reviews⁵⁵ indicate that the dominant contribution to energy dissipation arises from the bulk viscous flow, and the relative contribution of contact-line friction becomes subdominant. This regime is characterized by interface velocities and system sizes where the classical hydrodynamic description without considering contact-line dissipation provides a good approximation. However, we fully recognize that this simplification may not be valid at lower velocities or in other parameter regimes, particularly when the microchannel is rough.^62–64 (ii) This study does not consider the effect of velocity on the contact angle. The contact angle is assumed to be constant, as dynamic contact-angle effects are typically significant only at very high flow velocities⁶⁵ or in systems with rapid wetting/dewetting transitions,⁶⁶ which are not the focus of this work. (iii) The temperature and fluid composition are assumed to be uniform throughout the system, eliminating any Marangoni-driven flow effects. This simplification is valid for systems where thermal gradients and compositional variations are minimal,^67,68 ensuring that capillary pressure gradients remain the primary driving force. (iv) Gravitational effects are fairly negligible in the autonomous motion of the interface, as quantified by the Bond number, Bo = ρga²/γ, where ρ is the fluid density, g is gravitational acceleration, a is the characteristic length (e.g., microchannel radius) and γ is the interfacial tension. For oil–water systems in microchannels, the Bond number is on the order of 10⁻⁶, confirming the dominance of capillary forces.⁴⁹

The system is modeled using two control volumes, each containing a single-phase incompressible fluid, as depicted in Fig. 2. The interfacial tension force, acting at the boundary separating the two control volumes, is incorporated implicitly through the boundary conditions, thereby eliminating the need for its explicit inclusion in the governing equations. The motion of the interface is governed by a balance between the capillary pressure gradient, induced by surface tension, and the viscous resistance exerted by both fluids. These dynamics are described by the incompressible Navier–Stokes equations, which account for the interplay of these forces. The governing motion of the fluids is described by the incompressible N–S equations, given in compact vector form as

∇·u = 0, (v2)

where u = (u_r, u_z) is the velocity vector in cylindrical coordinates (r, z).

Fig. 2 Tapered microchannel's exaggerated control volumes (red and blue regions) depicted for analytical modeling. The individual dimension R₂ denotes the radius of the microchannel at the interface of the two fluids.

Incorporating the above assumptions, the governing equations in radial and axial forms can thus be described as

The incompressible Navier–Stokes eqn (1)and (2) are simplified by assuming that the radial component of velocity is negligible, u_r = 0. The additional term −u_r/r² in eqn (2) arises from the vector Laplacian in cylindrical coordinates. Furthermore, extensional dissipation is neglected, as characterized by ∂u_z/∂z = 0, implying that the axial velocity component u_z is uniform along the length of the microchannel. This assumption is validated by prior studies for tapered microchannels where the channel length substantially exceeds the average diameter and the taper angle is sufficiently small.^48,49 After incorporating these assumptions at a specific time instant during the flow, eqn (2) reduces to the following expression:

It is important to note here that under the applied assumptions, the left-hand side (LHS) terms in eqn (1)and (2) become zero, indicating that inertial effects are entirely neglected in deriving the reduced analytical model. This simplification is supported by the high Ohnesorge number, defined as , where R_avg is the average radius of the microchannel. The small value of R_avg results in a large Oh, indicating the dominance of viscous forces over inertial and cohesive forces. This assumption is further corroborated by experimental studies on surface-tension-driven flows, which demonstrate that inertial effects are significant only during a transient initial phase following flow onset.^48,49 Beyond this initial phase, viscous forces dominate the flow dynamics, and the system transitions into a viscous-dominated regime.

In the case of a fully developed flow, the axial velocity u_z = U(1 − r²/R²), where U and R are the mean velocity of the fluid and the radius of the microchannel at a distance z from the inlet of the microchannel, respectively. Upon substituting the value of u_z in eqn (3), the following expression is obtained:

It is evident that all the terms in the above equation are functions of z. Therefore, upon integrating this equation across the two control volumes and then summing the results, we obtain the following expression:

Assuming a constant flow rate at a given time and expressing U = Q(t)/A(z), where A(z) = πR²(z) is the local cross-sectional area, we can further simplify in the following manner to obtain the expression below:

The interface velocity is related to the volumetric flow rate through the expression dl/dt = Q(t)/πR₂². Since the flow is driven purely by the interfacial surface tension force, the boundary pressures at the inlet and outlet are set to zero, i.e., p₀ = p_L = 0. The pressure jump across the interface, governed by the Young–Laplace equation, is given by p_l⁺ − p_l⁻ = 2γ cos(θ_w − β)/R₂. Substituting these relations into eqn (6) and eliminating Q(t), we arrive at the following expression for the rate of change of the interface position:

The above analytical equation is verified in the asymptotic case for a straight microchannel with a uniform radius. In the case of uniform radius (R = R₁ = R₂ = R₃) and a water–air system, β = 0 and μ₂ = 0. Hence, eqn (7) reduces to

Integrating eqn (8) yields , which is the classical Lucas–Washburn equation for describing the surface-tension-driven movement of water in an air-filled microchannel.⁴⁷ Therefore, the correctness of our analytical model is established.

On the other hand, for the tapered case, the microchannel radius changes along its length. At any given position within the tapered tube, the radius can be expressed as R = R₁ − z tan β, where β = tan⁻¹(L/(R₁ − R₃)). Substituting this expression into eqn (7) and performing the integration yields the final equation as follows:

Eqn (9) describes the rate of change of the interface position or the velocity of the interface (u_interface = dl/dt). The velocity of the fluid at the inlet and the outlet can be obtained using the continuity equation, utilizing the equation for the velocity of the interface:

The final eqn (10) is employed for estimating the velocity of the fluid at the inlet and the outlet of the microchannel from the interface velocity.

Here, the viscous dissipation losses may offer elevated resistances resulting in slowing down the fluid motion along the way, which necessitates quantifying them. Therefore, the viscous dissipation is also estimated to explain capillary filling dynamics in different scenarios. The local viscous dissipation per unit volume, denoted as Φ, is given by⁶⁹

The total viscous dissipation E_diss(t) at any time t is calculated by integrating Φ over the entire domain Ω, wherein the factor 2πr accounts for the axisymmetric volume element in cylindrical coordinates and the expression for the same is read as

3.2 Numerical solution

The numerical solution is obtained by solving the partial differential equations (PDEs) using CFD techniques. As the present study focuses on modeling the multiphase flow of two immiscible, Newtonian, incompressible, and isothermal fluids, the continuity equation and the Navier–Stokes equation are applied to describe the system. These governing equations are presented as follows:

∇·u = 0, (13)

where F_st denotes the force due to surface tension. The level-set method, combined with a phase-field model, is employed to simulate multiphase flow in the microchannel. In this approach, the phase-field model utilizes a phase-field variable, ϕ, which represents the fluid regions and their interface, enabling smooth tracking of the multiphase dynamics. This variable ϕ assigns values of 1 for fluid 1, 0 for fluid 2, and 0.5 at the interface, marking the spatial boundary between the two fluids. The spatial distribution of ϕ evolves over time according to the level-set equation, accurately capturing the interface dynamics and fluid boundaries throughout the microchannel flow field. The level-set equation is expressed as⁷⁰

The LHS of eqn (15) represents the material derivative showing the advection of mass with the local velocity field. The right-hand side includes a regularization term with two components: (1) the diffusion term ε_ls∇ϕ, which maintains a finite interface thickness for numerical stability, and (2) the nonlinear term −ϕ(1 − ϕ)∇ϕ/|∇ϕ|, which reinitializes ϕ to approximate a signed-distance function, preventing excessive diffusion or distortion of the interface during advection. Here, ε_ls is a small parameter controlling the artificial interface thickness, while χ is known as the reinitialization parameter. Typically, the reinitialization parameter is set to match the maximum possible fluid velocity within the microchannel, while the interface width is approximately equal to the largest mesh element size.^71,72 In this regard, the artificial interface thickness is taken as 0.7 μm, corresponding to the size of the largest mesh element. Similarly, we have adopted 100 cm s⁻¹ as the reinitialization parameter, determined after conducting multiple test runs with extremely low contact angles (0 to 5 degrees). Based on the value of ϕ, the fluid properties are approximated using standard interpolation assumptions in the level-set method as:

ρ = ρ₁ + (ρ₂ − ρ₁)ϕ, (16)

μ = μ₁ + (μ₂ − μ₁)ϕ. (17)

The scalar function ϕ, obtained by solving the level-set transport equation (eqn (15)), is then used for tracking the fluid–fluid interface in the multiphase system and estimating the surface tension force F_st. This force, arising from interfacial tension at the boundary between two immiscible fluids, is modeled as:

F_st = γκδ(ϕ)∇ϕ + γδ(ϕ)(n_wall·n_int − cos(θ_w))n_int, (18)

where γ is the interfacial tension, δ(ϕ) is a regularized Dirac delta function, n_int is the interface unit normal, and n_wall is the wall normal. The first term, γκδ(ϕ)∇ϕ, represents the classical Laplace pressure jump due to interface curvature, where the curvature is defined as κ = ∇·n_int. This contribution drives fluid motion in response to variations in interfacial geometry, such that the interface tends to minimize its surface area.

The second term, γδ(ϕ)(n_wall·n_int − cos(θ_w))n_int, incorporates wall wettability effects. In particular, it enforces the prescribed static contact angle θ_w at the fluid–solid boundary by penalizing deviations between the actual orientation of the interface (n_wall·n_int) and the desired orientation defined by (cos(θ_w)). The regularized delta function ensures that the interfacial force F_st acts only in the vicinity of the fluid–fluid interface, thereby preventing spurious force contributions in the bulk regions.

The interface normal is obtained from the level-set function ϕ as

while the regularized Dirac delta function is expressed as

δ(ϕ) = 6|ϕ(1 − ϕ)||∇ϕ|, (20)

which smoothly distributes the surface-tension force across the finite interface thickness introduced by the level-set method.

The multiphase system is solved by coupling eqn (13)–(20) with the Navier–Stokes equations under the following conditions: zero pressure at both ends of the channel (p(r, z = 0) = p(r, z = L) = 0), no slip at the walls, and a quiescent initial velocity field. The phase-field variable is initialized such that fluid 1 (olive oil) enters from the inlet (ϕ = 1 at z = 0), while fluid 2 (water) fills the remainder of the channel (ϕ = 0 for z > 0). For uniform wettability, θ_w is constant, whereas for non-uniform wettability, it varies along z.

The primary fluid pair considered is olive oil (fluid 1) and water (fluid 2), with ρ₁ = 920 kg m⁻³, μ₁ = 0.083 Pa s, ρ₂ = 1000 kg m⁻³, μ₂ = 0.001 Pa s, and γ = 0.023 N m⁻¹. Additional cases involving olive oil–seawater, sunflower oil–water, and sunflower oil–seawater are examined in the SI to assess the role of fluid properties on capillary filling dynamics.

The coupled partial differential equations (PDEs) (13)–(20) are solved in the commercial CFD software COMSOL Multiphysics 5.3a. The software employs the finite-element method to solve the non-linear PDEs. The PDEs are discretized using the Galerkin least-squares (GLS) method, resulting in a set of algebraic equations. The velocity and pressure fields were discretized using Taylor–Hood elements, with second-order (P₂) interpolation for velocity and first-order (P₁) interpolation for pressure. For the phase-field variable, first-order (P₁) elements were employed to capture the interface evolution. Within COMSOL, the laminar-flow and level-set modules are specifically utilized to model and simulate the problem. In order to track the fluid–fluid interface, the position is determined at each time step by identifying the axial coordinate x and corresponding y-coordinate along the central axis of the microchannel where the phase-field variable ϕ = 0.5. The interface velocity is then extracted directly from COMSOL's velocity-field data at this (x, y)-coordinate pair, using a custom MATLAB script integrated with the COMSOL model via the COMSOL LiveLink interface. A parallel direct sparse solver (PARDISO) is employed for transient numerical computations. All the initial and boundary conditions discussed above are incorporated into the model. The maximum relative error between successive iterations is maintained at 0.005 to ensure solution convergence. The computational domain is discretized into 12 860 grid elements, determined through a grid independence study. Additional details of this study, including the mesh resolution and the distribution of grid element sizes, are provided in the SI. An adaptive time-stepping method, with a maximum step size of 10⁻² ms, is used for the transient solver. The time-step selection method and convergence tolerance are discussed in the SI. The temporal discretization is performed using the second-order backward difference formula (BDF) for enhanced accuracy.

3.3 Validation

Fig. 3(a)–(c) compare the analytical values of u_interface, u_inlet, and u_outlet (derived from eqn (9) and (10)) with the numerical results for capillary filling of olive oil in a water-filled tapered microchannel having uniform wall wettability over 38 ms. The numerical results closely align with the analytical predictions, yielding time-averaged root-mean-square deviations of 0.257 cm s⁻¹, 0.229 cm s⁻¹, and 0.225 cm s⁻¹ for u_interface, u_inlet, and u_outlet, respectively. This agreement confirms the accuracy of the present numerical model in reproducing the theoretical dynamics. The model is further assessed against experimental data from the literature,⁷³ where capillary filling of shale oil into an air-filled channel (dimensions are given in the caption) was studied. Interestingly, the minor dip at the final instant in Fig. 3(a) and (c) results from interface interaction with the outlet boundary as the front exits and is a finite-domain artifact, not a physical or numerical instability or sudden interfacial change. Rather, it is a result of the transition of the advancing interface very close to the microchannel outlet in the final simulation moments. A similar configuration is simulated here, and Fig. 3(d) shows the deviation of the simulated interface position l from the reported experimental values. The numerical model predicts a faster filling than observed experimentally. This deviation arises primarily due to additional resistance associated with dynamic contact angles, possible contact-line dissipation, and enhanced liquid–wall interactions, which are not captured by the present model.⁵⁵ Nonetheless, the numerical results show reasonable agreement with both analytical and experimental benchmarks, highlighting the capability of the model to accurately capture capillary filling dynamics under uniform wettability and its potential applicability to cases involving non-uniform contact angle distributions.

Fig. 3 Comparison of capillary filling dynamics from numerical, analytical, and experimental results. Plots (a)–(c) present the deviation of the numerical solutions for u_interface, u_inlet, and u_outlet from their corresponding analytical predictions given in eqn (9) and (10). The values of R₁, R₃, L, and θ_w are considered as 7.5 μm, 2.5 μm, 250 μm, and 45°, respectively. Plot (d) shows the deviation of the interface position l from experimental data reported in ref. 73. The values of R₁, R₃, L, and θ_w are considered as 0.078 μm, 0.078 μm, 100 μm, and 0°, respectively.

4 Results and discussion

4.1 Uniform wall wettability

This section discusses the results of microfluidic capillary filling within a tapered microchannel with a uniform contact angle. Fig. 4(a) presents the transient contour plots of ϕ during the filling of olive oil in the tapered microchannel with a uniform wall contact angle of 45°. The observations reveal that olive oil progressively displaces water over time, driven by the acute contact angle (wetting surface) and the additional pressure drop induced by the tapering effect. Fig. 4(b)–(d) depict the temporal variation of u_interface, u_inlet, u_outlet, l, and E_diss during the filling process. Furthermore, plot (c) shows the power-law relation between l and t, and similar curve fittings are performed in the following sections. The velocity is initially high and gradually decreases over time, with a similar trend observed for the inlet and outlet velocities. As part of the validation, the interface velocity aligns with the inlet velocity at the start and the outlet velocity at the conclusion of the capillary filling. The observed behavior of the interface velocity aligns with the findings from previous experimental^33,73 and numerical studies^48,49 on capillary filling. The initially high velocity of the interface is attributed to the tapering effect, which induces an additional pressure drop within the microchannel. The subsequent deceleration sets in due to viscous dissipation. This effect intensifies as the highly viscous olive oil advances into the narrowing regions of the tapered microchannel, where the local flow resistance becomes more significant. In such confined geometries, the hydraulic resistance scales as ∼1/R⁴, and the local viscous dissipation grows with increasing length of the fluid column and decreasing radius. As a result, shear stresses intensify due to the interaction of the high-viscosity fluid with the channel walls, thereby amplifying the overall flow resistance.³³ Although local viscous dissipation intensifies in the narrowing sections due to increased flow resistance and shear stresses, the total viscous dissipation across the channel gradually decreases, mirroring the declining interface velocity. This reduction arises because the overall capillary driving force weakens as the interface advances, lowering the flow velocity and hence the energy dissipated through viscous effects. Consequently, the interface position exhibits a near-parabolic variation over time (∝ t^0.71), highlighting the progressive decline in fluid velocity as the interface advances through the microchannel. This behavior deviates from the classical ∝ t^0.5 scaling observed in straight microchannels, where the tapering is absent.

Fig. 4 Capillary filling of olive oil in a water-filled tapered microchannel with a uniform wall contact angle. Plot (a) depicts the ϕ contours at various time intervals for t = 1–24 ms, with a zoomed view of the interface. Plot (b) illustrates the temporal variation of u_interface, u_inlet, and u_outlet, while plot (c) shows the temporal variation of l over the same time range. The colored spheres show the curve fitting in order to establish the power law relationship between l and t. Plot (d) shows the temporal evolution of E_diss over the entire simulation period. The values of R₁, R₃, L, and θ_w are considered as 10 μm, 2.5 μm, 200 μm, and 45°, respectively.

Fig. 5 plots the axial pressure variation along the central longitudinal axis of the microchannel at different time intervals during the filling process. A distinct pressure dip is observed at the interface location, followed by a pressure increase across the interface. This dip is responsible for driving the forward motion of the fluid within the microchannel, countering the opposing viscous forces. As time progresses, the position of the pressure dip advances, reflecting the movement of the fluid interface. The advantages of tapering are quite evident in this plot, as the magnitude of the negative pressure increases with the progression of the interface. This increasingly pronounced dip is a direct consequence of the rising capillary pressure gradient as the microchannel radius decreases (Δp(z) ∝ 1/R(z)). This amplified pressure dip effectively counteracts the growing viscous dissipation. This case is also compared for the further examination of a straight microchannel with a constant radius, and the results are provided in the SI, highlighting why tapered microchannels achieve superior performance. Fig. 6(a) and (b) plot the temporal variation of u_interface and l for different wall contact angles during the filling process. It is evident that the interface moves more slowly as the contact angle becomes less acute and moves more rapidly for sharper angles. This behavior arises from the F_st, whose horizontal component diminishes with increasing contact angle, despite the overall magnitude of the force remaining unchanged.

Fig. 5 Pressure variation along the central longitudinal axis of the microchannel during capillary filling with a uniform wall contact angle, shown at different time intervals. The values of R₁, R₃, L, and θ_w are considered as 10 μm, 2.5 μm, 200 μm, and 45°, respectively.

Fig. 6 Temporal variation of (a) u_interface and (b) l for capillary filling under different wall contact angles of 15°, 30°, 45°, 60°, and 75°. The values of R₁, R₃, and L are considered as 10 μm, 2.5 μm, and 200 μm, respectively.

A common trend observed across all cases is the reduction in initial velocity due to viscous dissipation as the interface progresses through the microchannel, eventually leading to saturation in fluid movement. In the following section, a varying contact angle, increasing along the microchannel length, is introduced to mitigate the effects of viscous dissipation.

4.2 Non-uniform wall wettability

This section examines the effect of non-uniform contact-angle variations along the tapered microchannel walls on capillary filling dynamics. Three types of contact-angle variations, including stepwise, linear, and quadratic, are analyzed, where the contact angle progressively decreases from 75° to 15° along the microchannel length. This reduction enhances the wall's wetting nature, facilitating capillary-driven fluid motion. Additionally, a reverse scenario is explored, where the contact angle increases from 15° to 75° along the microchannel. This configuration transforms the initially wetting surface near the inlet into a non-wetting surface near the outlet, leading to a gradual reduction in fluid velocity and eventually halting the flow. By tailoring the wall contact-angle variations, the findings demonstrate how the filling dynamics of olive oil within the microchannel can be effectively controlled and optimized. These contact-angle variations are experimentally achievable through various established surface-chemistry techniques. These include localized deposition of hydrophilic or hydrophobic coatings through plasma treatment⁷⁴ or UV-ozone treatment⁷⁵ to adjust the surface energy and achieve non-uniform contact angles; photolithography or mask-assisted plasma exposure to create spatially varying wettability;⁷⁶ dip-coating with controlled withdrawal in surfactant or silane solutions for smooth gradients;⁷⁷ and electrowetting to dynamically adjust contact angles via the Young–Lippmann relation.⁷⁸ Further, it is important to note that this section focuses primarily on results for tapered microchannels. The effects of linear and quadratic contact-angle variations in a straight microchannel, where no tapering is present, are discussed in the SI and compared separately with the tapered case.

4.2.1 Increasing wall wettability.
4.2.1.1 Stepwise contact-angle variation. This section investigates the influence of stepwise variations in wall contact angle on the capillary filling dynamics within the microchannel. Specifically, the effect of discrete changes in wall wettability on fluid motion is analyzed. The axial variation of the contact angle is mathematically defined in eqn (21) and given below:

The range of contact angle variation is carefully selected to ensure that the capillary walls remain in the wetting regime throughout the microchannel. The upper limit of 75° and the lower limit of 15° are chosen to maintain sufficient wettability for effective capillary-driven flow while avoiding non-wetting conditions that could hinder fluid motion.

Fig. 7(a) presents the ϕ contours at different time intervals, while Fig. 7(b) illustrates the axial variation of the wall contact angle within the microchannel. It is important to note that zoomed-in views of the contact line are not shown in this and subsequent figures, as the interface thickness remains constant, as pre-defined in the numerical model formulation. Only the shape of the meniscus changes in response to variations in wall wettability. A key observation in this case is the increasing concavity of the interface meniscus as it advances along the microchannel. This is a direct consequence of the stepwise reduction in the contact angle, which progressively enhances the wettability of the microchannel walls and promotes more rapid fluid displacement. Unlike the uniform-contact-angle case, where viscous dissipation causes significant slowing of the interface motion near the outlet, the stepwise variation in wettability sustains interface movement along the microchannel. This is achieved by counteracting the viscous resistance through successive increases in capillary pressure gradient associated with each stepwise decrease in contact angle. The combination of geometric tapering and enhanced wettability yields a progressively increasing capillary pressure gradient, which maintains forward interface motion. Further, in Fig. 7(c)–(e), the temporal evolution of u_interface, u_inlet, u_outlet, l, and E_diss are plotted to quantitatively describe the interface motion. The interface velocity exhibits step-like increases corresponding to each stepwise decrease in contact angle, reflecting the five distinct steps in the contact angle variation. Also, the interface position follows a temporal scaling of ∝ t^1.05, indicating an enhanced filling rate in this configuration. This scaling behavior reinforces the role of stepwise wettability combined with geometric tapering in accelerating interface motion beyond the classical capillary dynamics observed in uniformly wetting systems. The total viscous dissipation shows a marked increase over time, closely following the stepwise acceleration pattern of the interface. Each stepwise drop in contact angle raises the capillary pressure, momentarily boosting the driving force and increasing shear rates, which amplify losses due to viscous resistance.

Fig. 7 Visualization of capillary filling of olive oil in a water-filled tapered microchannel with a stepwise varying wall contact angle across 5 steps. Plot (a) depicts the ϕ contours at various time intervals for t = 0.5–28 ms. Plot (b) displays the axial variation of θ_w along the microchannel length. Plot (c) illustrates the temporal variation of u_interface, u_inlet, and u_outlet, while plot (d) shows the temporal variation of l over the same time range. The colored spheres show the curve fitting in order to establish the power law relationship between l and t. Plot (e) shows the temporal evolution of E_diss over the entire simulation period. Plot (f) presents the pressure variation along the central longitudinal axis of the microchannel at different time intervals. The values of R₁, R₃, and L are considered as 10 μm, 2.5 μm, and 200 μm, respectively.

Fig. 7(f) shows the axial pressure distribution along the central longitudinal axis of the microchannel at various time intervals. Similar to the case of a uniform contact angle, a distinct pressure drop is observed at the interface, followed by a subsequent pressure rise across it. However, in this scenario, the surface-tension force increases due to a simultaneous reduction in radius and contact angle, leading to a more pronounced pressure drop compared to the case of a uniform contact angle. Additionally, a more negative pressure value is observed as the interface progresses towards the outlet of the microchannel, highlighting the advantages of employing a non-uniform, decreasing contact angle variation along with a tapered geometry. This analysis is extended in the next section, where a linear variation in the contact angle is applied to explore its impact on capillary filling dynamics.

4.2.1.2 Linear contact-angle variation. This section introduces a continuous variation in the wall contact angle, represented as an infinite series of steps, effectively modeled as a linear gradient along the microchannel walls. The applied axial variation of the contact angle is mathematically defined as

θ_w(z) = −0.3z + 75, (22)

where the units of z and θ_w are μm and degrees, respectively. Similar to the previous case, the expression for θ_w is carefully chosen such that the contact angle is 75° at the inlet and 15° at the outlet of the microchannel. This ensures that the walls remain in the wetting regime throughout the capillary.

Fig. 8(a) presents the ϕ contours at different time intervals, while Fig. 8(b) shows the axial variation of the wall contact angle in the microchannel. Similar to the previous case of stepwise variation, the interface meniscus becomes increasingly concave as it progresses through the microchannel, a result of the gradual enhancement in the wetting nature of the walls. Contrary to the uniform contact angle case, where the interface motion slows down near the outlet, the linear variation in wettability ensures sustained interface progression along the microchannel. The gradual reduction in contact angle results in a continuously increasing capillary pressure gradient, which effectively balances the growing viscous resistance imposed by the narrowing geometry. The linear variation in the contact angle offers a notable advantage by ensuring that the interface advances almost equal distances over equal time intervals, as demonstrated by the contour plots. Additionally, a slightly higher velocity is observed in the microchannel compared to the previous case. This demonstrates the effectiveness of the linear variation in maintaining more consistent and controlled capillary filling dynamics, particularly when compared to the stepwise variation. Fig. 8(c)–(e) present the temporal evolution of u_interface, u_inlet, u_outlet, l, and E_diss providing a quantitative analysis of the interface motion within the microchannel. Notably, the stepwise jumps in the interface velocity plot, as observed in the previous case, are smoothed out under the linear variation in contact angle, resulting in a nearly constant velocity throughout the capillary, counteracting the effects of viscous dissipation. Furthermore, the interface position exhibits a temporal scaling of ∝ t^0.93, approaching a near-linear trend with a slope close to π/4, which is indicative of uniform filling. This is a marked improvement over the stepwise case, promoting steadier interface advancement and reducing fluctuations in filling rate. Moreover, similar to the previous cases, the total viscous dissipation increases almost linearly as the contact angle linearly drops, increasing capillary pressure and also viscous losses as the interface progresses in the microchannel.

Fig. 8 Linearly varying wall contact angle. Plot (a) depicts the ϕ contours at various time intervals for t = 0.5–25 ms. Plot (b) displays the axial variation of θ_w along the microchannel length. Plot (c) illustrates the temporal variation of u_interface, u_inlet, and u_outlet, while plot (d) shows the temporal variation of l over the same time range. The colored spheres show the curve fitting in order to establish the power law relationship between l and t. Plot (e) shows the temporal evolution of E_diss over the entire simulation period. Plot (f) presents the pressure variation along the central longitudinal axis of the microchannel at different time intervals. The values of R₁, R₃, and L are considered as 10 μm, 2.5 μm, and 200 μm, respectively.

Fig. 8(f) plots the axial pressure distribution along the central longitudinal axis of the microchannel at various time intervals. As observed in the case of non-uniform stepwise contact-angle variation, the magnitude of the most negative pressure intensifies as the interface advances. In this case, slightly more negative pressure values are recorded over time as the interface approaches the microchannel outlet. This behavior is attributed to the linear contact-angle variation, which leads to a continuous increase in surface tension force in contrast to the stepwise increments observed previously, effectively mitigating viscous dissipation. The following section introduces a quadratic variation in the contact angle, and its impact on the interface dynamics is analyzed in detail.

4.2.1.3 Quadratic contact-angle variation. This section explores the influence of a quadratic variation in the wall contact angle on capillary filling dynamics within the microchannel. While the linear variation case typically yields a nearly constant flow velocity, introducing a quadratic variation aims to determine whether it can further increase the fluid velocity as olive oil advances through the microchannel, thereby reducing the time required to traverse from the inlet to the outlet. The goal is to leverage nonlinear surface tension changes to accelerate fluid motion and improve capillary-driven transport efficiency. The applied axial variation of the contact angle is mathematically expressed aswhere z and θ_w are given in μm and degrees, respectively. This quadratic function is carefully selected to ensure that the contact angle is 75° at the inlet (z = 0 μm) and 15° at the outlet (z = 200 μm), keeping the walls within the wetting regime throughout the microchannel. Additionally, the function is constructed so that its tangent at the midpoint (z = 150 μm) is parallel to the linear variation defined in eqn (22). This ensures a smooth transition between linear and quadratic variations, enabling a systematic comparison of their effects on capillary filling dynamics.

Fig. 9(a) presents the ϕ contours at different time intervals, while Fig. 9(b) shows the axial variation of the wall contact angle in the microchannel. A notable increase in fluid velocity is observed within the microchannel, driven by the quadratic variation of the contact angle along the wall. This variation enables the interface to traverse from the inlet to the outlet in 20 ms. In contrast to the linear wettability case, where the capillary pressure gradient increased uniformly and maintained nearly constant interface velocity, the quadratic variation introduces a non-monotonic profile. Initially, the increasing wettability enhances the capillary pressure gradient and supports interface advancement. However, as the interface progresses toward the outlet, the contact angle begins to increase again, as shown in Fig. 9(b), leading to a reduction in the capillary driving force. This results in a slight decline in interface velocity, as the growing viscous resistance is no longer fully balanced by the surface tension force. The deceleration observed toward the end of the channel reflects this competing influence between viscous dissipation and the spatially varying capillary pressure gradient. This behavior is quantitatively shown in Fig. 9(c)–(e), which plot the temporal evolution of u_interface, u_inlet, u_outlet, l, and E_diss. The interface notably moves faster compared to the two previous cases of non-uniform wettability, where linear and stepwise variations in the contact angle were studied. Moreover, the interface position exhibits a temporal scaling of ∝ t^0.84, reflecting a parabolic trend that lies between the linear and stepwise wettability cases. This intermediate scaling behavior indicates a more balanced interplay between increasing viscous resistance and the spatially varying capillary pressure gradient. The total viscous dissipation closely follows the velocity profile. It increases as capillary forces accelerate the flow, reflecting greater work done against viscous resistance, and decreases as the interface slows down.

Fig. 9 Quadratically varying wall contact angle. Plot (a) depicts the ϕ contours at various time intervals for t = 0.5–24 ms. Plot (b) displays the axial variation of θ_w along the microchannel length. Plot (c) illustrates the temporal variation of u_interface, u_inlet, and u_outlet, while plot (d) shows the temporal variation of l over the same time range. The colored spheres show the curve fitting in order to establish the power law relationship between l and t. Plot (e) shows the temporal evolution of E_diss over the entire simulation period. Plot (f) presents the pressure variation along the central longitudinal axis of the microchannel at different time intervals. The values of R₁, R₃, and L are considered as 10 μm, 2.5 μm, and 200 μm, respectively.

Fig. 9(f) shows the axial pressure distribution along the central longitudinal axis of the microchannel at various time intervals. A similar trend was observed in previous cases. Still, the magnitude of the negative pressure at the interface position increases more rapidly with time, reaching higher values at a faster rate. This is due to the quadratic variation in the contact angle, which causes a faster rise in surface-tension force compared to other variations, effectively counteracting viscous dissipation. Therefore, based on specific requirements, the variation in the contact angle can be tailored to achieve the desired velocity of the autonomous capillary filling motion inside the microchannel.

4.3 Decreasing wall wettability

This section explores the influence of a linearly increasing wall contact angle variation on the microchannel capillary filling dynamics. Unlike the approach in the previous section, the focus here is on gradually decelerating the fluid motion within the capillary, ultimately bringing it to a halt at a specific position. The applied axial variation of the contact angle is mathematically defined as

θ_w(z) = 0.3z + 60, (24)

where the units of z and θ_w are μm and degrees, respectively. In this case, eqn (24) is chosen such that the value of θ_w increases linearly from 60° at the inlet to 120° at the outlet of the microchannel.

Fig. 10(a) presents the ϕ contours at different time intervals, while Fig. 10(b) shows the axial variation of the wall contact angle in the microchannel. A noticeable deceleration of the interface occurs as it progresses towards the midpoint of the microchannel, driven by viscous dissipation and the rising contact angle. This deceleration corresponds to the contact angle approaching 90°, transitioning the wall to a non-wetting state. Although the overall magnitude of F_st increases along the microchannel length due to the decreasing capillary radius, its horizontal component diminishes with the increasing contact angle. Consequently, the horizontal motion of the fluid slows down, eventually halting when the horizontal component of F_st reduces to zero as the contact angle nears a right angle, exacerbated by viscous dissipation. This phenomenon is quantitatively presented in Fig. 10(c)–(e), which plot the temporal evolution of u_interface, u_inlet, u_outlet, l, and E_diss. The velocities of the inlet, interface, and outlet decrease steadily, nearing zero as the interface reaches the midpoint of the microchannel. Similarly, the interface position follows a parabolic trend (∝ 0.30), plateauing as it approaches the midsection. The total viscous dissipation also diminishes over time, ultimately vanishing as the interface velocity approaches zero. The deceleration of fluid within the microchannel is quantitatively illustrated in the pressure variation plot along the central longitudinal axis (Fig. 10(f)). A pressure dip is observed at the interface position, consistent with previous cases. Initially, the magnitude of this negative pressure increases up to 10 ms due to the tapered capillary geometry, which enhances F_st. However, as the interface advances after 10 ms, the magnitude of the pressure dip diminishes, influenced by the increasing contact angle and viscous dissipation, which outweigh the effect of the decreasing capillary radius caused by tapering. These results demonstrate that a reverse linear variation in the wall contact angle effectively decelerates fluid motion within the microchannel, providing a robust mechanism to regulate and control capillary filling dynamics.

Fig. 10 Linearly increasing wall contact angle. Plot (a) depicts the ϕ contours at various time intervals for t = 1–28 ms. Plot (b) displays the axial variation of θ_w along the microchannel length. Plot (c) illustrates the temporal variation of u_interface, u_inlet, and u_outlet, while plot (d) shows the temporal variation of l over the same time range. The colored spheres show the curve fitting in order to establish the power law relationship between l and t. Plot (e) shows the temporal evolution of E_diss over the entire simulation period. Plot (f) presents the pressure variation along the central longitudinal axis of the microchannel at different time intervals. The values of R₁, R₃, and L are considered as 10 μm, 2.5 μm, and 200 μm, respectively.

5 Conclusions

This study investigates the dynamics of capillary filling of olive oil in a tapered microchannel filled with water under different wall wettability configurations. A CFD model is developed that integrates the continuity and Navier–Stokes equations with the level-set method to track the fluid–fluid interface over time. While the primary emphasis is on non-uniform wall wettability, a uniform wettability configuration is initially applied for comparative analysis. The interface dynamics, quantified by its velocity and position, as well as the pressure distribution within the microchannel, are presented graphically. In the case of uniform wall wettability, the combination of an acute contact angle and tapered geometry initially drives the olive oil forward. However, as the high-viscosity fluid advances into narrower regions of the channel, the resulting increase in local viscous dissipation leads to a gradual decline in interface velocity. This behavior manifests as a near-parabolic velocity profile, reflecting progressive deceleration despite the increasing capillary pressure gradient introduced by the tapering. These factors highlight the limitations of uniform wettability, where viscous forces impede sustained flow. To overcome this limitation, three non-uniform wettability variations including stepwise, linear, and quadratic, were applied at the walls. The stepwise variation triggers velocity spikes at sudden contact-angle drops, countering slowdown and maintaining flow. The linear variation delivers a steady velocity, promoting uniform interface movement and a gradual negative pressure rise, ideal for consistent capillary filling. The quadratic variation markedly speeds up flow, overcoming viscous dissipation with a rapid interface advance and a very high pressure dip, perfect for fast fluid transport. Finally, a reverse linear variation in contact angles is introduced, transitioning the wall from wetting to non-wetting along the microchannel length. This configuration effectively decelerates the fluid motion, eventually halting the interface at a desired location. These findings underscore the significant potential of integrating geometric modifications with tailored wettability to regulate autonomous capillary filling dynamics in microfluidic systems. Such control enables precise modulation of fluid motion, including acceleration, stabilization, deceleration, and halting, which can be applied to enhance the performance of passive microfluidic pumps and other capillary-driven systems. Future research can look into how outside factors like temperature changes, electric forces, or magnetic fields interact with changes in wettability, providing a deeper understanding of how this method can be used in real-world situations and scaled up. In addition, the incorporation of advanced numerical strategies, such as moving or adaptive mesh techniques, could further improve computational efficiency and accuracy in capturing interface dynamics under these complex conditions.

Author contributions

Soumadip Das: conceptualization (equal), methodology (equal), software (lead), validation (equal), visualization (equal), formal analysis (equal), writing – original draft (equal), and writing – review and editing (equal). Vinod B. Vanarse: conceptualization (equal), methodology (equal), software, validation (equal), visualization (equal), formal analysis (equal), writing – original draft (equal), writing – review and editing (equal), and supervision (equal). Omkar S. Deshmukh: conceptualization (equal), methodology (equal), writing – original draft (equal), and supervision (equal).

Conflicts of interest

There are no conflicts to declare.

Data availability

All the data utilized in this work are incorporated in the main manuscript. The required details to support results and discussion are given in the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm00771b.

Acknowledgements

We gratefully acknowledge the Department of Chemical Engineering, IIT Guwahati, for providing computational facilities. Vinod and Soumadip extend their sincere gratitude to Prof. Dipankar Bandyopadhyay, Department of Chemical Engineering, IIT Guwahati, for mentorship and invaluable lessons imparted in our previous projects toward independent research that helped us to successfully complete this project. We also thank Dr Siddhartha Mukherjee, Assistant Professor, and Ms Nidhi Shrivastava from the Department of Chemical Engineering, IIT Guwahati, for insightful discussions. Finally, we appreciate the anonymous reviewers whose thoughtful comments and suggestions greatly improved the quality and clarity of this work.

References

1. K. M. Wisdom, J. A. Watson, X. Qu, F. Liu, G. S. Watson and C.-H. Chen, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 7992–7997.
2. F. Sartor and Á. T. Kovács, mBio, 2022, 13, e0170322.
3. J. S. Guasto, R. Rusconi and R. Stocker, Annu. Rev. Fluid Mech., 2012, 44, 373–400.
4. E. Lauga, Annu. Rev. Fluid Mech., 2016, 48, 105–130.
5. X. Zhu, Z. Yin and Z. Ni, Microfluid. Nanofluid., 2011, 12, 529–544.
6. M. Prakash, D. Quéré and J. W. M. Bush, Science, 2008, 320, 931–934.
7. J. W. M. Bush, F. Peaudecerf, M. Prakash and D. Quéré, Adv. Colloid Interface Sci., 2010, 161, 10–14.
8. C. Li, H. Dai, C. Gao, T. Wang, Z. Dong and L. Jiang, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 12704–12709.
9. Y. Zheng, X. Jia, H. Shi, W. Xu, Z. Tan, Y. Cao, Z. Dan and Z. Dai, Microchem. J., 2024, 206, 111579.
10. S. Augustine, M. V. Chinnamani, C. W. Mun, J.-Y. Shin, T. Q. Trung, S. J. Hong, L. T. N. Huyen, E. H. Lee, S. H. Lee, J.-J. Rha, S. Jung, Y. Lee, S.-G. Park and N.-E. Lee, Biosens. Bioelectron., 2024, 248, 115987.
11. A. Olanrewaju, M. Beaugrand, M. Yafia and D. Juncker, Lab Chip, 2018, 18, 2323–2347.
12. N. S. Lynn and D. S. Dandy, Lab Chip, 2009, 9, 3422–3429.
13. N. Limjanthong, Y. Tohbaru, T. Okamoto, R. Okajima, Y. Kusama, H. Kojima, A. Fujimura, T. Miyazaki, T. Kanamori, S. Sugiura and K. Ohnuma, J. Biosci. Bioeng., 2023, 135, 151–159.
14. N. M. Reis, S. H. Needs, S. M. Jegouic, K. K. Gill, S. Sirivisoot, S. Howard, J. Kempe, S. Bola, K. Al-Hakeem, I. M. Jones, T. Prommool, P. Luangaram, P. Avirutnan, C. Puttikhunt and A. D. Edwards, ACS Sens., 2021, 6, 4338–4348.
15. I. Lunati and D. Or, Phys. Fluids, 2009, 21, 052003.
16. J. J. Davis, M. Padalino, A. S. Kaplitz, G. Murray, S. W. Foster, J. Maturano and J. P. Grinias, Anal. Chim. Acta, 2021, 1151, 338230.
17. Y. Bai, X. Yu, X. Han, Y. Liu and G. Li, Sens. Actuators, A, 2024, 377, 115683.
18. T. N. A. Vo, P.-C. Chen, P.-S. Chen, Y.-C. Jair and Y.-H. Wu, Sens. Actuators, A, 2024, 379, 115845.
19. P. Bacchin, J. Leng and J.-B. Salmon, Chem. Rev., 2022, 122, 6938–6985.
20. T. N.-D. Cao, C.-C. Chang, H. Mukhtar, Q. Sun, Y. Li and C.-P. Yu, Environ. Res., 2022, 215, 114347.
21. V. B. Vanarse, S. Thakur, A. Ghosh, P. R. Parmar and D. Bandyopadhyay, Phys. Fluids, 2024, 36, 022115.
22. P. Wang, S. Yuan, N. Yang and P. K. Oppong, J. Micromech. Microeng., 2021, 31, 093001.
23. Y.-N. Wang and L.-M. Fu, Microelectron. Eng., 2018, 195, 121–138.
24. P. Azizian, J. Casals-Terré, J. Ricart and J. M. Cabot, Analyst, 2023, 148, 2657–2675.
25. R. Lucas, Kolloid-Z., 1918, 23, 15–22.
26. E. W. Washburn, Phys. Rev., 1921, 17, 273–283.
27. F. Brochard, Langmuir, 1989, 5, 432–438.
28. M. K. Chaudhury and G. M. Whitesides, Science, 1992, 256, 1539–1541.
29. O. Sandre, L. Gorre-Talini, A. Ajdari, J. Prost and P. Silberzan, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 2964–2972.
30. J. Cai, Y. Chen, Y. Liu, S. Li and C. Sun, Adv. Colloid Interface Sci., 2022, 304, 102654.
31. J. McCraney, M. Weislogel and P. Steen, npj Microgravity, 2021, 7, 1–11.
32. R. Dangla, S. C. Kayi and C. N. Baroud, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 853–858.
33. J.-B. Gorce, I. J. Hewitt and D. Vella, Langmuir, 2016, 32, 1560–1567.
34. D. Baratian, A. Cavalli, D. V. D. Ende and F. Mugele, Soft Matter, 2015, 11, 7717–7721.
35. S. Xing, J. Jiang and T. Pan, Lab Chip, 2013, 13, 1937–1947.
36. A. K. Yetisen, M. S. Akram and C. R. Lowe, Lab Chip, 2013, 13, 2210–2251.
37. T. Tian, Y. Bi, X. Xu, Z. Zhu and C. Yang, Anal. Methods, 2018, 10, 3567–3581.
38. C. A. Holstein, A. Chevalier, S. Bennett, C. E. Anderson, K. Keniston, C. Olsen, B. Li, B. Bales, D. R. Moore, E. Fu, D. Baker and P. Yager, Anal. Bioanal. Chem., 2016, 408, 1335–1346.
39. A. Nilghaz, D. R. Ballerini and W. Shen, Biomicrofluidics, 2013, 7, 051501.
40. R. Safavieh, G. Z. Zhou and D. Juncker, Lab Chip, 2011, 11, 2618–2624.
41. S. B. Berry, J. J. Lee, J. Berthier, E. Berthier and A. B. Theberge, Anal. Methods, 2019, 11, 4528–4536.
42. L. Gervais, M. Hitzbleck and E. Delamarche, Biosens. Bioelectron., 2011, 27, 64–70.
43. L. Gervais and E. Delamarche, Lab Chip, 2009, 9, 3330–3337.
44. M. Zimmermann, H. Schmid, P. Hunziker and E. Delamarche, Lab Chip, 2006, 7, 119–125.
45. A. Javadi, M. Habibi, F. S. Taheri, S. Moulinet and D. Bonn, Sci. Rep., 2013, 3, 1412.
46. M. I. Khodabocus, M. Sellier and V. Nock, Adv. Math. Phys., 2016, 2016, e1234642.
47. J. Cai, T. Jin, J. Kou, S. Zou, J. Xiao and Q. Meng, Langmuir, 2021, 37, 1623–1636.
48. A. Salama, J. Fluid Mech., 2022, 947, A12.
49. A. Salama, J. Kou, B. Dawoud, M. Rady and S. El Morshedy, Colloids Surf., A, 2023, 664, 131151.
50. R. Shabani and H. J. Cho, Sens. Actuators, B, 2013, 180, 114–121.
51. S. Das and S. K. Mitra, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2013, 87, 063005.
52. D. H. Jun and S. S. Yang, Sens. Actuators, A, 2015, 234, 294–301.
53. I. Jang, H. Kang, S. Song, D. S. Dandy, B. J. Geiss and C. S. Henry, Analyst, 2021, 146, 1932–1939.
54. H. Kusumaatmaja, C. M. Pooley, S. Girardo, D. Pisignano and J. M. Yeomans, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2008, 77, 067301.
55. J. H. Snoeijer and B. Andreotti, Annu. Rev. Fluid Mech., 2013, 45, 269–292.
56. J. Cheng, Y. Tao, Y. Zhang, Z. Q. Cai, P. H. Pi, X. F. Wen, D. F. Zheng, Z. R. Yang, L. S. Lu and Y. Tang, Adv. Mater. Res., 2011, 233–235, 1152–1156.
57. E. Berthier and D. J. Beebe, Lab Chip, 2007, 7, 1475–1478.
58. A. Hooshanginejad and S. Lee, Phys. Rev. Fluids, 2022, 7, 033601.
59. A. Russo, M. Icardi, M. Elsharkawy, D. Ceglia, P. Asinari and C. M. Megaridis, Phys. Rev. Fluids, 2020, 5, 074002.
60. S. Prokopev, A. Vorobev and T. Lyubimova, Phys. Rev. E, 2019, 99, 033113.
61. M. J. Miksis, J. Geophys. Res.:Solid Earth, 1988, 93, 6624–6634.
62. E. Reyssat, J. Fluid Mech., 2014, 748, 641–662.
63. J. Delannoy, S. Lafon, Y. Koga, É. Reyssat and D. Quéré, Soft Matter, 2019, 15, 2757–2761.
64. B. Primkulov, J. Chui, A. Pahlavan, C. MacMinn and R. Juanes, Phys. Rev. Lett., 2020, 125, 174503.
65. A. Nabizadeh, H. Hassanzadeh, M. Sharifi and M. Keshavarz Moraveji, J. Mol. Liq., 2019, 292, 111457.
66. T. P. Allred, J. A. Weibel and S. V. Garimella, Int. J. Heat Mass Transfer, 2021, 172, 121167.
67. M. Gugliotti, M. S. Baptista, M. J. Politi, T. P. Silverstein and C. D. Slater, J. Chem. Educ., 2004, 81, 824.
68. A. Ebrahimi, Y. Kazemzadeh and A. Akbari, Heliyon, 2024, 10, e39919.
69. J. Dhar, S. Mukherjee, K. R. M and S. Chakraborty, EPL, 2019, 125, 14003.
70. V. H. Gada and A. Sharma, Numer. Heat Transfer, Part B, 2009, 56, 307–322.
71. M. Sussman, P. Smereka and S. Osher, J. Comput. Phys., 1994, 114, 146–159.
72. S. Xu and W. Ren, Commun. Comput. Phys., 2016, 20, 1163–1182.
73. H. Lu, Y. Xu, C. Duan, P. Jiang and R. Xu, Energy Fuels, 2022, 36, 5267–5275.
74. J. A. S. Ting, L. M. D. Rosario, H. V. Lee, H. J. Ramos and R. B. Tumlos, Surf. Coat. Technol., 2014, 259, 7–11.
75. Z. Fan, C. Zhi, L. Wu, P. Zhang, C. Feng, L. Deng, B. Yu and L. Qian, Coatings, 2019, 9, 762.
76. M. Zhu and Y. Mao, ACS Appl. Polym. Mater., 2025, 7, 2451–2458.
77. M. Faustini, D. R. Ceratti, B. Louis, M. Boudot, P.-A. Albouy, C. Boissière and D. Grosso, ACS Appl. Mater. Interfaces, 2014, 6, 17102–17110.
78. S. Adriano and M. Oliveira, Geophysics, 2003, 68, 672–676,  DOI:10.1190/1.1567237.

---

Footnote

† These authors contributed equally to this work.

---