Capillary-flow dynamics in open rectangular microchannels

Abstract Spontaneous capillary flow of liquids in narrow spaces plays a key role in a plethora of applications including lab-on-a-chip devices, heat pipes, propellant management devices in spacecrafts and flexible printed electronics manufacturing. In this work we use a combination of theory and experiment to examine capillary-flow dynamics in open rectangular microchannels, which are often found in these applications. Scanning electron microscopy and profilometry are used to highlight the complexity of the free-surface morphology. We develop a self-similar lubrication-theory-based model accounting for this complexity and compare model predictions to those from the widely used modified Lucas–Washburn model, as well as experimental observations over a wide range of channel aspect ratios $\lambda$ and equilibrium contact angles $\theta _0$. We demonstrate that for large $\lambda$ the two model predictions are indistinguishable, whereas for smaller $\lambda$ the lubrication-theory-based model agrees better with experiments. The lubrication-theory-based model is also shown to have better agreement with experiments at smaller $\theta _0$, although as $\theta _0\rightarrow {\rm \pi}/4$ it fails to account for important axial curvature contributions to the free surface and the agreement worsens. Finally, we show that the lubrication-theory-based model also quantitatively predicts the dynamics of fingers that extend ahead of the meniscus. These findings elucidate the limitations of the modified Lucas–Washburn model and demonstrate the importance of accounting for the effects of complex free-surface morphology on capillary-flow dynamics in open rectangular microchannels.


