Dynamic contact angle of a liquid spreading on an unsaturated wettable porous substrate

The modelling of the spreading of liquids over porous substrates in the framework of continuum mechanics requires, and is based on, the separation of scales between the ‘macroscopic’ (or ‘Darcy-scale’) and ‘microscopic’ (or ‘pore-scale’) processes (Barenblatt et al., 1990). In the continuum approximation, the description of the spreading phenomenon brings in the notion of an ‘effective’ smooth penetrable solid substrate, which is how the actual porous medium is represented, together with the notions of ‘effective’ contact lines and ‘effective’ contact angles that the free surface of the pure fluid above and the wetting front inside the substrate form with it. These notions, being averages in the sense of mechanics of multiphase media (Whitaker, 1999), are fundamentally different from the concepts of ‘contact line’ and ‘contact angle’ used in the modelling of dynamic wetting on the actual, as opposed to ‘effective’, impermeable solid surfaces. In particular, the notion of a ‘static’ (or ‘equilibrium’) contact angle, central to the modelling of dynamic wetting where it is a measure of wettability of the solid (Ralston et al., 2008), becomes meaningless for the Darcy-scale description of the liquid spreading over a porous medium: as experiments show (e.g.  Clarke et al., 2002; Starov et al., 2003; Markicevic et al., 2010), if, say, a drop of liquid is deposited onto an unsaturated wettable porous substrate, the eventual equilibrium state will be the drop completely imbibed into the porous medium, with no liquid left above it and hence no ‘effective’ static (equilibrium) contact angle between the, now non-existent, free surface of the pure liquid and the effective surface of the substrate.

The main theoretical implication of this absence of a meaningful effective static (equilibrium) contact angle is that the effective dynamic contact angle that the free surface of a spreading liquid forms with a porous substrate on the Darcy-scale cannot be regarded as essentially a perturbation of the static contact angle, which is what one invariably finds in all models dealing with the dynamic wetting of impermeable solid surfaces (see, Dussan V., 1979; de Gennes, 1985; Blake, 2006; Shikhmurzaev, 2011, for reviews). Therefore, it becomes important to consider the problem from first principles, without implying a-priori that the case of a porous substrate can be described by adjusting the concepts borrowed from the modelling of dynamic wetting of impermeable solid surfaces, such as the equilibrium contact angle.

In the present paper, we show that, unlike the situation one has in dynamic wetting, where the ‘microscopic’ dynamic contact angle has to be specified as an additional boundary condition³³3Here we are not discussing the so-called ‘apparent’ contact angle resulting from the free-surface bending near the contact line. This angle is not part of the mathematical problem formulation; it is an auxiliary concept introduced in some works and in different ways (see, e.g., Wilson-etal06) to interpret experimental results when the measured angle differs from the angle imposed in the problem formulation., for the process of the fluid spreading over a porous substrate one can find the ‘effective’ dynamic contact angle from the requirement of the absence of a nonintegrable singularity of the fluid’s pressure at the contact line, i.e. essentially from the requirement that a solution exists. (The flow inside the porous matrix is treated in the standard way, i.e. employing a dynamic wetting model for the impermeable solid surface to describe the propagation of menisci in the pores.) The limits confining the considered regime suggest certain experimentally verifiable predictions of the model which we briefly discuss in the light of the available experimental data.

Consider an incompressible Newtonian fluid of density $\rho$ and viscosity $\mu$ spreading at a speed $U$ over an isotropic homogeneous unsaturated porous substrate characterized by an effective pore size $a$. The gas displaced by the fluid from the substrate and from the inside of the porous matrix is assumed to be ideal and at a constant pressure with respect to which the pressure in the pure fluid and in the fluid inside the porous medium will be measured. In order to be able to model the process in the framework of continuum mechanics, we need a separation of scales between the macroscopic and the pore-scale phenomena, i.e. we have to consider the continuum limit

$$\epsilon=\frac{a}{L}\to 0,$$

(1)

