2.1. Model geometry

An open rectangular channel with width $\hat {W}$, height $\hat {H}$ and length $\hat {L}$, connected to a reservoir of radius $\hat {R}$ is depicted in figure 1. The amount of liquid in the reservoir is described using the contact angle on the reservoir sidewall $\theta _R$. The contact line is assumed to be pinned to the top of the reservoir sidewall. As described in Kolliopoulos et al. (Reference Kolliopoulos, Jochem, Johnson, Suszynski, Francis and Kumar2021), the free surface undergoes a transition from that observed in figure 1(a) ($\lambda \ge \lambda _c$) to that in figure 1(b) ($\lambda <\lambda _c$). Here, $\lambda =\hat {H}/\hat {W}$ is the channel aspect ratio and $\lambda _c=(1-\sin \theta _0)/2\cos \theta _0$ 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 an equilibrium contact angle $\theta _0$.

Figure 1. Schematic of liquid undergoing capillary flow in an open rectangular channel connected to a circular reservoir for aspect ratios (a) $\lambda \ge \lambda _c$ and (b) $\lambda <\lambda _c$. Here, $\beta =\arctan (\cos \theta /\cos \theta _0)$ and the finger length is $\tilde {z}_T-\tilde {z}_M$.

Each regime is described using the liquid height $\hat {a}$ and the contact angle $\theta$ on the channel sidewall. For $\lambda \ge \lambda _c$ (figure 1a), the free-surface morphology is divided into three regimes along the channel length: meniscus deformation $[0,\hat {z}_D]$, meniscus recession $[\hat {z}_D,\hat {z}_M]$ and corner flow $[\hat {z}_M,\hat {z}_T]$. In the meniscus-deformation regime, the upper meniscus curvature at the channel inlet is matched to the liquid–air interface curvature in the reservoir and the upper meniscus curvature increases down the channel length until $\theta (\hat {z}_D,\hat {t})=\theta _0$, while the liquid remains pinned to the top of the channel sidewall $\hat {a}=\hat {H}$.

In the meniscus-recession regime, the contact angle on the sidewall remains constant at $\theta =\theta _0$ and the liquid height recedes down the sidewall until the upper meniscus contacts the channel bottom $\hat {a}(\hat {z}_M,\hat {t})=\hat {W}\lambda _c$. At this point, the contact angle on the bottom is assumed to reach $\theta _0$ instantaneously. This is the simplest possible assumption, implies a neglect of contact-angle hysteresis, and works well in describing experiments in the absence of evaporation (Kolliopoulos et al. Reference Kolliopoulos, Jochem, Johnson, Suszynski, Francis and Kumar2021). Subsequently, the flow splits into the channel corners, leading to the corner-flow regime. Here, the contact angle on the sidewall and bottom remains constant at $\theta =\theta _0$ and the liquid height further recedes down the sidewall from $\hat {a}(\hat {z}_M,\hat {t})=\hat {W}\lambda _c$ to $\hat {a}(\hat {z}_T,\hat {t})=0$ at the finger tip.

For $\lambda <\lambda _c$ (figure 1b) the free-surface morphology is also divided into three regimes: meniscus deformation $[0,\hat {z}_M]$, corner transition $[\hat {z}_M,\hat {z}_C]$ and corner flow $[\hat {z}_C,\hat {z}_T]$. In the meniscus-deformation regime, the liquid is pinned to the top of the sidewall ($\hat {a}=\hat {H}$). The upper meniscus curvature at the channel inlet is matched to the liquid–air interface curvature in the reservoir and the upper meniscus curvature increases down the channel length until $\theta (\hat {z}_M,\hat {t})=\theta _C$ when the upper meniscus touches the channel bottom.

Upon contact with the channel bottom, the contact angle there is assumed to reach $\theta _0$ instantaneously. The upper meniscus splits into the channel corners, yielding the corner-transition regime. During this stage, the liquid remains pinned to the top of the sidewall ($\hat {a}=\hat {H}$). To conserve mass, the contact angle at the sidewall $\theta (\hat {z}_M,\hat {t})$ must change from $\theta _C$ to another value denoted $\theta _T$. The contact angle on the sidewall decreases down the channel length until $\theta (\hat {z}_C,\hat {t})=\theta _0$, at which point the flow transitions to the corner-flow regime. Here, $\theta =\theta _0$ and the liquid depins from the top of the sidewall and decreases until $\hat {a}(\hat {z}_T,\hat {t})=0$.

In the following sections we develop a mathematical model for capillary flow of a liquid solution with evaporation considering both $\lambda \ge \lambda _c$ (figure 1a) and $\lambda <\lambda _c$ (figure 1b), and account for the finite-size reservoir used in experiments (Lade et al. Reference Lade, Jochem, Macosko and Francis2018).

2.2. Hydrodynamics

Mass and momentum conservation for an incompressible Newtonian liquid are governed by

(2.1a)\begin{align} \boldsymbol{\widehat{\nabla}} \boldsymbol{\cdot} \boldsymbol{\hat{u}}=0 \quad \text{and} \quad \hat{\rho} \left[\frac{\partial \boldsymbol{\hat{u}}}{\partial \hat{t}}+ (\boldsymbol{\hat{u}} \boldsymbol{\cdot} \boldsymbol{\widehat{\nabla}})\boldsymbol{\hat{u}}\right]={-}\boldsymbol{\widehat{\nabla}}\hat{p} +\hat{\mu} \widehat{\nabla}^2 \boldsymbol{\hat{u}} + \hat{\rho} \boldsymbol{\hat{g}}, \end{align}

where $\boldsymbol {\hat {u}}=(\hat {u},\hat {v},\hat {w})$ is the velocity in Cartesian coordinates $(\hat {x},\hat {y},\hat {z})$, $\hat {p}$ is the liquid pressure and $\boldsymbol {\hat {g}}=(\hat {g}_x,\hat {g}_y,\hat {g}_z)$ is the gravitational acceleration. Along the channel walls, the no-slip and no-penetration conditions require $\boldsymbol {\hat {u}}=0$. The boundary conditions for the jump in normal, transverse tangential, and axial tangential stresses across the liquid–air interface $\hat {h}(\hat {x},\hat {z},\hat {t})$ are

(2.1b)\begin{align} [\![\boldsymbol{n}\boldsymbol{\cdot} \hat{\boldsymbol{T}}\boldsymbol{\cdot}\boldsymbol{n}]\!]=\hat{\sigma} (\boldsymbol{\widehat{\nabla}}_s \boldsymbol{\cdot} \boldsymbol{n}), \quad [\![\boldsymbol{t}_1\boldsymbol{\cdot} \hat{\boldsymbol{T}}\boldsymbol{\cdot}\boldsymbol{n}]\!]=\boldsymbol{t}_1 \boldsymbol{\cdot} \boldsymbol{\widehat{\nabla}}_s \hat{\sigma}, \quad \text{and} \quad [\![\boldsymbol{t}_2\boldsymbol{\cdot} \hat{\boldsymbol{T}}\boldsymbol{\cdot}\boldsymbol{n}]\!]=\boldsymbol{t}_2 \boldsymbol{\cdot} \boldsymbol{\widehat{\nabla}}_s \hat{\sigma}. \end{align}