Introduction
Capillary flow is the spontaneous wicking of liquid in narrow spaces without the assistance of, or even in opposition to, external forces such as gravity. This phenomenon has been investigated since the early twentieth century and has been exploited for a diverse range of applications including lab-on-a-chip devices (Olanrewaju et al. 2018), heat pipes (Faghri 1995), propellant management devices in spacecrafts (Levine et al. 2015) and fabrication of flexible printed electronics (Cao et al. 2018;Jochem et al. 2018).
Early studies focused on understanding the physical mechanism driving spontaneous capillary flow in capillary tubes. Lucas (1918) and Washburn (1921) appear to have been the first to propose theoretical models describing the meniscus positionz m as a function of timet for flow of a Newtonian liquid in cylindrical capillaries. Lucas (1918) assumed the flow is driven by the capillary pressure gradient caused by the circular-arc meniscus front, while Washburn (1921) also included hydrostatic pressure gradients and an imposed pressure difference between the two ends of the capillary. For a horizontal capillary tube open at both ends, an analytical solutionz m = kt is obtained, commonly referred to as the Lucas-Washburn relation, wherek is known as the mobility parameter and depends on the cylinder radius, liquid viscosity, surface tension and contact angle. The mobility parameterk can be thought of as a diffusion coefficient driving the growth of the liquid interface.
Due to breakthroughs in lithographic fabrication techniques, open microchannels with various cross-sectional geometries can be fabricated easily and inexpensively, including rectangular (Yang et al. 2011;Sowers et al. 2016;Lade et al. 2018;Kolliopoulos et al. 2019), trapezoidal (Chen 2014), U-shaped (Yang et al. 2011) and V-shaped (Rye, Yost & O'Toole 1998;Mann et al. 1995;Rye, Mann & Yost 1996;Yost, Rye & Mann 1997) cross-sections. The lack of a top provides access to the inside of the channel, and has been exploited in applications such as capillary micromoulding and microfluidics. Some studies have generalized the Lucas-Washburn relation to arbitrary cross-sectional geometries (Ouali et al. 2013;Berthier, Gosselin & Berthier 2015). However, predictions of the modified Lucas-Washburn models for open capillaries have resulted in varying agreement with experiments (Yang et al. 2011;Ouali et al. 2013;Chen 2014;Sowers et al. 2016;Kolliopoulos et al. 2019). This is because the mechanism for capillary flow in open channels is more complex than that for closed channels. While for closed channels the force driving the flow is due to the pressure gradient caused by the circular-arc meniscus front, for open channels the additional free surface also contributes to driving the flow (this will be discussed in more detail when presenting figure 4).
The additional contribution of the free-surface curvature to capillary flow has been theoretically and experimentally investigated primarily for V-shaped channels (Mann et al. 1995;Romero & Yost 1996;Rye et al. 1996;Yost et al. 1997;Rye et al. 1998;Weislogel & Lichter 1998;Weislogel 2012). However, while the most widely used open-channel cross-sectional geometry is rectangular (Olanrewaju et al. 2018), previous theoretical studies have only considered capillary flow in open rectangular channels for liquids with contact angles of θ 0 = 0 • and large channel aspect ratios λ = H/W (height/width) (Tchikanda, Nilson & Griffiths 2004;Nilson et al. 2006), or reported three-dimensional simulations using the volume-of-fluid method to study the effects of gravity on capillary rise in open rectangular channels (Gurumurthy et al. 2018).
In open rectangular channels the free-surface morphology is more complex than in V-shaped channels. From the channel inlet to the meniscus front the upper meniscus spans the entire channel width. However, at the meniscus front the flow splits into the channel corners provided the equilibrium contact angle θ 0 < π/4 (Concus & Finn 1969). This splitting of the flow leads to filaments or fingers extending ahead of the meniscus front and influencing its propagation. Such a transition is not observed in V-shaped channels.
In this work we use a combination of experiment ( § 2) and theory ( § 3) to study capillary-flow dynamics in open rectangular channels. This is achieved by developing a self-similar lubrication-theory-based model ( § 3.2), and comparing model predictions to the modified Lucas-Washburn (MLW) model ( § 3.1) and complementary flow visualization experiments. We investigate the effects of the complex free-surface morphology on the flow dynamics over a wide range of channel aspect ratios λ and equilibrium contact angles θ 0 ( § 4.1) and identify limitations of the MLW model ( § 4.2). Finally, we show good agreement between lubrication-theory-based model predictions of the finger dynamics and experiments ( § 4.3).

Capillary-flow experiments
Experiments with a non-volatile liquid are used to study capillary flow in open rectangular microchannels. Flow visualization is used to track the meniscus front and a combination of scanning electron microscopy (SEM) and profilometry is used to characterize the effect of channel aspect ratio on the free-surface morphology.

Channel fabrication and materials characterization
2.1.1. Fabrication of master pattern Traditional microfabrication techniques were used to form silicon master patterns of capillary channels. A 10.2 cm diameter silicon wafer was cleaned in an oxygen asher (Technics Oxygen Asher) for 5 min with an oxygen flow of 200 SCCM and RF power of 250 W. MicroChem SU-8 2010 negative tone photoresist was spin-coated onto the wafer at 300 rpm for 5 s and 1000 rpm for 30 s, followed by edge-bead removal with MicroChem EBR PG. These coating conditions target a 20 μm layer thickness. Fabrication of capillary channels using SU-8 was chosen because it gives smother sidewalls, sharper bottom corners and a flatter channel bottom than deep reactive-ion etching. The resist was soft-baked on a hot plate at 95 • C for 4 min. The photoresist was exposed through a photomask using a Karl Suss MA6 contact mask aligner in soft contact mode for 12.5 s with a 50 μm gap to define the capillary channels. Measurement gradient marks were included in the master pattern to facilitate tracking of the capillary flow. The wafer was then baked at 95 • C for 4 min. The exposed wafer was developed in propylene glycol momomethyl ether acetate (Sigma Aldrich) and rinsed with isopropanol. The resist was then hard baked at 150 • C for 30 min and an anti-stick fluorinated monolayer was formed by placing the dried wafer in a reduced-pressure chamber with trichloro(1H,1H,2H,2H-perfluorooctyl)silane (Sigma Aldrich) vapour overnight. The resulting microchannel height was 22.5 μm, measured with a KLA Tencor P16 surface profilometer.

Substrate fabrication
Capillary channels were prepared by first casting a silicone stamp (Sylgard 184) over the master pattern, curing the stamp and then using the stamp to imprint UV-curable adhesive (Norland Products NOA68 or NOA73) as explained in Kolliopoulos et al. (2019). Briefly, the UV-curable adhesive was coated on glass slides and then the silicone stamp was pressed into the adhesive. The adhesive was solidified by exposure to 365 nm UV light (Honle UV Spot 100) at 30 mW cm −2 for 270 s. The stamp was then delaminated from the prepared capillary channels and the channels were inspected with a digital microscope for defects. Any channels with defects were not used for capillary-flow experiments. The microchannel length and height were 30 mm and 22.5 μm, respectively. Microchannel widths were 17, 25, 50, 75, 100 and 200 μm. The reservoir radius was 3 mm. A schematic of the microchannel geometry and a SEM image of a 100 μm wide and 22.5 μm deep channel are shown in figure 1.

Materials characterization
The non-volatile test liquids chosen for capillary-flow experiments included UV-curable adhesive (NOA74, Norland Products), silicone oil (DC-704, Dow Corning Corporation), mineral oil (Sigma-Aldrich) and propylene glycol (Froggy's Fog). Shear viscosity μ was measured using a stress-controlled rheometer (AR-G2, TA Instruments) with a stainless steel cone-and-plate geometry (40 mm, 2 • cone angle). Surface tension σ was measured using a Krüss DSA-30 digital tensiometer. A Krüss goniometer was used to measure equilibrium contact angles θ 0 on flat test substrates prepared in the same way as the capillary channels. Density values were obtained from the manufacturer specifications. The physical properties and equilibrium contact angles of the test liquids are shown in  Nordson EFD tip mounted above the reservoir. Sufficient liquid was deposited at the centre of the reservoir to fully fill the reservoir. After deposition, a programmed microstepping motor (Automation Direct STP-MTRD-23042RE) moved the stage assembly and ensured the liquid front remained in the field of view, allowing for visualization of longer flow distances compared with prior studies (Yang et al. 2011;Sowers et al. 2016;Lade et al. 2018;Kolliopoulos et al. 2019). A high-speed camera (Photron Fastcam-Ultima APX) with a Micro-Nikkor 105 mm lens, Nikon PN-11, Nikon PK-13 and Kenko 20 mm and 36 mm extension tubes and a Kenko N-AFD 2× Teleplus MC7 lens was used to visualize the flow at 60 frames per second. Flow was recorded until the liquid meniscus reached the end of the 30 mm long channel or until the maximum recording time of the camera (∼400 s) was reached. Experiments were conducted at ambient conditions (23 ± 1 • C). Flow visualization experiments were analysed using ImageJ software. A minimum of four trials were conducted for a given channel aspect ratio and test liquid. The meniscus-position time evolutionz m (t) was averaged over all trials and the maximum and minimumz m (t) were used for the range of experimental results. The meniscus-position time evolutioñ z m (t) for different channel aspect ratios λ using NOA74 as the test liquid can be seen in figure 2(b). Results using the other test liquids are reported in § 4.2.

Free-surface profile characterization
The following experiments were conducted to investigate the effect of channel aspect ratio λ on the free-surface morphology. A UV-curable liquid (NOA74, Norland Optical Adhesives) was deposited in the reservoir connected to the microchannel and allowed to flow along the channel length. The flow was terminated at a desired time by exposing the NOA74 to a high-intensity UV light source (Omnicure S1500A with a custom light guide) at a UV intensity of approximately 1.6 W cm −2 . The liquid was fully cured in <2 s, but the flow terminates well before full solidification, so the process essentially creates a snapshot of the free-surface profile at a given time. The position of the meniscus front and the time at which curing occurred are represented by the filled symbols in figure 2(b) for two channel aspect ratios λ.
After solidification, the free-surface profile was measured with a stylus profilometer (KLA-Tencor P16) by making repeated scans across the channel width. The samples were then coated with a conductive gold film and the region near the meniscus front was imaged with a SEM instrument (JEOL JSM-6010PLUS/LA) in secondary electron imaging mode with a sample rotated 40 • about thez axis (figure 3f ). The SEM images and profilometry scans for channels of aspect ratio λ = 0.45 and λ = 0.225 corresponding to the filled symbols in figure 2(b) are shown in figure 3.

Free-surface morphology
We investigate the effect of channel aspect ratio λ on the free-surface morphology by initially examining the profilometry scans in figures 3(a) and 3(c). It is observed that for both channel aspect ratios the liquid heighth at the centre of the channel decreases down the channel length. This decrease inh at the centre of the channel results in an increase of the free-surface curvature, causing capillary pressure gradients that drive the flow. At a certain distance down the channel,h at the centre of the channel goes to zero and the liquid splits into two filaments along the bottom corners. The filament morphology can be seen in the SEM images in figures 3(b) and 3(d).
From figures 3(a) and 3(b), it appears that the free-surface morphology can be divided into three regimes. The first is a meniscus-deformation regime (I) (or accommodation regime) where the liquid is pinned to the top of sidewalls and the top meniscus curvature increases down the length of the channel. The second is a meniscus-recession regime (II) where the liquid depins from the top of the channel wall and the meniscus begins to recede down the channel walls. The third is a corner-flow regime (III) where the liquid splits and recedes into the corners.
Examination of figures 3(c) and 3(d) also suggests the presence of three regimes. The first is the meniscus-deformation regime (I) similar to that seen in figure 3(b). However, in this case the meniscus splits into filaments prior to the liquid depinning from the top of the channel wall (see figure 3e) so that the meniscus-recession regime (II) is absent. After the splitting of the meniscus a corner-transition regime (IV) is observed, where the liquid remains pinned to the top of the channel wall. This is followed by a corner-flow regime (III) similar to that seen in figure 3 The above visualizations suggest that there is a critical channel aspect ratio λ c at which the free-surface morphology transitions from that seen in figures 3(a) and 3(b) to that seen in figures 3(c) and 3(d). This is in agreement with experimental observations of Seemann et al. (2005). In their study, polystyrene droplets were deposited on grooves with rectangular cross-sections via vapour condensation. The polystyrene droplets flowed in the grooves and were solidified by lowering the temperature of the polymer below its glass transition temperature. The solidified samples were then characterized using atomic force microscopy. The expression that Seemann et al. (2005) used for λ c was by assuming a circular upper meniscus contacting the bottom of the rectangular channel while being attached to the top of the channel walls. For NOA74 with θ 0 = 10 • we obtain λ c = 0.42 from (2.1), which is consistent with the free-surface morphology transition observed in figure 3. While the corner-transition (IV) regime has been previously observed experimentally, it has not, to the best of our knowledge, been accounted for in theoretical studies.
In the following sections we evaluate the importance of the free-surface morphology in model predictions. This is achieved by comparing two theoretical models describing capillary flow, where one accounts for the complexity of the free-surface morphology ( § 3.2) whereas the other assumes that it is flat ( § 3.1). 911 A32-7

Mathematical modelling
Here, we describe two mathematical models for capillary flow of a non-volatile, isothermal Newtonian liquid in an open rectangular channel in contact with an ambient passive gas. We consider a liquid of density ρ, viscosity μ, surface tension σ and equilibrium contact angle θ 0 . The open rectangular channel has width W, height H and length L. In this work, we use the notationf to denote the dimensional version of a variable f .

Modified Lucas-Washburn model
In this model, the flow is assumed to be driven by the capillary pressure gradient caused by the circular-arc meniscus front, while viscous forces resist the flow and inertial and gravitational forces are neglected. The capillary driving force is obtained by assuming a flat upper liquid-air interface and a circular-arc meniscus front governed by fluid statics, while the viscous force is obtained by assuming a fully developed parallel flow. Through conservation of linear momentum in the axial direction an analytical expression is obtained for the position of the meniscus frontz m as a function of timet, which in dimensional form isz Here,k is the mobility parameter and has units of (length) 2 /time, λ = H/W is the channel aspect ratio and ζ o (λ) is an aspect-ratio function defined as Detailed derivation of ζ o (λ) can be found in the work of Ouali et al. (2013). Equation (3.1) will be referred to as the MLW model and has been used in several other studies (Baret et al. 2007;Yang et al. 2011;Ouali et al. 2013;Sowers et al. 2016;Kolliopoulos et al. 2019).

Lubrication-theory-based model
We develop a model describing capillary flow in open rectangular channels that accounts for the non-flat shape of the upper liquid-air interface, which results in capillary pressure gradients that drive flow. This is a more complex model than the MLW model (3.1), which assumes a flat upper liquid-air interface.

Model geometry
We begin by considering flow in an open rectangular channel as depicted in figure 4, motivated by the experiments in § 2.3. Recall that λ c is the aspect ratio at which the circular upper meniscus contacts the bottom of the rectangular channel while being attached to the top of the channel sidewalls with a contact angle of θ 0 .
For λ ≥ λ c (figure 4a) the free-surface morphology is divided into three regimes along thez axis as discussed in § 2.3: a meniscus-deformation regime [0,z d (t)], a meniscus-recession regime [z d (t),z m (t)] and a corner-flow regime [z m (t),z t (t)]. In the meniscus-deformation regime, the liquid is pinned to the top of the channel wall (ã = H). The channel inlet is assumed to be fully filled (θ(0,t) = π/2) and the upper meniscus  curvature increases down the channel length (θ(z,t) decreases) until the contact angle θ(z d ,t) = θ 0 , wherez d is the end of meniscus-deformation regime.
We then transition to the meniscus-recession regime, where the contact angle θ = θ 0 and the liquid height starts to recede down the channel sidewalls (ã(z,t) decreases) until the upper meniscus contacts the channel bottom. From (2.1), hereã(z m ,t) = Wλ c , wherẽ z m is the meniscus position. This results in a morphology transition where the flow splits into the channel corners, leading to the corner-flow regime. In the corner-flow regime, θ = θ 0 at the channel bottom and the sidewall, and the liquid height on the sidewallã(z,t) decreases fromã(z m ,t) = Wλ c toã(z t ,t) = 0, wherez t is the finger tip position.
For λ < λ c (figure 4b) the free-surface morphology is also divided into three regimes: as discussed in § 2.3. In the meniscus-deformation regime the liquid is pinned to the top of the channel wall (ã = H). The channel inlet is assumed to be fully filled (θ(0,t) = π/2) and the upper meniscus curvature increases down the channel length (θ(z,t) decreases) until the contact angle θ(z m ,t) = θ C , where θ C is the contact angle at the channel sidewall when the upper meniscus touches the channel bottom.
After the upper meniscus contacts the channel bottom, it splits into the channel corners, leading to the corner-transition regime. In this regime the liquid remains pinned to the top of the channel wall (ã = H), and we assume that the contact angle at the channel bottom reaches θ 0 instantaneously. To conserve mass, the contact angle at the sidewall must change from θ C to θ(z m ,t) = θ T , where θ T is defined in § 3.2.4 (note that θ T = θ C if θ 0 = 0). The upper meniscus curvature increases down the channel length (θ(z,t) at the sidewall decreases) until θ(z c ,t) = θ 0 , wherez c is the finger depinning position (the position at which the liquid depins from the top of the channel wall). Once θ = θ 0 on the channel bottom and sidewall, the morphology transitions to the corner-flow regime where the liquid height on the channel sidewallã(z,t) decreases fromã(z c ,t) = H tõ a(z t ,t) = 0.
In the following sections we develop a mathematical model for capillary flow considering both λ ≥ λ c (figure 4a) and λ < λ c (figure 4b) and accounting for the complex upper liquid-air interface morphology.

Governing equations
We consider mass and momentum conservation of an incompressible Newtonian liquid with constant density, given by∇ is the velocity field in Cartesian coordinates,p is the liquid pressure andg = (g x ,g y ,g z ) is the gravitational acceleration. The no-slip and no-penetration conditions are applied along the solid walls as The boundary conditions for the normal and tangential stresses at the liquid-air interfacẽ h(x,z,t) are given by Here,T = −pI + μ[∇ũ + (∇ũ) T ] is the stress tensor, I is the identity tensor,∇ s = ∇ − n(n ·∇) is the surface gradient operator, n is the unit outward normal vector and t 1 and t 2 are the two tangent vectors at the interface in the transverse and axial directions, respectively (expressions for these vectors can be found in the supplementary material available at https://doi.org/10.1017/jfm.2020.986). Equations (3.3a) and (3.3b) are rendered dimensionless using the following scalings: Additionally, the gravitational acceleration vector is scaled as (g x ,g y ,g z ) = (gg x , gg y , gg z ), where g is the magnitude of the gravitational acceleration. The dimensionless parameters that arise are the Reynolds number Re = ρUH/μ (ratio of inertial to viscous forces), the capillary number Ca = μU/ σ (ratio of viscous to surface tension forces) and the Bond number Bo = ρgH 2 /σ (ratio of gravitational to surface tension forces).
In the limits where 2 1, Re 1 and Bo/Ca , the governing equations reduce to The boundary conditions for the normal (3.5a), transverse tangential (3.5b) and axial tangential (3.5c) stresses at the free surface reduce to The normal stress balance in (3.5a) has as a special case the Young-Laplace equation p = −Ca −1 κ, where κ accounts for both transverse and axial curvature contributions. However, in the limit 2 1, axial curvature contributions are negligible and only the leading-order transverse curvature contributions are accounted for in (3.8a). Based on (3.7b) the O(1) term in p is only dependent on z and t, and thus the leading-order curvature term (term in brackets on the far right of (3.8a)) is actually independent of x and must only depend on z and t. The derivation of (3.7) and (3.8) can also be seen in Yang & Homsy (2006) and White & Troian (2019), who considered V-shaped channel cross-sections.
Up to this point no assumption has been made regarding the channel cross-sectional geometry. Here, we consider two geometries for the channel cross-section: (a) rectangular (figure 5a) and (b) V-shaped (figure 5b). Using these two geometries we can describe all the liquid cross-sections in figures 4(a) and 4(b) in terms of the liquid height on the solid wall a(z, t) and the contact angle θ(z, t). The meniscus-deformation (a = 1) and meniscus-recession (θ = θ 0 ) regimes are described using the rectangular cross-section, while the corner-transition (a = 1) and corner-flow (θ = θ 0 ) regimes are described using the V-shaped cross-section.
Each cross-sectional geometry requires three additional boundary conditions to obtain expressions for p(z, t) and h(x, z, t): the contact-line location on the solid wall, a symmetry condition and the definition of the contact angle θ. Expressions for these boundary conditions can be found in the supplementary material. We obtain expressions for p (z, t) and h(x, z, t) as a function of a(z, t) and θ(z, t) for each regime in figure 4 by integrating (3.8a) twice with respect to x and imposing the boundary conditions. The resulting O(1) expressions are Figure 5. Cross-sectional schematics of (a) meniscus-deformation (a = 1) and meniscus-recession (θ = θ 0 ) regimes, and (b) corner-transition (a = 1) and corner-flow (θ = θ 0 ) regimes.

Diffusion equations
(3.10) The dimensionless flux is defined as wherew i is a rescaled cross-sectional-averaged dimensionless velocity. Here, i is equal to D, T or C for the meniscus-deformation, corner-transition and corner-flow regimes, respectively. (As discussed below, the meniscus-recession regime is neglected.) Details of the calculation ofw i are discussed in the supplementary material. It is evident from (3.9a)-(3.9d) that the dimensionless streamwise flux q in (3.11) is a function of either the dimensionless liquid height a(z, t) on the sidewall or the liquid contact angle θ(z, t) depending on the regime. Rather than considering a(z, t) and θ(z, t) separately, we introduce the liquid saturation s =Ã(a, θ)/HW = λA (ratio of channel cross-sectional area filled with liquid to total channel cross-sectional area). For each regime the liquid saturation s is given by where the geometric functionsB andÂ can be found, respectively, in (A1) and (A2) in appendix A. Equations (3.12a), (3.12b) and (3.12d) are equivalent to expressions reported by Nilson et al. (2006). Since the pressure p in (3.9b) is constant in the meniscus-recession regime, the flux q = 0 for this regime based on (3.11). This is because the transverse curvature gradients are zero and the only contribution to q is from the O( 2 ) axial curvature gradients, for which we did not account. Nilson et al. (2006) estimated that the meniscus-recession regime size δ ≈ L( 2 λ/2) 1/3 . For the microchannel dimensions considered in our study δ ≈ 180-320 μm, which is negligible considering the channel length is 30 mm. This estimate for δ agrees with the observations in figure 3(b) where the meniscus-recession regime size is less than 50 μm. Therefore, effects of the meniscus-recession regime will be neglected (i.e.z d =z m in figure 4a) and the regime transition from meniscus deformation to corner flow (for λ/λ c > 1) will be treated as a saturation jump. By using (3.12a)-(3.12d) in (3.10) we obtain the following system of nonlinear partial differential equations governing the liquid saturation: where D D =w D 2λ 2 sin θ cos 2 θ 1 − tan θ arcsin(cos θ) , (3.14a) The quantities D D , D T and D C can be thought of as dimensionless diffusion coefficients describing the interface growth.
Recall that λ c (see (2.1)) is the aspect ratio at which the circular upper meniscus contacts the bottom of the rectangular channel while being attached to the top of the channel sidewalls with a contact angle of θ 0 . When λ ≥ λ c (figure 4a), the bounds of the meniscus-deformation and corner-flow regimes are (0, z m ) and (z m , z t ), respectively, where z m is the meniscus position and z t is the finger tip position. When λ < λ c (figure 4b), the bounds of the meniscus-deformation, corner-transition and corner-flow regimes are (0, z m ), (z m , z c ) and (z c , z t ), respectively, where z c is the finger depinning position.

Similarity transformation
where η 0 = z t / √ t is the rescaled finger tip position and δ 0 η 0 = z m / √ t is the rescaled meniscus position. The channel cross-section at the inlet is assumed to be fully filled and θ = π/2. At the end of the meniscus-deformation regime θ = θ 0 , which is used in (3.12a) to calculate s D .
At the beginning of the corner-flow regime a = λ c /λ, which is used in (3.12d) to determine s C . (Recall from § 3.2.1 that the corner-flow regime begins whenã = Wλ c , which in dimensionless form is a = λ c /λ.) Finally at the finger tip, the liquid height goes to zero. Note that s D = s C only for λ = λ c (meniscus contacts the channel bottom at end of the meniscus-deformation regime). Equation (2.1) is used to determine λ c , which depends only on θ 0 . For all λ, it is assumed in the corner-flow regime that the contact angle on the channel sidewall and bottom is always θ 0 , and thus independent of speed. This is the simplest assumption and allows us to focus on the influence of other problem parameters.
Two additional conditions are required to determine η 0 and δ 0 , which specify the bounds of each regime. The first condition is the flux matching condition given by where the second term in each bracket accounts for the potential discontinuity in s due to transitioning from the meniscus-deformation regime to the corner-flow regime. A derivation of (3.15d) can be seen in the appendix B. The second condition is that the flux approaches zero at the finger tip (i.e. D C s 1/2 ds/dη → 0, as η → η 0 ). Following Romero & Yost (1996) and using (3.15b), it can be shown that to satisfy this condition, the following must be true: For λ < λ c shown in figure 4(b), the self-similar governing equations are subject to where δ 1 η 0 = z c / √ t is the rescaled finger depinning position ( § 3.2.1). The channel cross-section at the inlet is assumed to be fully filled and θ = π/2. At the end of the meniscus-deformation regime θ = θ C (critical angle at which the upper meniscus touches the channel bottom, calculated from λ = (1 − sin θ C )/2 cos θ C ), which is used in (3.12a) to calculate s D , and the contact angle at the channel bottom is θ = 0.
At the transition from the meniscus-recession regime to the corner-transition regime, the liquid remains pinned to the top of the channel sidewall and the upper meniscus contacts the channel bottom with the flow splitting into the channel corners. At the beginning of the corner-transition regime we assume the liquid instantaneously attains θ 0 at the channel bottom and θ C → θ T at the channel sidewall. To conserve mass, we equate the amount of liquid in the channel cross-section on each side of this transition. This specifies θ T , which is calculated (via Newton's method) by setting (3.12c) equal to s D . If the calculated θ T ≤ π/4, then s T = s D . If θ T > π/4, then D T < 0 which makes the problem ill-posed (Romero & Yost 1996). In this case, we set θ T = π/4, leading to a saturation jump. Equation (3.12c) is then used to determine s T based on θ T . Note that in the corner-transition regime the contact line at the channel sidewall is assumed to be pinned while the contact line at the channel bottom is allowed to move with constant contact angle θ 0 . At the end of the corner-transition regime and the beginning of the corner-flow regime a = 1, which is used in (3.12d) to determine s C . (Recall from § 3.2.1 that the corner-flow regime begins wheñ a = H, which in dimensionless form is a = 1.) Finally at the finger tip, the liquid height goes to zero.

A32-15
Three additional conditions are required to determine η 0 , δ 0 and δ 1 , which are where (3.16a) and (3.16b) are flux matching conditions (see appendix B) and (3.16c) is the condition setting the flux to zero at the finger tip. For λ ≥ λ c the system of governing equations is (3.15), whereas for λ < λ c the system consists of (3.16). What is required to solve these systems are the cross-sectional-averaged dimensionless velocitiesw D (s),w T (s) andw C (s), which influence the values of D D , D T and D C through (3.14). The cross-sectional-averaged dimensionless velocities are calculated for a given cross-section by solving (3.7c) subject to no-slip and no-penetration conditions along the solid walls and no-stress condition (3.8c) at the liquid-air interface (see supplementary material for further details).

Numerical methods
Velocity fields (see supplementary material) are numerically solved for with a Galerkin finite-element method using quadratic basis functions. To validate our computations, our results forw D (s),w T (s) andw C are compared to results from prior studies. Results for w D (s) andw T (s) are in agreement with results of Tchikanda et al. (2004) and Weislogel & Nardin (2005), respectively. Results forw C agree with results of Ayyaswamy, Catton & Edwards (1974), Ransohoff & Radke (1988) and Yang & Homsy (2006). Note that these prior studies do not consider capillary flow in open rectangular channels over the range of contact angles θ 0 and aspect ratios λ examined in the present work.
Results forw D (s) andw T (s) from the finite-element method simulations are fitted using Chebyshev polynomials of the first kind using the least-squares method. These Chebyshev polynomials are then used in the system of equations (3.15) (λ ≥ λ c ) and (3.16) (λ < λ c ) to evaluate D D and D T . Sincew C does not depend on s, an exact expression for D C can be obtained via (3.14c). Both nonlinear systems of equations (3.15) and (3.16) are discretized using a second-order centred finite-difference method and solved using the fsolve solver in MATLAB.

Results and discussion
Similarity solutions for the liquid saturation profiles s(η) and their dependence on the channel aspect ratio λ and equilibrium contact angle θ 0 are presented first ( § 4.1). Using these similarity solutions, three-dimensional liquid height profiles are obtained to highlight the complex free-surface morphology similar to that seen in figure 3. Model predictions for the evolution of the meniscus positionz m (t) from the lubrication-theory-based and MLW models are then compared to experimental observations ( § 4.2). Finally, lubrication-theory-based model predictions of the finger length evolution l f (t) = z t (t) − z m (t) are compared to experimental results ( § 4.3).

Saturation profiles
Computed similarity solutions of the liquid saturation s(η) for different aspect ratios λ and an equilibrium contact angle θ 0 = 10 • are shown in figure 6(a). Solutions for λ > λ c and λ < λ c are obtained solving the system of equations (3.15) and (3.16), respectively. These similarity solutions are valid for intermediate times, when channel entrance and end effects can be neglected.
In figure 6(a) when λ > λ c (here λ c = 0.42), s(0) = 1 corresponds to a fully filled channel cross-section. Moving down the length of the channel, s decreases monotonically and at the meniscus position (filled symbol) the flow transitions from the meniscus-deformation regime to the corner-flow regime (see figure 4a). A jump in s (dashed lines) is observed because we neglected the meniscus-recession regime as discussed in § 3.2.3. In the corner-flow regime s continues to decrease until s(η 0 ) = 0 at the finger tip.
From figure 6(a), the s(η) profiles have a non-monotonic dependence on λ, suggesting that there is an optimal λ for capillary flow. The effect of the equilibrium contact angle θ 0 on s(η) for an aspect ratio λ = 0.75 is shown in figure 6(b). Decreasing θ 0 results in more capillary filling. Although in figure 6(b) we consider λ > λ c , the same trend is observed for λ < λ c .
The non-monotonic effect of λ on capillary flow becomes more clear in figure 7. Here, the rescaled finger tip position η 0 , meniscus position δ 0 η 0 and finger depinning position δ 1 η 0 (defined in § 3.2.4) are presented as a function of λ for different θ 0 . The shaded areas between the curves represent the sizes of the meniscus-deformation (I), corner-flow (III) and corner-transition (IV) regimes (seen in figures 4a and 4b), which depend on λ and θ 0 .
We first consider results in figure 7(a), where θ 0 = 10 • . When λ λ c the flow is dominated by the meniscus-deformation (I) regime. With decreasing λ, the size of the corner-flow (III) regime monotonically increases. However, the size of the meniscus-deformation (I) regime increases and then decreases, with decreasing λ. When λ drops below λ c , the corner-transition (IV) regime appears. As λ is further decreased, the sizes of the corner-flow (III) and corner-transition (IV) regimes increase, while the size of the meniscus-deformation (I) regime decreases. These trends are observed for the other values of θ 0 considered in figure 7(b-d).  Similarity solutions for s(η) are used to construct three-dimensional free-surface profiles. These solutions for s(η) are used in (3.12) to determine θ(η) (via Newton's method) and a(η) for each regime. The three-dimensional free-surface profiles h are determined using (3.9). Since the h expressions in (3.9c) and (3.9d) for the corner-transition and corner-flow regimes correspond to figure 5(b), they require rotation by angles β = arctan(cos θ/ cos θ 0 ) and β = π/4, respectively, to match the channel orientation seen in figure 4. By solving the system of ordinary differential equations (3.15) and (3.16) we can construct three-dimensional free-surface profiles for any timẽ t. Dimensional free-surface profilesh for λ > λ c and λ < λ c after 100 s of flow of NOA74 are depicted in figure 8.
Qualitative agreement is observed between the three-dimensional free-surface profiles in figure 8 and the profilometry measurements in figures 3(a) and 3(c). The channel is fully filled at the channel inlet and the upper meniscus bows down as we move down the length of the channel until it contacts the channel bottom and splits into the channel corners. A quantitative comparison between theory and experiment is made in the following section.

Comparison with experiments
We compare predictions of the meniscus-position time evolutionz m (t) from the lubrication-theory-based and MLW models to experimental observations. This comparison is made in figure 9 for the test liquids detailed in table 1. Solid and dashed lines represent the lubrication-theory-based and MLW model predictions, respectively, while experimental observations are shown as symbols. Each panel in figure 9 includes experiments for λ > λ c and λ < λ c , except for the case of propylene glycol (figure 9d) where only experiments for λ > λ c were conducted. When λ > λ c , the lubrication-theory-based and MLW models are in good agreement with experiments for NOA74 (figure 9a) and silicone oil (figure 9b). For mineral oil (figure 9c) and propylene glycol (figure 9d), the lubrication-theory-based model agrees well with the experiments but the MLW model underpredictsz m (t).
For NOA74 (figure 9a) and silicone oil (figure 9b), the lubrication-theory-based model agrees with experiments for λ < λ c , but the MLW model overpredictsz m (t). For mineral oil (figure 9c), both models overpredictz m (t) for λ < λ c . For propylene glycol (figure 9d), at larger λ the lubrication-theory-based model predictions agree better with experiments than predictions of the MLW model. However, as λ → λ c the MLW model predictions agree with experiments, whereas the lubrication-theory-based model overpredictsz m (t).
To further compare theory and experiment we consider the dimensionless mobility parameter k = z 2 m /t, which can be thought of as a diffusion coefficient driving the growth of the meniscus position. In the lubrication-theory-based model, wherek is defined in (3.1), U is the characteristic velocity (see § 3.2) and L is the channel length. Experimentally, k is determined by fitting the function z m = (kt) n to λ = 1.36 experiments similar to those shown in figure 9 using nonlinear regression. A comparison of theoretically predicted and experimentally determined k values as a function of λ and θ 0 is shown in figure 10. When λ λ c , the lubrication-theory-based and MLW model predictions are indistinguishable regardless of θ 0 . In this limit the finger contribution to the flow becomes negligible and the flow asymptotically approaches that between two parallel plates. Note that when λ = H/W λ c , the effects from the channel bottom become negligible compared with the effects from the sidewalls. As λ → 1 the two model predictions begin to deviate because the viscous resistance in the fingers becomes significant, and fingers are not accounted for in the MLW model. However, both model predictions for k are in reasonable agreement with experiments. As λ → λ c the lubrication-theory-based model agrees better with experiments, compared to the MLW model (except in figure 10d). The disagreement between the MLW model and experiments becomes more evident as λ decreases because the free-surface height profiles (figures 3b and 3d) become more non-uniform. Therefore, the assumption of a flat free-surface height profile used in the MLW model becomes less valid with decreasing λ.
However, the lubrication-theory-based model does not always perform better than the MLW model (figure 10d). As θ 0 → π/4 the transverse free-surface curvature in the corner-transition and corner-flow regimes vanishes and the flow in these regimes is then driven by axial curvature gradients (Yang & Homsy 2006). These axial curvature gradients are neglected in the lubrication-theory-based model (see (3.8a)) because they are O( 2 ) terms. When the axial curvature contributions are included in the pressure, then where C T (s) and C A (s, ∂ z s, ∂ 2 z s) are the transverse and axial curvature contributions, respectively. Including the axial curvature gradients in the model results in fourth-order spatial derivatives of s, for which we were unable to find similarity solutions. An additional reason for the deviation between the lubrication-theory-based model predictions and the experimental observations as θ 0 → π/4 is that the lubrication approximation itself is expected to become less accurate as the contact angle increases.
For predicting z m (t), the choice of using the lubrication-theory-based or MLW model depends on λ and θ 0 . For large λ both models give similar results, so the MLW model is preferred because of its simplicity, but for small λ the lubrication-theory-based model is more accurate. Additionally, the lubrication-theory-based model performs better for smaller θ 0 compared to the MLW model. However, it neglects key physical contributions as θ 0 → π/4 which results in poorer agreement. Nevertheless, when additional information is desired, such as the free-surface morphology (figure 8) or finger dynamics ( § 4.3), then the lubrication-theory-based model must be chosen.

Finger dynamics
A key advantage of the lubrication-theory-based model compared with the MLW model is that it describes finger dynamics. The dimensionless finger length l f (t) is defined as the distance between the dimensionless finger tip position z t (t) and the dimensionless meniscus position z m (t) which can be seen in figure 4. Predictions of the dimensionless finger length time evolution l f (t) by the lubrication-theory-based model are compared to experimental observations in figure 11. Model predictions for l f were computed using where η 0 = z t / √ t and δ 0 η 0 = z m / √ t; values for η 0 and δ 0 η 0 are shown in figure 7. The effects of λ and θ 0 on l f (t) are shown in figures 11(a) and 11(b). Solid lines and symbols represent lubrication-theory-based model predictions and experimental observations, respectively. It is important to note the shorter dimensionless times compared to the meniscus-tracking experiments (figure 9) were caused by the finger tip moving out of the field of view. In figure 11(a), experimental observations are compared to model predictions for mineral oil (λ c = 0.28), where increasing channel aspect ratio λ results in a decrease in l f (t). Good agreement between theory and experiments is observed, with model predictions being within the range of experimental observations for all but one trial (i.e. λ = 0.3 in figure 11a). A possible reason for this discrepancy is the error in identifying the position of the finger tip where the liquid height goes to zero, which is challenging during the meniscus-tracking experiments where the stage is moving (figure 2a).
In figure 11(b), experimental observations for mineral oil (θ 0 = 32 • , λ c = 0.28) and silicone oil (θ 0 = 18 • , λ c = 0.36) are compared to model predictions for λ = 0.23 (λ < λ c ). Increasing θ 0 leads to a decrease in l f (t). Model predictions are in good agreement with experimental observations in this case as well. Hence, the lubrication-theory-based model accurately captures effects of λ and θ 0 on l f (t). Additionally, the good agreement between theory and experiment in figure 11 emphasizes the importance of accounting for the corner-transition regime (figure 4b) when describing the finger dynamics in open rectangular channels when λ < λ c .
We make two additional comments about the fingers. First, although curvature of the channel corners may affect finger length (Chen et al. 2006;Gerlach et al. 2020), the good agreement observed here suggests that the corner curvature does not play a significant role in our experiments. Second, although the lubrication-theory-based model overpredicts the meniscus position for mineral oil (figure 9d), it accurately predicts the finger length (figure 11) which is the difference between the positions of the meniscus and finger tip (4.3). This indicates that while the lubrication-theory-based model may not predict absolute positions accurately, it can predict relative positions accurately. A deeper understanding of the reasons for this likely requires a comparison of the lubrication-theory-based model and direct numerical simulations, which is beyond the scope of the present work.

Conclusions
In this work we combine theory and experiment to examine capillary-flow dynamics in open rectangular microchannels. For open microchannels, the free surface greatly influences the capillary-flow dynamics. We visualize the free-surface morphology and its dependence on the channel aspect ratio λ using SEM and profilometry. The SEM images suggest a qualitative difference in the free-surface morphology at λ = λ c and highlight the significance of a corner-transition regime when λ < λ c , which was not accounted for in previous studies of capillary flow in open rectangular channels.
Effects of the free-surface morphology on capillary-flow dynamics were examined using two theoretical models. The first model is a MLW model, which assumes a flat free surface and has been extensively used in prior studies. The second model is a self-similar lubrication-theory-based model, which was developed to account for the complex free-surface morphology observed in the experiments.
Predictions of the lubrication-theory-based and MLW models were compared to complementary flow visualization experiments over a wide range of channel aspect ratios λ and equilibrium contact angles θ 0 . For large λ, predictions from the two models are indistinguishable, since free-surface morphology effects are negligible. However, for smaller λ the lubrication-theory-based model is in better agreement with experiments since the influence of fingers is significant. Additionally, the lubrication-theory-based model agrees better with experiments for smaller θ 0 , although as θ 0 → π/4 the agreement worsens because important axial curvature contributions are neglected. Finally, the lubrication-theory-based model predictions accurately capture the finger length time evolution seen in experiments over a range of λ and θ 0 , further highlighting the importance of accounting for the corner-transition regime in the model formulation when λ < λ c .
Our lubrication-theory-based model and flow visualization experiments reveal the limitations of the widely used MLW model. In addition, our results significantly advance physical understanding of capillary-flow dynamics in open rectangular microchannels, which play a key role in a number of technological applications ( § 1) but are less well studied than V-shaped channels. Finally, the systematic development of our model allows for it to be readily extended to incorporate other important phenomena such as gravitational effects, solute transport and evaporation.