where $L$ is the characteristic length scale on which the phenomenon is described. The resulting model will be applicable to experiments if $a/L\ll 1$, with the actual value of $a/L$ determining its accuracy. In the 0th approximation in the above limit, one has (see Fig. 1) (a) a macroscopic ‘wetting front’ as a sharp interface $S_{w}$ separating the saturated porous medium $\Omega_{2}$ from the unsaturated matrix $\Omega_{3}$, (b) an effective ‘contact line’ at which the free surface $S_{f}$, confining the domain $\Omega_{1}$ occupied by the pure fluid, and the wetting front meet, and (c) well-defined ‘contact angles’ $\theta_{D}$ and $\theta_{w}$ that the free surface and the wetting front form with the ‘effective’ surface $S_{0}$ of the solid. The reference frame in which the problem will be considered and the directions of unit normals $\mathbf{n}$, $\tilde{\mathbf{n}}$ and $\mathbf{n}_{0}$ to, respectively, $S_{w}$, $S_{f}$ and $S_{0}$ are shown in the figure. Importantly, in the scheme outlined above we have already made an assumption that the two contact lines, i.e. the contact line CL1 formed with the substrate by the free surface $S_{f}$ and the contact line CL2 formed with it by the wetting front $S_{w}$, coincide. This is always what happens when the fluid is first brought in contact with the substrate, and we will examine what follows from this initial situation and, later, the conditions when the assumption that CL1 and CL2 coincide no longer holds.

We will consider the case of small Reynolds and capillary numbers for the flow of the pure fluid, i.e. the limit $Re=\rho LU/\mu\to 0$ and

$$Ca=\frac{\mu U}{\sigma}\to 0,$$

(2)

where $\sigma$ is the surface tension of the fluid-gas interface. In this limit, to leading order, inertial effects can be neglected, and, as in the case of the dynamic wetting of an impermeable solid surface (Huh & Scriven, 1971; Shikhmurzaev, 1993), from the normal-stress boundary condition on the free surface we have that near the contact line the free surface ($S_{f}$) is locally planar, so that locally $\Omega_{1}$ is a wedge region. For simplicity we will also neglect gravity, though its inclusion would not change the main results of the analysis below.

Scaling the pressure and velocity in $\Omega_{1}$ with $\mu U/L$ and $U$ respectively, one has that in $\Omega_{1}$ the dimensionless pressure $\tilde{p}$ and velocity $\tilde{\mathbf{u}}$ obey the Stokes equations

$$\nabla\cdot\tilde{\mathbf{u}}=0,\qquad\nabla\tilde{p}=\nabla^{2}\tilde{\mathbf{u}},\qquad(\mathbf{r}\in\Omega_{1}),$$

(3)

and on the free surface satisfy the standard kinematic and tangential-stress boundary conditions,

$$\tilde{\mathbf{u}}\cdot\tilde{\mathbf{n}}=0,\qquad\tilde{\mathbf{n}}\cdot[\nabla\tilde{\mathbf{u}}+(\nabla\tilde{\mathbf{u}})^{T}]\cdot(\mathbf{I}-\tilde{\mathbf{n}}\tilde{\mathbf{n}})=0,\qquad(\mathbf{r}\in S_{f}),$$

(4)

where $\mathbf{I}$ is the metric tensor.

In the porous medium, the flow is driven by the capillary pressure in the menisci that collectively form the wetting front, and, for the problem in question, the characteristic pressure and velocity are $2\sigma/a$ and $U$, respectively. Using the notation $p$ and $\mathbf{u}$ for the dimensionless pressure and velocity in $\Omega_{2}$, in a frame moving with the contact line we have an equation of motion in the form of Darcy’s law,

$$\mathbf{u}-\hat{\mathbf{U}}=-K\nabla p,\qquad(\mathbf{r}\in\Omega_{2}),$$

(5)

where $K=2\sigma\kappa/(\mu aLU)$ is the nondimensionalized permittivity of the porous matrix ($\kappa$ is the actual permittivity) and $\hat{\mathbf{U}}=\mathbf{U}/U$ is a unit vector directed along the velocity of the porous substrate. Given that in porous media $\kappa\propto a^{2}$ (Probstein, 1989), we have that $K=O(\epsilon/Ca)$ and, for the problem to be nontrivial, it is assumed to be finite in the limits (1), (2). In the present context, it is convenient to define $L$ by setting $K=1$.

Since the porosity $\phi$ in a homogeneous matrix is constant, the mass balance equation has the standard form $\nabla\cdot\mathbf{u}=0$, so that, after substituting (5) into it, one arrives at Laplace’s equation for $p$ in $\Omega_{2}$:

$$\nabla^{2}p=0,\qquad(\mathbf{r}\in\Omega_{2}).$$

(6)