Here, $\hat {\boldsymbol {T}}=-\hat {p}\boldsymbol {I}+\hat {\mu }[\boldsymbol {\widehat {\nabla }}\boldsymbol {\hat {u}}+(\boldsymbol {\widehat {\nabla }}\boldsymbol {\hat {u}})^{\rm T}]$ is the stress tensor, $\boldsymbol {I}$ is the identity tensor, $\boldsymbol {\widehat {\nabla }}_s=\boldsymbol {\widehat {\nabla }}-\boldsymbol {n}(\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {\widehat {\nabla }})$ is the surface-gradient operator, $\boldsymbol {n}$ is the unit outward normal vector and $\boldsymbol {t}_1$, $\boldsymbol {t}_2$ are the two unit tangent vectors to the interface in the transverse and axial directions, respectively. These vectors are given by

(2.2a–c)\begin{equation} \boldsymbol{n}=\frac{(-\hat{h}_{\hat{x}},1,-\hat{h}_{\hat{z}})}{\sqrt{1+\hat{h}_{\hat{x}}^2+\hat{h}_{\hat{z}}^2}},\quad \boldsymbol{t}_1=\frac{(1+\hat{h}_{\hat{z}}^2,\hat{h}_{\hat{x}},-\hat{h}_{\hat{x}}\hat{h}_{\hat{z}})}{\sqrt{(1+\hat{h}_{\hat{z}}^2)^2 +\hat{h}_{\hat{x}}^2 +\hat{h}_{\hat{x}}^2\hat{h}_{\hat{z}}^2}},\quad \boldsymbol{t}_2=\frac{(0,\hat{h}_{\hat{z}},1)}{\sqrt{1+\hat{h}_{\hat{z}}^2}}. \end{equation}

Equations (2.1a) are rendered dimensionless using the following scalings:

(2.3)\begin{gather} \left. \begin{gathered} (\hat{x},\hat{y},\hat{z})=\hat{L}(\epsilon x,\epsilon y, z), \enspace (\hat{u},\hat{v},\hat{w})=\hat{U}(\epsilon u,\epsilon v,w), \enspace (\hat{g}_x,\hat{g}_y,\hat{g}_z)=\hat{g}(g_x,g_y,g_z),\\ \hat{t}=(\hat{L}/\hat{U})t, \enspace \epsilon=\hat{H}/\hat{L}, \quad \hat{U}=\epsilon\hat{\sigma}_0/\hat{\mu}_0, \enspace \hat{p}=(\hat{\mu}_0 \hat{U}/\epsilon \hat{H})p, \enspace \hat{\mu}=\hat{\mu}_0M, \enspace \hat{\sigma}=\hat{\sigma}_0\varSigma,\\ \end{gathered} \right\} \end{gather}

where $\epsilon$ is the ‘slenderness’ parameter, $\hat {U}$ is the viscocapillary velocity, $\hat {\mu }_0$ is the solvent viscosity, $\hat {\sigma }_0$ is the solvent surface tension and $M$ and $\varSigma$ are the dimensionless solution viscosity and surface tension, respectively, defined in § 2.6. The dimensionless numbers that arise are the Reynolds number $Re=\hat {\rho } \hat {U}\hat {H}/\hat {\mu }_0$ (ratio of inertial to viscous forces), the capillary number $Ca=\hat {\mu }_0\hat {U}/\epsilon \hat {\sigma }_0$ (ratio of viscous to surface-tension forces) and the Bond number $Bo=\hat {\rho } \hat {g}\hat {H}^2/\hat {\sigma }_0$ (ratio of gravitational to surface-tension forces). Note that $Ca=1$ based on our choice of $\hat {U}$.

In the limits where $\epsilon ^2\ll 1$, $\epsilon Re\ll 1$ and $Bo/Ca\ll \epsilon$, (2.1a) reduces to

(2.4a)\begin{gather} u_x + v_y + w_z =0, \end{gather}

(2.4b)\begin{gather} p_x = p_y=0, \end{gather}

(2.4c)\begin{gather} p_z= M(w_{xx} +w_{yy}). \end{gather}

The boundary conditions at the free surface in (2.1b) reduce to

(2.5a)\begin{gather} p-p_{V,T}={-}\frac{\varSigma}{Ca} \left[\frac{h_{x}}{(1+h_{x}^2)^{1/2}}\right]_x ={-}\frac{\varSigma}{Ca} \kappa, \end{gather}

(2.5b)\begin{gather} (1-h_x^2)(v_x+u_y) +2h_x({-}u_x+v_y)\,{-}\,(1-h_x^2)h_z w_x \,{-}\, 2h_xh_z w_y\nonumber\\ \hspace{-8pc}={-}h_xh_z(1+h_x^2)^{1/2}\frac{\varSigma_z}{MCa}, \end{gather}

(2.5c)\begin{gather} w_y-h_x w_x = \frac{\varSigma_z}{MCa}, \end{gather}

where $p_{V,T}$ is the total pressure in the vapour phase and is assumed constant. Since the leading-order pressure term is independent of $(x,y)$ according to (2.4b), the leading-order curvature term $\kappa$ only depends on $(z,t)$ according to (2.5a).

We note two things about these reduced equations. First, the interface shape at a given $z$-position is determined by capillary statics (Yang & Homsy Reference Yang and Homsy2006), as can be seen from (2.5a). As a consequence, contact angles can be imposed as boundary conditions without having to specify a slip law. The interface shape slowly varies in the $z$-direction via an evolution equation to be presented in § 2.7. This approach is expected to hold for $\hat {\mu }_0 \hat {U}/\hat {\sigma }_0 \ll 1$. Second, because of the problem we are considering and our choice of variables, $h_x$ is ${O}(1)$, so the $h_x^2$ terms are retained.

As shown by Kolliopoulos et al. (Reference Kolliopoulos, Jochem, Johnson, Suszynski, Francis and Kumar2021), using a condition for the contact-line location on the solid wall, a symmetry condition, and the definition of the contact angle on the sidewall, expressions for $p(z,t)$ and $h(x,z,t)$ are obtained by integrating (2.5a) twice with respect to $x$, leading to

(2.6a)\begin{gather} \text{meniscus deformation} \quad \begin{cases} p={-}2\varSigma\lambda\cos\theta(z,t) + p_{V,T}, \\ h=1+\dfrac{\tan\theta(z,t)}{2\lambda}-\left[\dfrac{1}{4\lambda^2\cos^2\theta(z,t)}-x^2\right]^{1/2} , \\ A = \dfrac{1}{4\lambda^2}\left[4\lambda-\dfrac{{\rm \pi}/2-\theta(z,t)}{\cos^2\theta(z,t)}+\tan\theta(z,t)\right], \end{cases} \end{gather}

(2.6b)\begin{gather} \text{meniscus recession} \quad \begin{cases} p={-}2\varSigma\lambda\cos\theta_0 + p_{V,T}, \\ h=a(z,t)+\dfrac{\tan\theta_0}{2\lambda}-\left[\dfrac{1}{4\lambda^2\cos^2\theta_0}-x^2\right]^{1/2}, \\ A = \dfrac{1}{4\lambda^2}\left[4\lambda a(z,t)-\dfrac{{\rm \pi}/2-\theta_0}{\cos^2\theta_0}+\tan\theta_0\right], \end{cases} \end{gather}

(2.6c)\begin{gather} \text{corner transition} \quad \begin{cases} p={-}\varSigma[\cos\theta_0-\sin\theta(z,t)] + p_{V,T}, \\ h=\dfrac{\cos\theta\cos\beta}{\cos(\theta+\beta)}-\left[\left(\dfrac{\sin\beta}{\cos(\theta+\beta)}\right)^2-\left(\dfrac{\cos\theta\sin\beta}{\cos(\theta+\beta)}-x-\dfrac{1}{2\lambda}\right)^2\right]^{1/2}, \\ A = \dfrac{B(\theta,\theta_0)}{(\cos\theta_0-\sin\theta)^2}, \end{cases} \end{gather}

(2.6d)\begin{gather} \text{corner flow} \quad \begin{cases} p={-}\dfrac{\varSigma[\cos\theta_0-\sin\theta_0]}{a(z,t)} + p_{V,T}, \\ h=\dfrac{a\cos\theta_0\cos\frac{\rm \pi}{4}}{\cos(\theta_0+\frac{\rm \pi}{4})}-\left[\left(\dfrac{a\sin\frac{\rm \pi}{4}}{\cos(\theta_0+\frac{\rm \pi}{4})}\right)^2-\left(\dfrac{a\cos\theta_0\sin\frac{\rm \pi}{4}}{\cos(\theta_0+\frac{\rm \pi}{4})}-x-\dfrac{1}{2\lambda}\right)^2\right]^{1/2}, \\ A = \dfrac{a^2B(\theta_0,\theta_0)}{(\cos\theta_0-\sin\theta_0)^2}, \end{cases} \end{gather}

where $\beta =\arctan (\cos \theta /\cos \theta _0)$ and $A=2\int _{x_1}^{x_2} h\,\mathrm {d}\kern0.7pt x$ is the dimensionless liquid cross-sectional area. Note that $\theta$ and $a$ in (2.6) are dependent on $(z,t)$. The integration bounds for each regime are

(2.7)\begin{equation} \left. \begin{aligned} & \text{meniscus deformation:} \quad x_1 =0 \quad \text{and} \quad x_2=\dfrac{1}{2\lambda},\\ & \text{meniscus recession:} \quad x_1 =0 \quad \text{and} \quad x_2=\dfrac{1}{2\lambda},\\ & \text{corner transition:} \quad x_1 = \dfrac{1}{2\lambda} -\dfrac{\cos\theta-\sin\theta_0}{\cos\theta_0-\sin\theta} \quad \text{and} \quad x_2=\dfrac{1}{2\lambda},\\ & \text{corner flow:} \quad x_1 =\dfrac{1}{2\lambda} -a\dfrac{\cos\theta-\sin\theta_0}{\cos\theta_0-\sin\theta} \quad \text{and} \quad x_2=\dfrac{1}{2\lambda}. \end{aligned} \right\} \end{equation}

The geometric function $B(\theta,\theta _0)$ in the expressions for $A$ is given by

(2.8)\begin{equation} B(\theta,\theta_0)=\theta+\theta_0 -{\rm \pi}/2-\cos\theta(\sin\theta-\cos\theta_0)+\cos\theta_0(\cos\theta-\sin\theta_0). \end{equation}

Equations (2.6) were also used by Kolliopoulos et al. (Reference Kolliopoulos, Jochem, Johnson, Suszynski, Francis and Kumar2021) and in other previous studies (Romero & Yost Reference Romero and Yost1996; Weislogel & Lichter Reference Weislogel and Lichter1998; Tchikanda, Nilson & Griffiths Reference Tchikanda, Nilson and Griffiths2004; Weislogel & Nardin Reference Weislogel and Nardin2005; Nilson et al. Reference Nilson, Tchikanda, Griffiths and Martinez2006; Yang & Homsy Reference Yang and Homsy2006; White & Troian Reference White and Troian2019).

2.3. Evaporation

We assume the liquid is in contact with a gas phase having ambient temperature $\hat {T}_A$ and relative humidity $R_H$. The gas phase consists of saturated vapour or a mixture of solvent vapour and inert gas (e.g. air), and its velocity is assumed to be negligible. We assume the gas phase density $\hat {\rho }_V$, viscosity $\hat {\mu }_V$, and thermal conductivity $\hat {k}_V$, are much smaller than their liquid counterparts (Burelbach, Bankoff & Davis Reference Burelbach, Bankoff and Davis1988), and the temperature across the liquid–air interface is continuous (Sefiane & Ward Reference Sefiane and Ward2007). These assumptions allow us to describe evaporation using a kinetically limited model focusing only on the liquid phase.

A kinetically limited model is chosen instead of a diffusion-limited model (Gambaryan- Roisman Reference Gambaryan-Roisman2019) for three reasons. First, a kinetically limited model is expected to be more accurate at describing evaporating flow of liquid solutions since the presence of solute at the free surface makes the rate-limiting step more likely to be in the liquid phase (Cazabat & Guéna Reference Cazabat and Guéna2010). Second, kinetically limited models have proven useful in interpreting experiments involving evaporation into an unsaturated environment when the vapour diffuses rapidly away from the evaporating interface even though they were originally developed for evaporation into a saturated environment (Murisic & Kondic Reference Murisic and Kondic2011). Third, a kinetically limited model does not require keeping track of the solvent concentration in the vapour phase, which makes it computationally simpler while still providing qualitatively accurate predictions (Ajaev Reference Ajaev2005; Murisic & Kondic Reference Murisic and Kondic2011). We note that, although a kinetically limited model is expected to provide qualitatively accurate predictions, it is unlikely to produce quantitatively accurate predictions for situations where evaporation is diffusion limited.

Evaporation is modelled using the non-equilibrium Hertz–Knudsen relation based on the kinetic theory of gases (Plesset & Prosperetti Reference Plesset and Prosperetti1976; Moosman & Homsy Reference Moosman and Homsy1980). The evaporative mass flux is described using

(2.9)\begin{equation} \hat{j}=\alpha_E\hat{p}_{V}\left(\frac{\hat{R}_g \hat{T}|_{\hat{h}}}{2 {\rm \pi}}\right)^{1/2}\left(\frac{\hat{p}_{V,e}}{\hat{p}_{V}}-1\right), \end{equation}

where $\alpha _E$ is the thermal accommodation coefficient, $\hat {R}_g$ is the gas constant per unit mass, $\hat {T}|_{\hat {h}}$ is the local liquid–air interface temperature, $\hat {p}_{V}$ is the partial pressure of solvent in the vapour phase and $\hat {p}_{V,e}$ is the equilibrium solvent vapour pressure. Note that we have assumed the thermal accommodation coefficients for evaporation and condensation are equal to $\alpha _E$. Physically, $\alpha _E$ can be thought of as a barrier to phase change, with $\alpha _E=1$ corresponding to no barrier and $\alpha _E=0$ corresponding to no phase change (Persad & Ward Reference Persad and Ward2016). Prior studies have reported values of $\alpha _E$ that vary over several orders of magnitude from ${O}(10^{-6})$ to ${O}(1)$ (Marek & Straub Reference Marek and Straub2001; Murisic & Kondic Reference Murisic and Kondic2011). In this work, we use the accommodation coefficient as a fitting parameter when comparing with experiments, similar to Murisic & Kondic (Reference Murisic and Kondic2011).