On the wetting front, one has the kinematic condition that the front propagates with the velocity of the fluid, i.e. in the reference frame described above where the process is steady,

$$\mathbf{u}\cdot\mathbf{n}=0,\qquad(\mathbf{r}\in S_{w}),$$

(7)

and needs to specify the dynamic condition on the pressure $p$. To formulate this condition in a general case, one has to consider the modes of motion the menisci go through on the pore scale as the macroscopic wetting front propagates (Shikhmurzaev & Sprittles, 2012; anomalous-12). Here we will be considering the simplest case involving the main, ‘wetting’, mode. In this case, for an ‘effective’ pore with a circular cross-section the dimensionless pressure at the front is given by

$$p=-\cos\theta_{d},$$

(8)

where $\theta_{d}$ is the contact angle formed on the pore scale by the representative meniscus with the pore wall (Fig. 1). Unlike the simplified approach pioneered by Washburn (1921), where it is assumed that $\theta_{d}$ is constant and equal to the prescribed static contact angle $\theta_{s}$, we need to take into account that, as demonstrated by numerous experiment (e.g. see Ch. 3 of Shikhmurzaev, 2007, for a review), $\theta_{d}$ depends on the wetting speed, i.e. to consider the dependence

$$(\mathbf{u}-\hat{\mathbf{U}})\cdot\mathbf{n}=U^{*}_{cl}f(\theta_{d}),\qquad(\mathbf{r}\in S_{w}),$$

(9)

where $U^{*}_{cl}$ is the appropriate velocity scale depending on the material parameters of the system (as all velocities above, we have it nondimensionalized using $U$). The function $f(\theta_{d})$ has to be determined theoretically or empirically. For example, one can apply $f(\theta_{d})$ derived using the theory of dynamic wetting as a process of interface formation (Shikhmurzaev, 2007), which has been shown to reliably describe experimental data, though, in the present context, any function $f(\theta_{d})$ representing the experimentally observed dependence of the form (9) could be used, such as, for instance, the one that comes from the molecular-kinetic theory of wetting (Blake & Haynes, 1969). For more information about the dynamic wetting modelling we refer the reader to a recent review (Shikhmurzaev, 2011).

From the theory of flows with forming interfaces one has

$$f(\theta_{d})=\left(\frac{(1+(1-\rho^{s}_{1e})\cos\theta_{s})(\cos\theta_{s}-\cos\theta_{d})^{2}}{4(\cos\theta_{s}+B)(\cos\theta_{d}+B)}\right)^{1/2},\qquad U^{*}_{cl}\equiv\frac{U_{cl}}{U}=\frac{1}{U}\left(\frac{\gamma\rho^{s}_{(0)}(1+4\alpha\beta)}{\tau\beta}\right)^{1/2}$$

(10)

where $B=(1-\rho^{s}_{1e})^{-1}(1+\rho^{s}_{1e}u_{0}(\theta_{d}))$,

$$u_{0}(\theta_{d})=\frac{\sin\theta_{d}-\theta_{d}\cos\theta_{d}}{\sin\theta_{d}\cos\theta_{d}-\theta_{d}},$$

and $\rho^{s}_{(0)}$, $\rho^{s}_{1e}$, $\alpha$, $\beta$, $\gamma$, $\tau$ are material constants characterizing the contacting media. Their values for some systems can be found elsewhere (Shikhmurzaev, 2007; Blake & Shikhmurzaev, 2002).

So far, the flow in the pure fluid and in the porous medium were considered separately, and to link them one has to specify three boundary conditions at $S_{0}$. One condition that we obviously have on this surface is the continuity of mass flux:

$$(\mathbf{\tilde{u}}-\phi\mathbf{u})\cdot\mathbf{n}_{0}=0,\qquad(\mathbf{r}\in S_{0}).$$

(11)

For the velocity components parallel to $S_{0}$, i.e. for $\tilde{\mathbf{u}}_{\parallel}\equiv\tilde{\mathbf{u}}\cdot(\mathbf{I}-\mathbf{n}_{0}\mathbf{n}_{0})$ and $\mathbf{u}_{\parallel}\equiv\mathbf{u}\cdot(\mathbf{I}-\mathbf{n}_{0}\mathbf{n}_{0})$, a number of boundary conditions have been discussed in the literature (e.g. Saffman, 1971; Jones, 1973; Murdoch & Soliman, 1999; Nield, 2009; Auriault, 2010), following the experiments reported by Beavers & Joseph (1967) and an empirical condition these authors proposed. As noted, for example, by Auriault (2010), these conditions are aiming at capturing the effects of order $O(\epsilon)$, i.e. go beyond the classical, i.e. 0th-order, approximation in the continuum limit (1). In the 0th approximation, all these conditions reduce to no-slip for the pure fluid,

$$\tilde{\mathbf{u}}_{\parallel}=\hat{\mathbf{U}}_{\parallel},\qquad(\mathbf{r}\in S_{0}),$$

(12)

and it is this condition that we will be using here.

The condition of continuity of pressure on $S_{0}$ in the dimensionless form yields that

$$p=\tilde{p}\,\frac{\epsilon\,Ca}{2},$$

and hence, to leading order in the limits (1) and (2), one has

$$p=0,\qquad(\mathbf{r}\in S_{0}).$$

(13)

Importantly, unlike the case of an impermeable solid substrate, where one has to specify the dynamic contact angle (as we need to specify $\theta_{d}$ on the pore scale), here, for the effective contact angle $\theta_{D}$, we will require only that the flow parameters in the porous medium remain regular at the contact line and that the solution in the pure fluid exists.

Consider the asymptotic behaviour of the solution to (3)–(13) in the case of a two-dimensional flow as the distance to the contact line $r\to 0$. In the polar coordinates $(r,\theta)$ shown in Fig. 1, to leading order, for the pressure $p$ in the porous medium one has

$$\frac{\partial^{2}p}{\partial r^{2}}+\frac{1}{r}\frac{\partial p}{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}p}{\partial\theta^{2}}=0,\qquad(-\theta_{w}<\theta<0),$$

$$p(r,0)=0,\qquad\frac{\partial p}{\partial\theta}(r,-\theta_{w})=r\sin\theta_{w},$$

where the last condition follows from (7) and (5). The separable solution to this problem is obviously given by

$$p=r\tan\theta_{w}\sin\theta,$$

(14)

so that, using the dynamic boundary condition (8), we obtain that

$$\theta_{d}=\arccos(-r\tan\theta_{w}\sin\theta_{w})\to\frac{\pi}{2},\qquad\hbox{as }r\to 0.$$

Then, from (9), where now $(\mathbf{u}-\hat{\mathbf{U}})\cdot\mathbf{n}=\sin\theta_{w}$, we have an equation determining $\theta_{w}$:

$$\sin\theta_{w}=U^{*}_{cl}f(\pi/2).$$

(15)

In order to consider the flow in the pure fluid, we introduce a stream function

$$\tilde{u}_{r}=\frac{1}{r}\frac{\partial\tilde{\psi}}{\partial\theta},\qquad\tilde{u}_{\theta}=-\frac{\partial\tilde{\psi}}{\partial r},$$

(16)

so that for the leading-order term $\tilde{\psi}_{1}$ of the asymptotic expansion of $\tilde{\psi}$ as $r\to 0$ we have a biharmonic equation

$$\left(\frac{\partial^{2}}{\partial r^{2}}+\frac{1}{r}\frac{\partial}{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\theta^{2}}\right)^{2}\tilde{\psi}_{1}=0,\qquad(0<r,0<\theta<\theta_{D}),$$

(17)

together with conditions (4), i.e.

$$\tilde{\psi}_{1}(r,\theta_{D})=0,\qquad\frac{\partial^{2}\tilde{\psi}_{1}}{\partial\theta^{2}}(r,\theta_{D})=0,$$

(18)

on the free surface, and conditions (11) and (12), i.e.

$$\tilde{\psi}_{1}(r,0)=r\phi\tan\theta_{w},\qquad\frac{\partial\tilde{\psi}_{1}}{\partial\theta}(r,0)=r,$$

(19)

on the surface of the solid substrate. In writing down the first of conditions (19) we made use of (14) and integrated along $S_{0}$.

The separable solution to the problem (17)–(19) has the form

$$\tilde{\psi}_{1}=r(A_{1}\sin\theta+A_{2}\theta\sin\theta+A_{3}\cos\theta+A_{4}\theta\cos\theta),$$

(20)

where the constants $A_{1},\dots,A_{4}$ are given by

$$A_{1}=-\frac{\theta_{D}}{\sin\theta_{D}\cos\theta_{D}-\theta_{D}}-A_{3}\frac{\cos^{2}\theta_{D}}{\sin\theta_{D}\cos\theta_{D}-\theta_{D}},$$