According to equilibrium thermodynamics (Moosman & Homsy Reference Moosman and Homsy1980; Ajaev & Homsy Reference Ajaev and Homsy2001; Ajaev Reference Ajaev2005) we can write

(2.10)\begin{equation} \ln\left(\frac{\hat{p}_{V,e}}{\hat{p}_V}\right)=\frac{\hat{\mathcal{L}}}{\hat{R}_g}\left(\frac{1}{\hat{T}_V}-\frac{1}{\hat{T}|_{\hat{h}}}\right)+\frac{\hat{p}-\hat{p}_V}{\hat{\rho} \hat{R}_g \hat{T}|_{\hat{h}}}, \end{equation}

where $\hat {\mathcal {L}}$ is the latent heat and $\hat {T}_V$ is the vapour temperature. The partial pressure of solvent in the vapour phase is calculated using $\hat {p}_V = R_H \hat {p}_S(\hat {T}_A)$ (Cazabat & Guéna Reference Cazabat and Guéna2010; Murisic & Kondic Reference Murisic and Kondic2011), where $R_H$ is the relative humidity, which ranges from 0 to 1, and $\hat {p}_S(\hat {T}_{A})$ is the saturation pressure corresponding to the ambient temperature $\hat {T}_{A}$. Both $\hat {p}_V$ and $\hat {T}_{A}$ in the vapour phase are assumed uniform and constant. The partial pressure of the solvent $\hat {p}_V$ is then used to calculate the vapour temperature $\hat {T}_V$ using the Clausius–Clapeyron equation (Murisic & Kondic Reference Murisic and Kondic2011). Note that for $R_H=1$ (which corresponds to a vapour phase saturated with solvent), $\hat {p}_V=\hat {p}_S$ and $\hat {T}_V=\hat {T}_A$.

We rescale (2.10) and (2.9) using

(2.11a–d)\begin{equation} \hat{p}=(\hat{\mu}_0 \hat{U}/\epsilon \hat{H})p, \quad \hat{T} = \hat{T}_V+T\Delta\hat{T}, \quad \Delta\hat{T}=\hat{T}_W-\hat{T}_V, \quad \text{and} \quad \hat{j}=(\hat{k}\Delta\hat{T}/\hat{\mathcal{L}}\hat{H})j, \end{equation}

where $\hat {T}_W$ is the temperature of the channel walls, and substitute (2.10) into (2.9), assuming $\ln (\hat {p}_{V,e}/\hat {p}_V)\approx \hat {p}_{V,e}/\hat {p}_V -1$ and $\sqrt {\hat {T}|_{\hat {h}}} \approx \sqrt {\hat {T}_V}$. The resulting expression for the evaporative mass flux is

(2.12)\begin{equation} Kj=\delta (p-p_V)+T|_h, \end{equation}

where $K=\hat {k}\sqrt {2{\rm \pi} \hat {R}_g^3\hat {T}_V^5}/\hat {p}_V\hat {\mathcal {L}}^2\hat {H}\alpha _E$ is the Knudsen number (ratio of interfacial to bulk heat transfer resistance), which is essentially the inverse of the Biot number, and $\delta =\hat {\mu }_0\hat {U}\hat {T}_V/\hat {\rho }\epsilon \hat {H}\hat {\mathcal {L}}\Delta \hat {T}$ accounts for the effects of pressure variation on the local interface temperature (Ajaev Reference Ajaev2005). From (2.12) it can be seen that the evaporative mass flux $j$ is proportional to deviations from $p_V$ and $\hat {T}_V$. Note that $\delta$ is typically ${O}(10^{-4}-10^{-6})$ compared with $T|_h$ which is ${O}(1)$, so contributions of the $\delta (p-p_V)$ term in (2.12) are typically neglected (Markos et al. Reference Markos, Ajaev and Homsy2006; Murisic & Kondic Reference Murisic and Kondic2011). However, these contributions become significant near the finger tip as $a\rightarrow 0$, where $p\rightarrow -\infty$ as seen in (2.6d).

2.4. Energy transport

Energy conservation is governed by

(2.13)\begin{equation} \hat{\rho}\hat{C}_p\left[\frac{\partial \hat{T}}{\partial \hat{t}} + \boldsymbol{\hat{u}} \boldsymbol{\cdot} \boldsymbol{\widehat{\nabla}}\hat{T}\right] = \hat{k} \widehat{\nabla}^2 \hat{T}. \end{equation}

For simplicity, we consider the limiting case where heat conduction in the channel walls is neglected, and assume that the channel walls are held at a constant temperature $\hat {T}_W$. The energy balance at the liquid–air interface is

(2.14)\begin{equation} \hat{\mathcal{L}}\hat{j}={-}\boldsymbol{n} \boldsymbol{\cdot} \hat{k} \boldsymbol{\widehat{\nabla}}\hat{T}. \end{equation}

We render (2.13) dimensionless using the scalings in (2.3) and (2.11a–d). The dimensionless numbers that arise are the Reynolds number $Re$ (see § 2.2) and the Prandtl number $Pr=\hat {\mu }_0\hat {C}_p/\hat {k}$ (ratio of momentum to thermal diffusivity).

In the limits of $\epsilon ^2\ll 1$ and $\epsilon Re Pr\ll 1$, the leading-order energy-transport equation is

(2.15a)\begin{equation} T_{xx} + T_{yy} = 0, \end{equation}

subject to

(2.15b)\begin{equation} T = 1 \quad \text{at solid boundaries} \quad \text{and} \quad j (1+h_x^2)^{1/2} = h_xT_x - T_y \quad \text{at free surface}, \end{equation}

where $j$ is determined using (2.12). The total evaporative mass flux for a given channel cross-section is defined as

(2.16)\begin{equation} \tilde{J}=\int_S j\,\mathrm{d}S = 2\int_{x_1}^{x_2} j (1+h_x^2)^{1/2} \,\mathrm{d} x, \end{equation}

where $x_1$ and $x_2$ are given by (2.7). The liquid–air interface arc length is given by

(2.17)\begin{equation} S = \begin{cases} \dfrac{{\rm \pi} - 2\theta}{\lambda\cos\theta} \quad \text{in meniscus-deformation regime}, \\ \dfrac{{\rm \pi} - 2\theta_0}{\lambda\cos\theta_0} \quad \text{in meniscus-recession regime}, \\ \dfrac{(2{\rm \pi} - 4\theta-4\beta)\sin\beta}{\cos(\theta+\beta)} \quad \text{in corner-transition regime}, \\ \dfrac{({\rm \pi} - 4\theta)a\sin{\rm \pi}/4}{\cos(\theta+{\rm \pi}/4)} \quad \text{in corner-flow regime}, \end{cases} \end{equation}

where $\theta$ and $a$ are dependent on $(z,t)$ and the cross-sectional-averaged dimensionless temperature is defined as

(2.18)\begin{equation} \bar{T}=\frac{1}{A}\int_A T\,\mathrm{d}A, \end{equation}

where $A$ is the liquid cross-sectional area from (2.6).

2.5. Solute transport

A convection–diffusion equation governs the transport of solute

(2.19a)\begin{equation} \frac{\partial c}{\partial \hat{t}} + \boldsymbol{\hat{u}} \boldsymbol{\cdot} \boldsymbol{\widehat{\nabla}} c= \hat{D} \widehat{\nabla}^2 c, \end{equation}

where $c$ is the solute concentration (mass fraction) and $\hat {D}$ is the diffusion coefficient, which is assumed to be constant. We impose no-flux boundary conditions

(2.19b)\begin{gather} \hat{D}(\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\widehat{\nabla}})c = 0 \quad \text{at solid boundaries}, \end{gather}

(2.19c)\begin{gather} \hat{D}(\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\widehat{\nabla}})c = c(\boldsymbol{\hat{u}}-\boldsymbol{\hat{u}_I}) \boldsymbol{\cdot} \boldsymbol{n} = c \hat{j}/\hat{\rho} \quad \text{at free surface}, \end{gather}

where $\boldsymbol {\hat {u}_I}$ is the liquid–air interface velocity and $\hat {j}$ is the local evaporative mass flux. Using the scalings in (2.3), (2.19a) becomes

(2.20a)\begin{equation} \epsilon^2 Pe (c_t + uc_x + v c_y + wc_z) = c_{xx} + c_{yy}+\epsilon^2 c_{zz}, \end{equation}

where $Pe=\hat {U}\hat {L}/\hat {D}$ is the Péclet number (ratio of convective to diffusive transport rates). Similarly, the no-flux boundary conditions become

(2.20b)\begin{gather} (\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\nabla})c = 0 \quad \text{at solid boundaries}, \end{gather}

(2.20c)\begin{gather} -h_x c_x + c_y - \epsilon^2 h_z c_z = \epsilon^2 Pe cEj(1+h_x^2)^{1/2} \quad \text{at free surface}, \end{gather}

where $E=\hat {k}\Delta \hat {T}/\hat {\rho }\hat {\mathcal {L}}\hat {H}\epsilon \hat {U}$ is the evaporation number (ratio of characteristic capillary to evaporation times).

To simplify (2.20a) further, we assume solute transport in the $x$ and $y$ directions is dominated by diffusion ($\epsilon ^2 Pe \ll 1$). As proposed by Jensen & Grotberg (Reference Jensen and Grotberg1993), we asymptotically expand the concentration in terms of $\epsilon ^2 Pe$, obtaining $c(x,y,z,t) = c_0(z,t) + \epsilon ^2 Pe c_1(x,y,z,t)$, where we assume $\int _A c_1 \, \mathrm {d}A=0$. This allows us to define the cross-sectional-averaged concentration as

(2.21)\begin{equation} \bar{c}= \frac{1}{A}\int_A c \,\mathrm{d}A = c_0(z,t). \end{equation}

We apply cross-sectional averaging to (2.20a) and replace the $c_1$ terms using the no-flux condition at the free surface in (2.20c). In the limit of $\epsilon ^2 \ll 1$, we obtain the following evolution equation for $\bar {c}$:

(2.22)\begin{equation} \bar{c}_{t} + \bar{w} \bar{c}_{z} = \frac{1}{PeA}(A\bar{c}_{z})_z + \frac{\bar{c}}{A}E\tilde{J}, \end{equation}

where $\bar {w}$ is the cross-sectional-averaged velocity (which will be obtained from (2.28)), $A$ is the dimensionless cross-sectional area from (2.6), and $\tilde {J}$ is the total cross-sectional evaporative mass flux in a given channel cross-section from (2.16).

2.6. Constitutive equations for viscosity and surface tension

The constitutive equations for viscosity $M$ and surface tension $\varSigma$ depend on the liquid solution we choose to study. In this work, we use aqueous poly(vinyl alcohol) (PVA) solutions and compare model predictions with capillary-flow experiments conducted by Lade et al. (Reference Lade, Jochem, Macosko and Francis2018). An empirical model proposed by Patton (Reference Patton1964) is used to capture the dependence of the viscosity on $\bar {T}$ and $\bar {c}$ through

(2.23a)\begin{equation} \log M = \dfrac{\bar{c}}{k_a(\bar{T})+\bar{c}k_b(\bar{T})}, \end{equation}

where

(2.23b)\begin{equation} \left. \begin{aligned} k_a(\bar{T}) & = 1.28\times10^{{-}5}(\hat{T}_V+\bar{T}\Delta\hat{T})+1.59\times10^{{-}2},\\ k_b(\bar{T}) & = 3.83\times10^{{-}4}(\hat{T}_V+\bar{T}\Delta\hat{T})-2.47\times10^{{-}2}. \end{aligned} \right\} \end{equation}

The $k_a(\bar {T})$ and $k_b(\bar {T})$ functions were reported by Lade et al. (Reference Lade, Jochem, Macosko and Francis2018) after fitting the empirical model to rheological data of PVA solutions for a range of temperatures and concentrations. In (2.23a), increasing the solute concentration can increase the viscosity by orders of magnitude, and increasing the temperature decreases the viscosity but does not change its order of magnitude. Similar models can be used to describe any solution or colloidal suspension where the shear viscosity is the dominant rheological parameter.

The effects of $\bar {T}$ and $\bar {c}$ on the surface tension are modelled using

(2.24)\begin{equation} \varSigma = 1 - Ma_c\bar{c}^{1/2} - Ma_T \bar{T}, \end{equation}

where $Ma_c=\hat {\gamma }_c \epsilon /\hat {\mu }_0 \hat {U}$ is the solutal Marangoni number and $Ma_T=\hat {\gamma }_T\Delta \hat {T}\epsilon /\hat {\mu }_0 \hat {U}$ is the thermal Marangoni number, which are ratios of surface-tension-gradient forces to viscous forces, and $\hat {\gamma }_c$ and $\hat {\gamma }_T$ are experimentally obtained constants. We assume the temperature at the liquid–air interface does not deviate much from the vapour temperature, which allows us to write the surface tension as a linear function of $\bar {T}$ (Burelbach et al. Reference Burelbach, Bankoff and Davis1988; Gramlich et al. Reference Gramlich, Kalliadasis, Homsy and Messer2002; Ajaev Reference Ajaev2005; Craster, Matar & Sefiane Reference Craster, Matar and Sefiane2009). Many prior studies assume a dilute solution and use a linearized surface-tension dependence on the concentration (e.g. Lam & Benson Reference Lam and Benson1970; Pham, Cheng & Kumar Reference Pham, Cheng and Kumar2017); this would yield results that are qualitatively similar to those obtained using (2.24). Comparison of the empirical models in (2.23a) and (2.24) with the experimental results of Lade et al. (Reference Lade, Jochem, Macosko and Francis2018) is found in the supplementary material available at https://doi.org/10.1017/jfm.2022.140.

2.7. Liquid height evolution

We begin with the no-flux boundary condition at the free surface given by

(2.25)\begin{equation} (\boldsymbol{\hat{u}}-\boldsymbol{\hat{u}_I}) \boldsymbol{\cdot} \boldsymbol{n} = \hat{j}/\hat{\rho}, \end{equation}

where the velocity of the liquid–air interface is $\boldsymbol {\hat {u}_I}=(0,\hat {h}_{\hat {t}},0)$. Using the scalings in (2.3) and (2.11a–d), (2.25) becomes