$$A_{2}=\frac{\sin^{2}\theta_{D}}{\sin\theta_{D}\cos\theta_{D}-\theta_{D}}+A_{3}\frac{\sin\theta_{D}\cos\theta_{D}}{\sin\theta_{D}\cos\theta_{D}-\theta_{D}},$$

(21)

$$A_{3}=\phi\tan\theta_{w},\qquad A_{4}=A_{2}\cot\theta_{D}.$$

(22)

Using (16), we can write down the radial projection of the second equation in (3) in the form

$$\frac{\partial\tilde{p}}{\partial r}=\left(\frac{1}{r}\frac{\partial^{3}}{\partial r^{2}\partial\theta}+\frac{1}{r^{3}}\frac{\partial^{3}}{\partial\theta^{3}}+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial r\partial\theta}\right)\tilde{\psi}_{1},$$

and, after substituting the solution (20), arrive at

$$\frac{\partial\tilde{p}}{\partial r}=-\frac{2}{r^{2}}(\sin\theta+\cot\theta_{D}\cos\theta)A_{2}.$$

Thus, the leading term in the coordinate expansion of the stream function will not give rise to a nonintegrable singularity of pressure if and only if $A_{2}=0$, or, using (21) and (22),

$$\tan\theta_{D}=-\phi\tan\theta_{w}.$$

Given the expression (15) for $\sin\theta_{w}$ and introducing the dimensionless contact-line speed $U_{*}\equiv 1/(f(\pi/2)U^{*}_{cl})=U/(f(\pi/2)U_{cl})$, we can also write this equation down as

$$\theta_{D}=\pi+\arctan\left(-\frac{\phi}{\sqrt{U^{2}_{*}-1}}\right).$$

(23)

This equation specifies $\theta_{D}$ in terms of the speed $U$ of the contact line with respect to the substrate, the porosity of the matrix $\phi$ and the material properties of the contacting media accumulated in $U_{cl}$ and $f(\pi/2)$. The velocity-dependence of $\theta_{D}$ for different porosities is shown in Fig. 2. As one can see, this dependence is much steeper than the velocity-dependence of the dynamic contact angle for the dynamic wetting of an impermeable substrate (e.g. see Ch. 3 of Shikhmurzaev, 2007, for a review). In particular, $d\theta_{D}/dU_{*}\to+\infty$ as $U_{*}\to 1+$.

The dynamic contact angle given by (23) ensures integrability of the normal stress on the free surface and hence the existence of a solution. The next term in the asymptotic expansion of $\tilde{\psi}$ has the form $\tilde{\psi}_{2}=r^{2}F(\theta)$ and therefore can give rise, at most, to a logarithmic (i.e. integrable) singularity of $\tilde{p}$ at the contact line, which does not affect the existence of the solution.

The obtained result has a clear physical meaning. In the case we are considering, the dynamics of imbibition determines the normal to the substrate component of the pure fluid’s velocity independently of the pure fluid’s bulk flow, and, given that the tangential component of the fluid’s velocity on the surface of the substrate is prescribed, as it satisfies the no-slip condition with a known speed of the substrate, we have a moving contact-line problem where on the solid surface both components of velocity are set. In this situation, in the pure fluid the solution exists only if from the family of stream functions described by (20) we choose the one that corresponds to a uniform flow ($A_{2}=0$, $A_{4}=0$) with the projections of velocity on the normal and tangential to the substrate directions equal to the speed of imbibition and the speed of the solid substrate, respectively.

The solution (23) for a steady spreading of a fluid over an unsaturated substrate ceases to exist when the dimensionless velocity $U_{*}$ becomes equal to $1$ and both $\theta_{D}$ and $\theta_{w}$ reach $\pi/2$. For $U_{*}<1$ we have a different regime with the two contact lines CL1 and CL2 no longer coinciding as the wetting front moves ahead (dashed line in Fig. 1), and the pure fluid finds itself on a saturated substrate. Then, in the vicinity of CL1 the pressure $p$ satisfies Laplace’s equation $\nabla^{2}p=0$ together with the boundary conditions

$$p(r,0)=0,\qquad\frac{\partial p}{\partial\theta}(r,-\pi)=0,$$

and the local solution with the regular pressure gradient has the form

$$p=\sum\limits_{n=1}^{\infty}C_{n}r^{1/2+n}\sin\left(\frac{1}{2}+n\right)\theta,$$