(2.26)\begin{equation} {-}uh_x+v-wh_z-h_t=Ej(1+h_x^2+\epsilon^2h_z^2)^{1/2}. \end{equation}

We apply cross-sectional averaging to the mass conservation equation (2.4a) and replace the $u$ and $v$ terms using (2.26), thus obtaining

(2.27)\begin{equation} A_t={-}(\bar{w} A)_z - E\tilde{J}, \end{equation}

where for each regime, the dimensionless liquid cross-sectional area $A=2\int _{x_1}^{x_2} h\,\mathrm {d} x$ is given by (2.6), $\bar {w}=A^{-1}\int _A w\,\mathrm {d}A$ is the cross-sectional-averaged velocity and $\tilde {J}$ is the total evaporative mass flux in a given channel cross-section from (2.16). Equation (2.27) is the mass-balance equation derived by Lenormand & Zarcone (Reference Lenormand and Zarcone1984), relating the time derivative of the dimensionless liquid cross-sectional area $A$ to the gradient in the dimensionless flux $Q=\int _A w\,\mathrm {d}A =\bar {w} A$, with an additional term accounting for mass lost due to solvent evaporation.

The velocity in (2.27) can be decomposed as follows:

(2.28)\begin{equation} \bar{w} = \bar{w}_{ca} + \bar{w}_{cg} + \bar{w}_{tg}, \end{equation}

where each contribution in each regime (figure 1) is expressed in (2.29). Each component corresponds to a different mechanism acting on the liquid to drive flow. Here, $\bar {w}_{ca}$ is the velocity due to capillary effects, while $\bar {w}_{cg}$ and $\bar {w}_{tg}$ correspond to the effects of Marangoni stresses. Specifically, $\bar {w}_{cg}$ is due to solute concentration gradients and $\bar {w}_{tg}$ is due to thermal gradients. For each regime these contributions are given by

(2.29a)$$\begin{gather} \bar{w}_{ca} ={-}\dfrac{\bar{U}^c_{D}}{M}p_z, \quad \bar{w}_{cg} = \dfrac{\bar{U}^g_{D}}{M}\dfrac{\partial\varSigma}{\partial\bar{c}}\bar{c}_z, \quad \bar{w}_{tg} = \dfrac{\bar{U}^g_{D}}{M}\dfrac{\partial\varSigma}{\partial\bar{T}}\bar{T}_z, \quad \text{meniscus deformation}, \end{gather}$$

(2.29b)$$\begin{gather}\bar{w}_{ca} ={-}\dfrac{\bar{U}^c_{R}}{M}p_z, \quad \bar{w}_{cg} = \dfrac{\bar{U}^g_{R}}{M}\dfrac{\partial\varSigma}{\partial\bar{c}}\bar{c}_z, \quad \bar{w}_{tg} = \dfrac{\bar{U}^g_{R}}{M}\dfrac{\partial\varSigma}{\partial\bar{T}}\bar{T}_z, \quad \text{meniscus recession}, \end{gather}$$

(2.29c)$$\begin{gather}\bar{w}_{ca} ={-}\dfrac{\bar{U}^c_{T}}{M}p_z, \quad \bar{w}_{cg} = \dfrac{\bar{U}^g_{T}}{M}\dfrac{\partial\varSigma}{\partial\bar{c}}\bar{c}_z, \quad \bar{w}_{tg} = \dfrac{\bar{U}^g_{T}}{M}\dfrac{\partial\varSigma}{\partial\bar{T}}\bar{T}_z, \quad \text{corner transition}, \end{gather}$$

(2.29d)$$\begin{gather}\bar{w}_{ca} ={-}a^2\dfrac{\bar{U}^c_{C}}{M}p_z, \quad \bar{w}_{cg} = \dfrac{\bar{U}^g_{C}}{M}\dfrac{\partial\varSigma}{\partial\bar{c}}\bar{c}_z, \quad \bar{w}_{tg} = a\dfrac{\bar{U}^g_{C}}{M}\dfrac{\partial\varSigma}{\partial\bar{T}}\bar{T}_z, \quad \text{corner flow}, \end{gather}$$

where $\bar {U}^c_{i}$ and $\bar {U}^g_{i}$ are rescaled cross-sectional-averaged dimensionless velocities, with the subscript $i$ being equal to $D,R,T$ or $C$ for the meniscus-deformation, meniscus-recession, corner-transition and corner-flow regimes, respectively. Details of the calculation of $\bar {U}^c_{i}$ and $\bar {U}^g_{i}$ can be found in the supplementary material. The expressions for $M$ and $\varSigma$ are given by (2.23a) and (2.24), respectively.

Consistent with prior studies considering horizontal rectangular channels, we neglect the meniscus-recession regime (i.e. $z_D=z_M$) where $\bar {w}_{ca}=0$. This is because the transverse curvature gradients are zero (constant $p$ in (2.6b)) and the only contribution to $\bar {w}_{ca}$ is from ${O}(\epsilon ^2)$ axial curvature gradients, which we did not account for. The transition from the meniscus-deformation regime to the corner-flow regime (for $\lambda >\lambda _c$) is treated as a jump in the dimensionless liquid height $a(z,t)$ (Nilson et al. Reference Nilson, Tchikanda, Griffiths and Martinez2006; Kolliopoulos et al. Reference Kolliopoulos, Jochem, Johnson, Suszynski, Francis and Kumar2021).

2.8. Reservoir

Liquid is supplied to the microchannel from a cylindrical reservoir of radius $\hat {R}$ and height $\hat {H}$, and the liquid–air interface in the reservoir is assumed to be axisymmetric. The channel inlet is assumed to have a negligible influence on the liquid–air interface in the reservoir. The reservoir aspect ratio is $\lambda _R=\hat {H}/\hat {R}$ and the cylindrical coordinates of the reservoir are scaled using $(\hat {r},\hat {y},\phi )=(\hat {R}r,\hat {H} y,\phi )$.

2.8.1. Hydrodynamics

Following a similar procedure to the one described in § 2.2, the normal stress balance in (2.1b) reduces to the Young–Laplace equation,

(2.30)\begin{equation} p-p_{V,T}={-}\frac{\varSigma \lambda_R^2}{Ca} \left[\frac{h_{rr}}{(1+\lambda_R^2h_{r}^2)^{3/2}}+\frac{h_r}{r(1+\lambda_R^2h_{r}^2)^{1/2}}\right]. \end{equation}

Using (2.30) in combination with a condition for the contact-line location on the solid wall ($h=1$, at $r=1$), a symmetry condition at the reservoir centre ($h_r=0$, at $r=0$) and the definition of the contact angle $\theta _R$ on the solid wall ($\lambda _Rh_r/[1+\lambda _R^2 h_r^2]^{1/2}=\cos \theta _R$, at $r=1$), we obtain the following leading-order expressions for pressure $p(t)$ and liquid–air interface profile $h(r,t)$ in the reservoir:

(2.31a)\begin{gather} p={-}2\varSigma\lambda_R\cos\theta_R(t)+p_{V,T}, \end{gather}

(2.31b)\begin{gather} h=1+\dfrac{\tan\theta_R(t)}{\lambda_R}-\frac{1}{\lambda_R}\left[\dfrac{1}{\cos^2\theta_R(t)}-r^2\right]^{1/2}. \end{gather}

2.8.2. Energy transport

Similar to § 2.4, we consider energy conservation in the limit of $\lambda _R Re Pr\ll 1$, resulting in

(2.32a)\begin{equation} \frac{\lambda_R^2}{r}\left(rT_r\right)_r+T_{yy}=0, \end{equation}

subject to

(2.32b)\begin{equation} T = 1 \quad \text{at solid boundaries}, \quad \text{and} \quad j (1+\lambda_R^2h_r^2)^{1/2} = \lambda_R^2 h_rT_r - T_y \quad \text{at free surface}. \end{equation}

Note that in the limit of $\lambda _R\rightarrow 0$, (2.32) solved along with (2.12) results in the evaporation models used for axisymmetric droplets and thin films, where $j = (1+\delta (p-p_V))/(K + h)$ (Ajaev & Homsy Reference Ajaev and Homsy2001; Ajaev Reference Ajaev2005; Pham & Kumar Reference Pham and Kumar2017, Reference Pham and Kumar2019).

The dimensionless total evaporative mass flux for a given reservoir cross-section is defined as

(2.33)\begin{equation} \tilde{J}_R=\int_{S_R} j\,\mathrm{d}S_R=\int_0^1 j (1+\lambda_R^2h_r^2)^{1/2}r\,\mathrm{d}r. \end{equation}

The liquid–air interface arc length of the reservoir cross-section is given by

(2.34)\begin{equation} S_R = \dfrac{1}{2}+\dfrac{1}{2}\left(\dfrac{1-\sin\theta_R(t)}{\cos\theta_R(t)}\right)^2. \end{equation}

2.8.3. Liquid volume evolution

The dimensionless liquid volume in the reservoir $V_R=2\int hr\,\mathrm {d}r$ is

(2.35)\begin{equation} V_R = 1 - \frac{1-\sin\theta_R(t)}{6\lambda_R \cos\theta_R(t)}\left[3+\left(\frac{1-\sin\theta_R(t)}{\cos\theta_R(t)}\right)^2\right], \end{equation}

where $V_R$ is scaled by the reservoir volume ${\rm \pi} \hat {R}^2 \hat {H}$. The ratio of channel to reservoir volume is $f_R = \lambda _R^2/ {\rm \pi}\epsilon \lambda$. The reservoir is considered depleted when the liquid–air interface contacts the reservoir bottom (i.e. $h(r=0)=0$ which corresponds to $\theta _R=\arcsin ((1-\lambda _R^2)/(1+\lambda _R^2))$ and $V_R=(3-{\rm \pi} \epsilon \lambda f_R)/6$).

We model the evolution of $V_R$ through the following total mass balance:

(2.36)\begin{equation} \left(V_R\right)_t=\left.-f_R\lambda\bar{w} A\right|_{z=0}-2E\tilde{J}_R, \end{equation}

where the rate of change of liquid volume in the reservoir is equal to the liquid flux into the channel plus the liquid lost to evaporation. Evaporation is accounted for through the total evaporative mass flux for a given reservoir cross-section $\tilde {J}_R$ using (2.33).

2.8.4. Solute transport

The solute concentration in the reservoir $c_R$ is assumed to be spatially uniform and its evolution is modelled using the following species mass balance:

(2.37)\begin{equation} \left(c_RV_R\right)_t=\left.-f_R\lambda\bar{w} \bar{c}A\right|_{z=0}, \end{equation}

where the rate of change of solute mass in the reservoir is equal to the solute mass flux into the channel.

3.1. Boundary conditions

At the channel inlet ($z=0$; see figure 1), the pressure and solute concentration are matched to the reservoir pressure and solute concentration, respectively, and the liquid height is assumed to be pinned to the top of the channel sidewall. This results in the following conditions at the channel inlet:

(3.1a–c)\begin{equation} \theta(0,t) = \cos^{{-}1}\left[\lambda_R\cos\theta_R(t)/\lambda\right], \quad a(0,t) = 1, \quad \text{and} \quad \bar{c}(0,t) = c_R(t). \end{equation}

Note that the condition on the contact angle at the channel inlet is obtained by matching the pressure in (2.6a) and (2.31a).

For $\lambda \geq \lambda _c$, at the transition from the meniscus-deformation to corner-flow regime ($z=z_M$), we impose the following boundary conditions:

(3.2)\begin{equation} \left. \begin{aligned} \theta(z_M^-,t) & =\theta(z_M^+,t)=\theta_0, \quad a(z_M^-,t)=1, \quad a(z_M^+,t)=\lambda_c/\lambda,\\ \bar{c}(z_M^-,t) & =\bar{c}(z_M^+,t), \quad \text{and} \quad \left. A\bar{c}_z\right|_{z=z_M^-}=\left.A\bar{c}_z\right|_{z=z_M^+}, \end{aligned} \right\} \end{equation}

where the boundary conditions on $\theta$ and $a$ are discussed in § 2.1 and the boundary conditions on $\bar {c}$ physically represent mass continuity with no accumulation at the interface between the two regimes.

For $\lambda <\lambda _c$, at the transition from the meniscus-deformation to corner-transition regime ($z=z_M$), we impose the following boundary conditions

(3.3)\begin{equation} \left. \begin{aligned} \theta(z_M^-,t) & =\theta_C, \quad\theta(z_M^+,t)=\theta_T , \quad a(z_M^-,t)=a(z_M^+,t)=1,\\ \bar{c}(z_M^-,t) & =\bar{c}(z_M^+,t), \quad \text{and} \quad \left.A\bar{c}_z\right|_{z=z_M^-}=\left.A\bar{c}_z\right|_{z=z_M^+}. \end{aligned} \right\} \end{equation}

Here, $\theta _C= \arcsin [(1-4\lambda ^2)/(1+4\lambda ^2)]$ is the critical angle at which the upper meniscus touches the channel bottom and $\theta _T$ (§ 2.1) is the angle determined (via Newton's method) by setting $A(z_M^-,t)=A(z_M^+,t)$ to conserve mass. At the transition from the corner-transition to corner-flow regime ($z=z_C$), we impose the following boundary conditions:

(3.4)\begin{equation} \left. \begin{aligned} \theta(z_C^-,t) & =\theta(z_C^+,t)=\theta_0 , \quad a(z_C^-,t)=a(z_C^+,t)=1,\\ \bar{c}(z_C^-,t) & =\bar{c}(z_C^+,t), \quad \text{and} \quad \left. A\bar{c}_z\right|_{z=z_C^-}=\left.A\bar{c}_z\right|_{z=z_C^+}. \end{aligned} \right\} \end{equation}

Finally, at the end of the corner-flow regime ($z=z_T$) we impose

(3.5a–c)\begin{equation} \theta(z_T,t) = \theta_0, \quad a(z_T,t) = 0, \quad \text{and} \quad \bar{c}_z(z_T,t)=0. \end{equation}

The boundary condition on $\bar {c}$ in (3.5a–c) corresponds to no flux. The boundary conditions in (3.2)–(3.5a–c) for the contact angle $\theta$ and the liquid height $a$ on the channel sidewall were also used by Kolliopoulos et al. (Reference Kolliopoulos, Jochem, Johnson, Suszynski, Francis and Kumar2021).