(24)

where $C_{n}$, $n=1,2,\dots$ are constants. Then, according to (5) and (11), one has $\tilde{u}_{\theta}(r,0)\propto r^{1/2}\to 0$ as $r\to 0$. As a result, if $\tilde{u}_{r}(r,0)=O(1)$ as $r\to 0$, then, to leading order as $r\to 0$, we will have the ‘classical’ moving contact-line problem, with no imbibition and the no-slip boundary condition on the solid. This problem, as is well known (see, e.g., Ch. 3 of Shikhmurzaev, 2007), has no solution. The only way out of this situation is to conclude that, as CL1 and CL2 separate, the contact line CL1 stops moving. Then, for the pure fluid one simply has a static contact line with the imbibition velocity $\propto r^{1/2}$ near it and $\theta_{D}=\pi/2$. It is easy to verify that the solution to this problem exists. As CL1 and CL2 separate and CL2 moves ahead, the imbibition process near CL1 slows down, leading to $C_{1}=0$ in (24) and hence $\tilde{u}_{\theta}\propto r^{3/2}$, so that $\theta_{D}$ no longer has to be equal to $\pi/2$ and can go down as the imbibition continues.

The described scenario is in agreement with the available experimental data. In experiments on the spreading drops, it has been observed that the regime where the drop base and the saturated area underneath it expand together is followed by the regime where the two contact lines, CL1 and CL2, separate, as drop’s base stops expanding whereas the saturated area continues to grow (Clarke et al., 2002; Starov et al., 2003; Keshav & Basu, 2007; Markicevic et al., 2010). This is usually attributed to the ‘competition’ between spreading and imbibition, and the above analysis shows what this actually means.

Markicevic et al. (2010) report that in their experiment the spreading droplet maintained the shape of a spherical cap throughout the process and that its base stopped expanding when “the droplet is a half of sphere”. In other words, it stopped expanding when the dynamic contact angle $\theta_{D}$ became equal to $90^{\circ}$. This is exactly what follows from the result of our analysis.

The dependence of the effective contact angle $\theta_{D}$ on the contact-line speed for $U_{*}>1$ given by (23) and illustrated in Fig. 2 awaits its experimental verification. The issue here is that, to extract this dependence from experiments with unsteady flows, one has to deal with very short time intervals ($\ll 40$ ms according to Markicevic et al. (2010)) with a well-controlled spatial resolution and conditions of the Darcy-scale description satisfied. The problem is made more complicated by the fact that, as one can see in Fig. 2, the velocity-dependence of $\theta_{D}$ is very steep, which brings in additional conditions on the temporal resolution of the experiments. On the other hand, experiments with steady flows of the kind commonly used to study dynamic wetting of impermeable solid surface require sizeable substrates with well-reproducible properties. To date, no systematic measurements of $\theta_{D}$ with controlled spatial resolution have been reported.

As we have shown, the requirement of the absence of a nonintegrable singularity of pressure in the pure fluid at the contact line, which is equivalent to the requirement of the existence of a solution, uniquely specifies the velocity-dependence of the ‘effective’ dynamic contact angle $\theta_{D}$ formed by the free surface and the ‘effective’ surface of the porous substrate. In practical computations, it may be convenient to set (23) as a boundary condition while considering the normal-stress condition on the free surface as the equation determining the free surface shape. One implication of the obtained result is that, if one specifies $\theta_{D}$ as an additional boundary condition different from (23), then, for a solution to exist, it will become necessary to use a slip boundary condition instead of the no-slip condition (12) employed here (Davis & Hocking, 2000). As mentioned earlier, this would mean bringing in effects of $O(\epsilon)$ into the essentially classical, i.e. $O(1)$ as $\epsilon\to 0$, formulation. If the velocity-dependence of $\theta_{D}$ different from (23) is imposed together with the no-slip condition (12) for the pure fluid (Alleborn & Raszillier, 2004; Reis Jr. et al., 2004), then there will be no solution to the problem, although this fact can be hidden behind, and as a result masked by, simplifying assumptions made in the process of finding the solution (e.g. lubrication approximation) and the numerical implementation.

From the theoretical viewpoint, it would be interesting to consider the essentially unsteady process of separation of CL1 and CL2, which marks the transition between the two regimes described above. As with every finite-time transition, this is a challenging problem deserving a detailed investigation.