3.2. Initial conditions

Initial conditions for $\theta (z,t_0)$, $a(z,t_0)$, $z_M(t_0)$, $z_C(t_0)$ and $z_T(t_0)$ are generated using the similarity solutions reported by Kolliopoulos et al. (Reference Kolliopoulos, Jochem, Johnson, Suszynski, Francis and Kumar2021) in the absence of evaporation ($E=0$). The reported solutions are in terms of the self-similar variable $\eta =z/\sqrt {t}$, so we determine the initial interface profile and its axial coordinates $z=\eta \sqrt {t_0}$ by setting $t_0=10^{-4}$. The chosen $t_0$ does not influence our results since total flow times are ${O}(10-10^2)$. Additionally, we assume the reservoir is initially fully filled and the solute concentration is uniform in the reservoir and channel. Hence,

(3.6a–c)\begin{equation} V_R(t_0) = 1, \quad c_R(t_0) = C_0, \quad \text{and} \quad \bar{c}(z,t_0) = C_0. \end{equation}

3.3. Solution procedure

The rescaled cross-sectional-averaged dimensionless velocities $\bar {U}_i^c$ and $\bar {U}_i^g$ in (2.29) are calculated as discussed in the supplementary material. Velocity fields are numerically solved for with a Galerkin finite-element method using quadratic basis functions (Kolliopoulos et al. Reference Kolliopoulos, Jochem, Johnson, Suszynski, Francis and Kumar2021). Results for $\bar {U}_D^c$ and $\bar {U}_T^c$ are found to be in agreement with results of Tchikanda et al. (Reference Tchikanda, Nilson and Griffiths2004) and Weislogel & Nardin (Reference Weislogel and Nardin2005), respectively. Results for $\bar {U}_C^c$ agree with results of Ayyaswamy, Catton & Edwards (Reference Ayyaswamy, Catton and Edwards1974), Ransohoff & Radke (Reference Ransohoff and Radke1988) and Yang & Homsy (Reference Yang and Homsy2006). Additionally, results for $\bar {U}_D^g$ and $\bar {U}_C^g$ are in agreement with results of Tchikanda et al. (Reference Tchikanda, Nilson and Griffiths2004) and Yang & Homsy (Reference Yang and Homsy2006), respectively.

Rather than consider effects of $a(z,t)$ and $\theta (z,t)$ on $\bar {U}_i^c$ and $\bar {U}_i^g$ separately, we consider the dependence of $\bar {U}_i^c$ and $\bar {U}_i^g$ on the liquid saturation $\lambda A$ (ratio of channel cross-sectional area filled with liquid to total channel cross-sectional area). These rescaled cross-sectional-averaged dimensionless velocities $\bar {U}_i^c(\lambda A)$ and $\bar {U}_i^g(\lambda A)$ are fitted with Chebyshev polynomials of the first kind using the least-squares method and then applied to calculate the cross-sectional-averaged dimensionless velocity components of $\bar {w}$ in (2.29).

The dimensionless total evaporative mass fluxes and cross-sectional-averaged dimensionless temperatures for the channel and reservoir are calculated as discussed in § 2.4 and § 2.8.2, respectively. Temperature fields are numerically solved for with a Galerkin finite-element method using quadratic basis functions. Results for $\tilde {J}$ in the corner-flow regime agree with results of Markos et al. (Reference Markos, Ajaev and Homsy2006). Results for $\tilde {J}(\lambda A)$ and $\bar {T}(\lambda A)$ in each regime, and $\tilde {J}_R(V_R)$ and $\bar {T}(V_R)$ for the reservoir, are fitted with Chebyshev polynomials of the first kind using the least-squares method and then applied in (2.27), (2.22), (2.36) and (2.37).

For $\lambda \ge \lambda _c$, the system of equations consists of (2.27) and (2.22) for each regime in the channel (meniscus-deformation and corner-flow regimes), and (2.36) and (2.37) for the reservoir. Additionally, the positions of the moving regime boundaries are determined using the global continuity equation assuming no accumulation at the interface between two regimes, as discussed by Kolliopoulos et al. (Reference Kolliopoulos, Jochem, Johnson, Suszynski, Francis and Kumar2021). The general form of the global continuity equation at the interface at position $z_i$ is

(3.7)\begin{equation} \left[\int_A \boldsymbol{n} \boldsymbol{\cdot} (\boldsymbol{u}-\boldsymbol{u_I}) \,\mathrm{d}A \right]_{z=z_i^-}=\left[\int_A \boldsymbol{n} \boldsymbol{\cdot} (\boldsymbol{u}-\boldsymbol{u_I}) \,\mathrm{d}A \right]_{z=z_i^+} \quad \text{at } z = z_i, \end{equation}

where $\boldsymbol {n}=(0,0,1)$ is the unit normal to the interface, $\boldsymbol {u} = (u,v,w)$ is the dimensionless liquid velocity, $\boldsymbol {u_I} = (0, 0, \textrm {d}z_i/\textrm {d}t)$ is the interface velocity at position $z_i$ and $A$ is the liquid cross-sectional area at position $z_i$. For the meniscus position $z_M$ and finger tip position $z_T$, (3.7) reduces to

(3.8a)\begin{equation} \left[A\bar{w} - A\frac{{\rm d}z_M}{{\rm d}t}\right]_{z=z_M^-} = \left[A\bar{w} - A\frac{{\rm d}z_M}{{\rm d}t}\right]_{z=z_M^+} \quad \text{at } z = z_M, \end{equation}

and

(3.8b)\begin{equation} \left.\bar{w}\right|_{z=z_T^-} = \frac{{\rm d}z_T}{{\rm d}t} \quad \text{at } z = z_T. \end{equation}

For $\lambda <\lambda _c$, the system of equations consists of (2.27) and (2.22) for each regime in the channel (meniscus-deformation, corner-transition and corner-flow regimes), along with (2.36) and (2.37) for the reservoir. In addition to the global continuity equations in (3.8) to determine $z_M$ and $z_T$, we use

(3.9)\begin{equation} \left[A\bar{w} - A\frac{{\rm d}z_C}{{\rm d}t}\right]_{z=z_C^-} = \left[A\bar{w} - A\frac{{\rm d}z_C}{{\rm d}t}\right]_{z=z_C^+} \quad \text{at } z = z_C, \end{equation}

to determine $z_C$.

The spatial derivatives in both systems of equations are approximated with second-order centred finite differences. The resulting discretized systems of ordinary differential equations are solved using the fully implicit, variable-step and variable-order ode15i solver in MATLAB.

Numerical solutions for the dimensionless total evaporative mass flux $\tilde {J}(\lambda A)$ and cross-sectional-averaged dimensionless temperature $\bar {T}(\lambda A)$ for the channel and their dependence on the channel aspect ratio $\lambda$ and equilibrium contact angle $\theta _0$ are presented first (§ 4.1). Using these numerical solutions we then consider capillary flow of an evaporating pure solvent (§ 4.2) and liquid solution (§ 5) for $\lambda \ge \lambda _c$ and $\lambda <\lambda _c$. Finally, model predictions for aqueous PVA solutions are compared with experimental results of Lade et al. (Reference Lade, Jochem, Macosko and Francis2018) (§ 6).