2.1. Full hydrodynamic model

Motivated by the system used in Vandre et al. (Reference Vandre, Carvalho and Kumar2012), which is representative of an experimental system, we consider two-dimensional flow between two parallel plates, as shown in figure 1. In the following discussion we denote dimensional quantities using an asterisk. Two fluids of viscosity $\mu _{1,2}^*$, and density $\rho _{1,2}^*$, fill the channel bounded by two rigid plates which are separated by a fixed distance $H^*$; subscript 1 indicates the upper fluid (the fluid phase) and subscript 2 indicates the lower fluid (the liquid phase). In our system the left plate moves with constant speed in the $y$ direction $U^*$ and the right plate is stationary. For a receding contact line $U^*>0$ and an advancing contact line $U^*<0$. The fluid flow of each phase is governed by the two-dimensional Navier–Stokes equations. All speeds, lengths, pressures and times are scaled by $U^*$, $H^*$, $\mu _2^*U^*/H^*$ and $H^*/U^*$ respectively. Finally, the viscosity ratio, denoted $\chi$, is defined with respect to the liquid phase, i.e.

(2.1)\begin{equation} \chi = \mu_{1}^*/\mu_{2}^*, \end{equation}

and we assume that the upper fluid is less viscous, i.e. $\chi < 1$.

As in previous studies (Sprittles & Shikhmurzaev Reference Sprittles and Shikhmurzaev2011a,Reference Sprittles and Shikhmurzaevb; Vandre et al. Reference Vandre, Carvalho and Kumar2012; Sprittles & Shikhmurzaev Reference Sprittles and Shikhmurzaev2013; Vandre et al. Reference Vandre, Carvalho and Kumar2013; Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016a,Reference Liu, Vandre, Carvalho and Kumarb, Reference Liu, Carvalho and Kumar2017, Reference Liu, Carvalho and Kumar2019) we apply the Stokes-flow approximation so that the Reynolds number, $Re = U^*H^*\rho _2^*/\mu _2^*$, is negligible and assumed zero; results in Vandre et al. (Reference Vandre, Carvalho and Kumar2013) show $Re$ can have an influence at sufficiently high values but it does not qualitatively alter the conclusions. We assume that gravitational effects are negligible throughout. The non-dimensional computational domain is shown in figure 3(c). The coordinate system is centred on the contact point between the two fluids and the left (moving) plate. The boundary corresponding to the left plate is denoted $\varGamma _1$, the right plate $\varGamma _2$, the bottom $\varGamma _3$, the free surface $\varGamma _4$ and the top $\varGamma _5$. The fluid and liquid domains are denoted by $\varOmega _1$ and $\varOmega _2$ respectively. The Stokes-flow equations for the fluid velocity, $\boldsymbol {u}_i = (u_i,v_i)$, and pressure, $p_i$, in each phase can be written as

(2.2)\begin{align} \chi \nabla^2\boldsymbol{u}_1 &= \boldsymbol{\nabla} p_1,\quad\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_1 = 0,\quad \boldsymbol{x}\in \varOmega_1\quad\text{(fluid)}, \end{align}

(2.3)\begin{align} \nabla^2\boldsymbol{u}_2 &= \boldsymbol{\nabla} p_2,\quad\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_2 = 0,\quad \boldsymbol{x}\in \varOmega_2\quad\text{(liquid)}. \end{align}

Figure 3. The computational domain for the hybrid and full models for an advancing contact line. The boundaries are denoted by $\varGamma _i$ (labelled in c) with the origin centred on the contact line. (a) A typical streamline pattern for a steady solution in the hybrid model. (b) The computational domain for the hybrid model. In this model we solve for the velocity and pressure in the liquid domain, but only solve for the pressure of the fluid on the interface boundary, $\varGamma _4$. (c) The computational domain for the full model where the velocity and pressure fields are also solved in the fluid domain. (d) The vertical component of velocity along the moving plate, i.e. $\varGamma _1$. (e) The enlargement near the contact line shows the mesh refinement required to ensure that the flow field is sufficiently resolved.

On the left (moving) and right (stationary) plates, $\varGamma _1$ and $\varGamma _2$ respectively, we implement a Navier slip condition written as

(2.4)\begin{align} \lambda ( \boldsymbol{\tau}_i \boldsymbol{\cdot} \boldsymbol{n}_{{p}})\boldsymbol{\cdot}\boldsymbol{t}_{{p}} &= (\boldsymbol{u}_{i} - \boldsymbol{U})\boldsymbol{\cdot} \boldsymbol{t}_{{p}},\quad\boldsymbol{x}\in\varGamma_{1},\quad i = 1,2, \end{align}

(2.5)\begin{align} \lambda ( \boldsymbol{\tau}_i \boldsymbol{\cdot} \boldsymbol{n}_{{p}})\boldsymbol{\cdot}\boldsymbol{t}_{{p}} &= \boldsymbol{u}_{i} \boldsymbol{\cdot} \boldsymbol{t}_{{p}},\quad\boldsymbol{x}\in\varGamma_{2},\quad i = 1,2, \end{align}

where $\boldsymbol {n}_{{p}}$ and $\boldsymbol {t}_{{p}}$ are the vectors normal and tangential to each plate respectively, $\boldsymbol {U} = (0,U)$ is the non-dimensional speed of the (moving) left plate (where $U$, in our dimensionless system, is $\pm 1$, corresponding to the receding (+)/advancing (−) problem) and $\lambda$ is the non-dimensional slip length which, for simplicity, we assume to be the same in each phase (see Sprittles (Reference Sprittles2017) for potential extensions). We could choose different conditions that regularise the singularity at the contact line (Shikhmurzaev Reference Shikhmurzaev2006), but, assuming the actual contact angle is unchanged, the details of the solution (i.e. the value of $Ca_{{crit}}$) are more sensitive to values of the slip-length parameter that arises in these models (in our case $\lambda$) than the actual form of the model (Dussan Reference Dussan1976). Thus, we would expect the results we obtain to be qualitatively similar to those obtained for a different slip model, although such slip models are easy to implement, if required. We choose a moderate value of $\lambda = 0.1$, and although quantitative details, i.e. the values of $Ca_{{crit}}$ and other solution measures, will differ as we vary $\lambda$, we find that the transient behaviour and solution structure, as we describe later in the article, are qualitatively the same (cf. figure 7).

We choose to implement a Navier slip condition on the stationary plate for consistency and to ensure that the contact point on the stationary plate is allowed to move. However, we could fix $\boldsymbol {u}_i = 0$ on $\varGamma _2$ and get similar results (see e.g. Vandre et al. Reference Vandre, Carvalho and Kumar2013; Liu et al. Reference Liu, Carvalho and Kumar2017) which corresponds to a pinned contact line. We also implement a no-penetration condition on each plate, i.e.

(2.6)\begin{equation} \boldsymbol{u}_i\boldsymbol{\cdot} \boldsymbol{n}_{{p}} = 0,\quad\boldsymbol{x}\in\varGamma_1\cup\varGamma_2,\quad i = 1,2. \end{equation}

The stress tensor in each phase $\boldsymbol {\tau }_i$ is defined as

(2.7)\begin{equation} \boldsymbol{\tau}_i ={-}p_i\boldsymbol{I} + \delta_i(\boldsymbol{\nabla} \boldsymbol{u}_i + (\boldsymbol{\nabla}\boldsymbol{u}_i)^{\rm T}), \end{equation}

where $\boldsymbol {I}$ is the identity matrix and $\delta _1 = \chi, \delta _2 = 1$. We denote the unknown position of the interface, $\varGamma _4$, as $\boldsymbol {r} = (x_s,y_s)$ (see figure 3b) and assume a constant surface tension, $\gamma ^*$, so that the dynamic boundary condition can be written as

(2.8)\begin{equation} \boldsymbol{\tau}_1\boldsymbol{\cdot}\boldsymbol{n} - \boldsymbol{\tau}_2\boldsymbol{\cdot}\boldsymbol{n} = \frac{1}{Ca}\kappa \boldsymbol{n},\quad\boldsymbol{x}\in\varGamma_{4}, \end{equation}

where $\boldsymbol {n}$ is the normal of the interface pointing towards the fluid phase (see figure 3b), $\kappa = \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {n}$ is the curvature of the interface and $Ca = \mu _2^*|U|/\gamma ^*$ is the capillary number. In addition to (2.8), we impose a kinematic condition on $\varGamma _4$, written as

(2.9)\begin{equation} \frac{\partial \boldsymbol{r}}{\partial t}\boldsymbol{\cdot}\boldsymbol{n} = \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n}, \quad\boldsymbol{x}\in\varGamma_{4}. \end{equation}

We have to specify the angle that the interface makes on the left and right plates. These angles can be allowed to vary with the capillary number, slip length or other quantities but we choose the simplest approach and take these to be constant values, i.e.

(2.10)$$\begin{gather} \theta = \theta_1,\quad \text{on}\ \varGamma_1\cap\varGamma_4, \end{gather}$$

(2.11)$$\begin{gather}\theta = \theta_2,\quad \text{on}\ \varGamma_2\cap\varGamma_4. \end{gather}$$

It is straightforward to replace these conditions with equations involving $Ca$ and other quantities, but this is not the focus of this article.

Finally, we implement fully developed flow conditions on the inflow and outflow boundaries:

(2.12)\begin{equation} \boldsymbol{u}_i\boldsymbol{\cdot}\boldsymbol{t}_{{inflow}} = 0,\,\boldsymbol{x}\in\varGamma_{3}\cup\varGamma_{5}, \end{equation}

where $\boldsymbol {t}_{{inflow}} = (1,0)^{\rm T}$, alongside a pressure drop across the domain so that

(2.13)$$\begin{gather} p_1 = 0,\,\boldsymbol{x}\in\varGamma_{5}, \end{gather}$$

(2.14)$$\begin{gather}p_2 = p_{{out}},\,\boldsymbol{x}\in\varGamma_{3}. \end{gather}$$

The full hydrodynamic system is defined in (2.2)–(2.14) with the following infinite-dimensional state vector (denoted $\boldsymbol {w}$) of unknowns:

(2.15)\begin{equation} \boldsymbol{w} = [\boldsymbol{u}_1,\boldsymbol{u}_2,p_1,p_2,\boldsymbol{r}]^{\rm T}. \end{equation}

It is worth noting that we model the effect of varying the speed of the plate by varying $Ca$ and that the non-dimensional slip length, $\lambda$, can be varied to investigate changes in physical channel width. Finally, as our primary interest is in understanding the transient behaviour, for simplicity, we set $\theta _1 = \theta _2 = {\rm \pi}/2$ in all simulations. Therefore we have a set of control parameters,

(2.16)\begin{equation} Ca,\lambda,\chi,p_{{out}}, \end{equation}

that need to be specified in order to solve (2.2)–(2.14).

2.2. Hybrid model

The computational cost of the full model can be significantly reduced by solving the thin-film equations where they are valid (Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Jacqmin Reference Jacqmin2004; Sbragaglia, Sugiyama & Biferale Reference Sbragaglia, Sugiyama and Biferale2008), leading to a hybrid model (see Stay & Barocas Reference Stay and Barocas2003; Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016a,Reference Liu, Vandre, Carvalho and Kumarb, Reference Liu, Carvalho and Kumar2017, Reference Liu, Carvalho and Kumar2019; Charitatos et al. Reference Charitatos, Suszynski, Carvalho and Kumar2020) which approximately halves the complexity of the problem, as unknowns in the fluid phase are only computed on the interface. The difference of our approach from previous implementations is that our hybrid model takes into account time dependence so that stability can be probed and IVP calculations can be performed. A key assumption is that a typical horizontal distance in the fluid phase, when air entrainment occurs, is small when compared to the vertical height of the meniscus (i.e. $\hat {h}\ll Y$ in figure 1) so that the flow is approximately parallel, i.e. the cross-stream component of $\boldsymbol {u}_1$ is small. The full derivation is discussed in Appendix A and the computational domain is shown in figure 3(b).

The effect of this reduction in the fluid phase is to replace a full two-dimensional description, given in (2.2), by a one-dimensional equation for the fluid pressure, $p_1$, on the interface only. This equation can be stated as

(2.17a,b)\begin{equation} \frac{\partial \boldsymbol{r}}{\partial t}\boldsymbol{\cdot}\boldsymbol{n} \pm \frac{1}{\chi}\frac{\partial Q_1}{\partial s} = 0,\quad Q_1 = \frac{1}{6}\frac{\partial p_1}{\partial s}h^3 + \frac{1}{2}Ah^2 + Bh, \end{equation}

where $h$ is the horizontal distance (i.e. $x_s$) from the left plate to the interface (see figure 3), $Q_1$ is the flux and the constants $A$ and $B$ are functions of $\chi,\lambda$ and $\boldsymbol {u}_2$ and are given in Appendix A. The $\pm$ sign is used for the advancing ($+$)/receding ($-$) contact line. The fluid phase is coupled to the liquid phase through the applied traction given in (A6). We now have a system of partial differential equations (PDEs) described by (2.3)–(2.14) and (2.17a,b) with the infinite-dimensional state vector of unknowns:

(2.18)\begin{equation} \boldsymbol{w} = [\boldsymbol{u}_2,p_2,p_1(\boldsymbol{r}),\boldsymbol{r}]^{\rm T}. \end{equation}

The hybrid model has been extensively validated for steady advancing contact-line problems against the full system and experiment (Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016b, Reference Liu, Carvalho and Kumar2019). We validate the time-dependent hybrid model by comparing to the full hydrodynamic model in Appendix A. Finally we note that this approach is strictly only valid for the advancing contact-line problem but, as shown later, the receding contact-line problem is effectively a one-phase problem (cf. figure 7), and implementing the hybrid model for a receding contact line does not significantly change the value of $Ca_{{crit}}$ (Appendix A).

2.3. System parameters and integral measures

We now describe additional system parameters and measures that will be useful in computing and describing the steady and time-dependent solutions. The pressure at the outflow boundary, $p_{{out}}$, is determined implicitly by an integral volume constraint acting on the liquid phase, i.e.

(2.19)\begin{equation} \iint_{{\varOmega_2}}\,{\rm d}V = V, \end{equation}

where $V$ is the volume of the liquid domain (corresponding to the area of the computational domain). In our numerical calculations the positions of the contact points on the moving and stationary plates are both allowed to move so that (2.19) can be satisfied. For ease of presentation, we post-process and rescale the solution so that the origin is always at the contact point of the moving plate. In the calculations that follow we choose a value of $V$ which is large enough for fully developed flow to occur near the outflow boundary, $\varGamma _1$. After careful experimentation we find that the solutions are independent of values of $V\geq 5$ that we choose.

For a fixed set of parameter values, as defined in (2.16), we calculate steady solution curves by setting the time derivatives in the governing equations to zero and then solving the resulting steady system. As we vary $Ca$, and then subsequently calculate a solution, a solution curve will be traced and a fold bifurcation will occur at the critical value of the capillary number, denoted $Ca_{{crit}}$. Whilst it is possible to trace a solution branch around the fold numerically by a pseudo-arclength continuation method (see e.g. Doedel Reference Doedel2007), we implement an alternative, bespoke approach. We expect the interface length, $L$, to increase monotonically as the curve is traced out around the fold and therefore it is a suitable candidate for a continuation parameter that allows us to calculate solutions smoothly around the fold. To achieve this, we let $Ca$ become an unknown parameter that is determined implicitly by setting the total length of the interface, i.e.

(2.20)\begin{equation} \int_{\varGamma_{4}}\,{\rm d}s = L. \end{equation}

This approach enables us to trace solution curves around the fold by incrementally increasing $L$ and solving the system of equations, with $Ca$ effectively determined by (2.20). We also emphasise that (2.19) and (2.20) are only implemented in steady calculations; the former constraint is unnecessary in time-dependent calculations as the second equation in (2.3) ensures volume is conserved, whilst the latter constraint is used as a means of tracing the solution curve.

Finally, when describing the steady solutions and time-dependent solutions we use the meniscus rise (more specifically, the vertical distance between the two contact lines) defined as

(2.21)\begin{equation} Y(t) = |y_s(s=L) - y_s(s=0)| \end{equation}

as a convenient solution measure (as previously considered, for example, in Kamal et al. (Reference Kamal, Sprittles, Snoeijer and Eggers2019); see figure 3 (in this article) for reference).

2.4. Numerical method

The governing equations are solved using the finite-element method from within the open-source oomph-lib object-orientated multi-physics library (Heil & Hazel Reference Heil and Hazel2006). The structure and implementation of our equations follow that of Sprittles & Shikhmurzaev (Reference Sprittles and Shikhmurzaev2013). Following multiplication of the equations by a test function, $\psi$, and then an integration over the domain, the boundary integrals that result from integration by parts require the traction to be specified on each of the boundaries. The dynamic condition, (2.8), and the Navier slip condition, (2.5), therefore can be implemented as a natural condition by these boundary integrals.

Special care has to be taken at the contact point. In other studies (Vandre et al. Reference Vandre, Carvalho and Kumar2012, Reference Vandre, Carvalho and Kumar2013; Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016a,Reference Liu, Vandre, Carvalho and Kumarb, Reference Liu, Carvalho and Kumar2017, Reference Liu, Carvalho and Kumar2019) the contact angle is imposed as an essential boundary condition at the expense of solving a component of the momentum equations at the contact point. We adopt the approach of Sprittles & Shikhmurzaev (Reference Sprittles and Shikhmurzaev2013) and impose the contact angle as a natural boundary condition on both the intersection of the free surface with the left plate ($\varGamma _1$) and the symmetry plate ($\varGamma _2$). We therefore introduce a field of Lagrange multiplier unknowns on $\varGamma _1$ and $\varGamma _2$ which are determined from the weak form of the no-penetration condition, (2.6). We refer the reader to Sprittles & Shikhmurzaev (Reference Sprittles and Shikhmurzaev2013) for a detailed description of this implementation (we adopt approach (B) in their nomenclature).

As is standard, the fluid velocities are interpolated using bi-quadratic shape functions and the pressure using linear continuous shape functions with Taylor–Hood triangular elements. We choose to mesh the liquid domain using an unstructured triangular grid (see figure 3b,c). The mesh is considered to be a fictitious pseudo-solid with the position of the nodes coming as part of the solution. The weak form of the kinematic condition, (2.9), is imposed as an essential condition and determines a field of Lagrange multipliers (not to be confused with the Lagrange multipliers in the previous paragraph) that act on the solid deformation equations which in turn determines the shape of the unknown interface, $\boldsymbol {r}$; see Sackinger, Schunk & Rao (Reference Sackinger, Schunk and Rao1996) for more details. We note that this approach results in a large system of equations which is disadvantageous, but it allows for the interface to become highly deformed, as well as naturally handling time-dependent flow (where the domain could significantly change shape; cf. figure 1), which is difficult to achieve if the mesh is structured.

To solve the hybrid equation, (2.17a,b), it is convenient to introduce two fields of unknowns on the fluid interface, the pressure $p_1$ and flux $Q_1$, interpolated using quadratic shape functions. We solve two equations in their weak form:

(2.22a,b)\begin{equation} \int_{S}(Q_1 - Q_{C})\psi\,{\rm d}S = 0,\quad Q_{C} = \frac{1}{\chi}\left(Ah + \frac{1}{2}Bh^2 + \frac{1}{6}h^3\frac{\partial p}{\partial s}\right) \end{equation}

and

(2.23)\begin{equation} \int_{S} \left(\frac{\partial \boldsymbol{r}}{\partial t}\boldsymbol{\cdot}\boldsymbol{n} \pm \frac{\partial Q_1}{\partial s}\right)\psi\,{\rm d}s = 0. \end{equation}

Equation (2.22a,b) projects the flux from the lubrication equation onto the finite-element space and then (2.23) ensures mass is conserved in the fluid phase.

The resulting discretised equations are solved with Newton's method using the SuperLu numerical algebra package (Li Reference Li2005). For time-dependent calculations the solution is updated in time using a backwards-difference second-order Euler method (BDF2). A typical streamline pattern is shown in figure 3(a).

Around the contact line the interface becomes highly deformed due to viscous bending and the pressure and velocity gradients are large (see figure 3d,e). In steady calculations, as $Ca\to Ca_{{crit}}$, we expect the number of elements required in the vicinity of the contact line to increase to ensure a smooth converged solution. We re-mesh the domain according to a ZZ error estimator (Zienkiewicz & Zhu Reference Zienkiewicz and Zhu1992), which measures the discontinuity of strain-rate gradients between adjacent elements and interprets this as a measure for the local error. Typically we set a minimum error as $10^{-6}$ and a maximum error $10^{-3}$, so that elements with error above this range get refined and those with error below this range get unrefined. We allow element sizes from $10^{-12}$ to $10^{-2}$ to accommodate these error estimates. We do not adapt the mesh at each calculation; rather we adapt the mesh based on the condition that

(2.24a,b)\begin{equation} |\theta_1 - \theta_c| < 1.0^{{\circ}},\quad \theta_c = \text{atan} ((y_2 - y_1)/(x_2 - x_1)), \end{equation}

where $\theta _c$ is the computed angle based on $(x_1,y_1)$ and $(x_2,y_2)$, the coordinates of the nodes on $\varGamma _4$ directly at the contact line and immediately adjacent, respectively. The number of elements and their sizes are highly dependent on $\lambda$ and $Ca$. As an illustrative example, for steady solutions at $\lambda = 0.1$, $\chi =0.1$ and $V=5$, the resulting mesh has ${\sim}10^3$ triangular elements and ${\sim}10^5$ discretised unknowns at $Ca = Ca_{{crit}}$.

3.1. Stability of the solution branches

We now discuss the bifurcation structure and the corresponding stability results of the advancing and receding dynamic contact-line problems. Figure 4 shows the bifurcation structures in a typical advancing case (figure 4a: $\chi = 0.1, \lambda = 0.1$) and receding case (figure 4b: $\chi = 0.0, \lambda = 0.1$). Notably, our focus here is on providing insight into the stability structure, rather than necessarily probing the precise values from experimental analyses, where the slip length could be far smaller and therefore typically require more computational resources. Previous works (Vandre et al. Reference Vandre, Carvalho and Kumar2012) have shown that whilst changes in slip length can have a weak effect on $Ca_{{crit}}$, they do not qualitatively alter the physical mechanisms at play (similarly for smaller viscosity ratios, e.g. with a glycerol–air system).

Figure 4. Solution curves for the advancing (a) and receding (b) contact lines. The values of the control parameters are $\chi = 0.1$ (a) and $0$ (b) and $\lambda = 0.1$ for both panels. Individual solutions are labelled on the curve and correspond to the inset panels with the same label. The solid/dashed lines sections of the curve are the stable/unstable branches, respectively. The smaller insets show the eigenspectra for the $\boldsymbol {A}_{1}/\boldsymbol {R}_{1}$ and $\boldsymbol {A}_{2}/\boldsymbol {R}_{2}$ solutions and at the fold. Solutions $\boldsymbol {A}_{3}/\boldsymbol {R}_{3}$ are the solutions where the inflection point on the interface first becomes parallel to the plate, i.e. when ${\rm d}\kern0.7pt x_s/{\rm d}y_s = 0$. In (a), $\boldsymbol {A}_{4}$ is the solution just before the numerical calculations cease to converge. In (b), $\boldsymbol {R}_{4}$ is the solution where the interface starts to ‘overhang’ the right plate. The top-left inset in (b) shows the bifurcation curve for larger $Y$.

The solution curves are shown in the $(Ca,Y)$ projection of the solution space (see (2.21) for a definition of $Y$). The markers on the curve indicate specific solutions which are shown in the inset panels labelled $\boldsymbol {A}_{1}$–$\boldsymbol {A}_{4}$ for the advancing contact line, and $\boldsymbol {R}_{1}$–$\boldsymbol {R}_{4}$ for the receding contact line. In both the advancing and receding cases, as $Y$ increases, the solution curve experiences a fold which separates the lower branch and upper branch. The eigenspectra are real and at the fold a single eigenvalue crosses the imaginary axis, as expected. The eigenspectra also indicate that the $\boldsymbol {A}_{1}/\boldsymbol {R}_{1}$ states (all eigenvalues in the left-hand Argand plane) are ‘attractors’ of the system and $\boldsymbol {A}_{2}/\boldsymbol {R}_{2}$ states (a single eigenvalue in the right-hand Argand plane, akin to a ‘saddle-node’ state in a two-dimensional dynamical system) are weakly unstable, as seen in figure 2(b), thus numerically confirming that the lower branch is stable (solid curve) and the upper branch is unstable (dashed curve).

The interface has an inflection point near the contact point. We measure the angle at the interface inflection point (to the downwards vertical) and, as is common, define this as $\theta _{{app}}$; the apparent contact angle (Vandre et al. Reference Vandre, Carvalho and Kumar2012; Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016b). Notably, as can be seen from solutions $\boldsymbol {A}_{1},\boldsymbol {A}_{2}$ and $\boldsymbol {R}_{1},\boldsymbol {R}_{2}$ in figure 4, $\theta _{{app}}<180^\circ$ not only on the stable branch, but also immediately after the fold on the unstable branch. This is consistent with the asymptotic results of Eggers (Reference Eggers2004a), who shows that the solution at $Ca_{{crit}}$ has $\theta _{{app}} \to 180^\circ$ as $\lambda \to 0$. Consequently, for given $\lambda > 0$ we would expect to find that $\theta _{{app}} < 180^\circ$ at $Ca_{{crit}}$. These findings are consistent with Sbragaglia et al. (Reference Sbragaglia, Sugiyama and Biferale2008).

Further up the unstable branch, see $\boldsymbol {A}_{3},\boldsymbol {A}_{4}$ and $\boldsymbol {R}_{3},\boldsymbol {R}_{4}$ in figure 4, $\theta _{{app}} = 180^{\circ }$ and the interface develops a stationary point (i.e. ${\rm d}\kern0.7pt x_s/{\rm d}y_s = 0$). In the advancing case, the solution curve terminates when the interface touches the right plate and effectively appears to ‘pinch off’ off the fluid domain; see $\boldsymbol {A}_{4}$. For the receding case, the calculations stop when the size of the computational domain no longer allows a converged solution.

The steady-solution curves of the advancing and receding contact-line problems, although treated separately in figure 4, are actually two halves of the same solution space. Figure 5 shows the connection for $\lambda = 0.1,\chi = 0\text { and }0.1$, where the signed meniscus rise, $y_s(s=L) - y_s(s=0)$, is plotted against $Ca\times U$, where $U=\pm 1$ with $\pm$ corresponding to the receding ($+$)/advancing ($-$) problem. In this projection the advancing and receding curves occupy the second and fourth quadrants, respectively, and the location of the respective folds in each quadrant highlights that the advancing contact line only becomes unstable if $\chi \neq 0$, in which case $Ca_{{crit,rec}} < Ca_{{crit,adv}}$ (Marchand et al. Reference Marchand, Chan, Snoeijer and Andreotti2012; Chan et al. Reference Chan, Srivastava, Marchand, Andreotti, Biferale, Toschi and Snoeijer2013) and that, in general, $Ca_{{adv,crit}}$ and $Ca_{{rec,crit}}$ increase as $\chi \to 0$.

Figure 5. The advancing and receding steady solution curves for $\lambda = 0.1$ and $\chi = 0$ (red curves), $\chi = 0.1$ (blue curves). The horizontal axis is $Ca$ scaled by $U=\pm 1$ and the vertical axis is $y_s(s=L) - y_s(s=0)$, i.e. the signed version of $Y$. The critical $Ca$ for each problem is denoted by a dotted line. The two problems are connected through the origin.

Finally, we note that in both cases the solution curve does not experience additional bifurcations as $Y$ increases along the unstable branch. For a system where gravity is included it is known that within the lubrication approximation, the solution curve (for the receding contact line, at least) oscillates around a fixed value of $Ca = Ca^*$ (see Chan et al. Reference Chan, Snoeijer and Eggers2012), experiencing multiple saddle-node bifurcations as $Y\to \infty$. Preliminary calculations show that, if gravity is included, the oscillations are also present in the advancing/receding hybrid system, although, for brevity, we do not show the results here.

3.2. Physical interpretation of the bifurcation

As discussed in Vandre et al. (Reference Vandre, Carvalho and Kumar2013) for the advancing contact line, the fold occurs when the fluid-pressure gradients (fluid 1) are comparable to the capillary-stress gradients (see Vandre Reference Vandre2013) near the contact line, i.e. when

(3.7)\begin{equation} \frac{\partial p_1}{\partial s} \sim \frac{1}{Ca}\frac{\partial \kappa}{\partial s}. \end{equation}

This is because as $Ca\to Ca_{{crit}}$ the air-pressure gradients near the contact line will increase as the system seeks to ‘pump’ air out of the region near the contact line to maintain a steady state. Eventually these air-pressure gradients will exceed the capillary-stress gradients and the system will be unable to maintain a stable steady equilibrium. Figure 6 shows the evolution of the quantities on either side of (3.7) calculated at the inflection point for the advancing case. Here, one can see that the air-pressure and capillary-stress gradients balance close to $Ca_{{crit}}$, as seen by the intersection of the curves, which confirms the ideas of Vandre et al. (Reference Vandre, Carvalho and Kumar2013).

Figure 6. Comparison of the air-pressure gradient and the capillary-stress gradient for the steady solutions at the inflection point (IP) on the interface for the advancing contact-line problem. Parameter values are $\chi =0.1,\lambda = 0.1$.

3.3. Fold tracking

We can take advantage of the fact that at the fold bifurcation the leading eigenvalue crosses the imaginary axis to develop an algorithm for finding $Ca_{{crit}}$. We augment the system with the additional constraint

(3.8)\begin{equation} \mathrm{Re}(\sigma_1) = 0, \end{equation}

and let another control parameter come as part of the solution. It is convenient to let the interface length, $L$, be determined by (3.8) so we are able to track the evolution of $Ca_{{crit}}$ as another parameter, the viscosity ratio $\chi$, for example, is varied. This is a robust way of tracking the fold without having to recalculate the solution curve for every set of parameters, as previously considered in Kamal et al. (Reference Kamal, Sprittles, Snoeijer and Eggers2019) and Vandre et al. (Reference Vandre, Carvalho and Kumar2012).

If we vary $\chi$ and calculate $Ca_{{crit}}$ we observe that the curve of the loci of $Ca_{{crit}}$ does not itself experience any bifurcation, as seen in figure 7(a,b). In addition, we observe that the bifurcation structure also remains intact when the slip length, $\lambda$, is varied (see figure 7c) . In fact, we note that the value of $Ca_{{crit}}$ only changes relatively weakly as $\lambda$ changes from $O(1)$ to $O(10^{-4})$. Therefore we expect the dynamics to be qualitatively similar (from a dynamical systems perspective) regardless of the slip length and provided $\chi \neq 0$.

Figure 7. The evolution of the fold, and hence $Ca_{{crit}}$, as the relative viscosity $\chi$ (a,b) and the slip length $\lambda$ (c) are varied. (a) The advancing contact line, $\lambda = 0.1$. As $\chi \to 0$, $Ca_{{crit}}\to \infty$ so that the fold in the bifurcation structure ceases to exist. The inset shows the generic bifurcation structure in the $\chi =0$ and $\chi \neq 0$ cases. In the former case there is always a stable solution for the system to be attracted to. (b) The receding contact line, $\lambda = 0.1$. The fold exists for all viscosity ratios. For a given $\chi$ the bifurcation structure is shown in the inset. (c) Variation of $Ca_{{crit}}$ as $\lambda$ is varied for the advancing and receding problem (different colours) with $\chi = 0.1$.

An important observation is that the advancing and receding cases differ significantly as $\chi \to 0$. For the advancing case, $Ca_{{crit}}\to \infty$ in this limit, whilst for the receding case it tends to a finite value. This indicates that the viscosity of the fluid phase has to be taken into account for the advancing contact line in order to describe the bifurcation structure. This feature has been identified before, in driven liquid filaments (Ledesma-Aguilar et al. Reference Ledesma-Aguilar, Nistal, Hernández-Machado and Pagonabarraga2011) and in plate-plunging experiments (Marchand et al. Reference Marchand, Chan, Snoeijer and Andreotti2012). In contrast the receding contact line is essentially a one-phase problem, and the qualitative features of the bifurcation structure are the same regardless of the viscosity of the fluid.

3.4. Eigenmode perturbations

We now discuss the nature of the eigenmodes resulting from the stability analysis. The modes corresponding to the three leading eigenvalues of the unstable branch, i.e. $\sigma _1,\sigma _2\text { and }\sigma _3$, are shown in figure 8. These eigenmodes correspond to the base states $\boldsymbol {w}_{\star } = \boldsymbol {A}_2,\boldsymbol {R}_2$ in figure 4. In this figure the dotted profile indicates the steady interface shape and the coloured lines indicate the shape of the interface when it is perturbed by a single eigenmode, i.e.

(3.9)\begin{equation} \boldsymbol{w} = \boldsymbol{w}_{{\star}} \pm \rho \boldsymbol{g}_i,\quad i = 1,2,3. \end{equation}

The dashed/solid curves correspond to the $+/-$ sign, respectively. The amplitude of the perturbation, $\rho$, is constrained so that the meniscus rise of the perturbation is no more than 10 % of the rise of the steady solution. Each successive mode intersects the steady interface at precisely one more location, in similarity to the form of the (sinusoidal) eigenmodes in a related lubrication model (see figure 20). Thus, the effect of adding higher-order eigenmodes to the steady state is to add extra corrugations to the interface and increase the overall length of the interface.

Figure 8. The eigenmodes of the unstable branch. The steady profile is shown as the dotted line. (a–c) The leading eigenmodes, $\boldsymbol {g}_1,\boldsymbol {g}_2,\boldsymbol {g}_3$, corresponding to the advancing $\boldsymbol {w}_{\star } = \boldsymbol {A}_2$ ($Ca=0.873,\chi =0.1$) solution in figure 4(a). (d–f) The leading eigenmodes corresponding to the receding $\boldsymbol {w}_{\star } = \boldsymbol {R}_2$ ($Ca=0.3,\chi =0$) solution in figure 4(b). In each case the eigenmodes cross the steady interface in successively more locations. The slip length is $\lambda = 0.1$.

Concentrating on the leading eigenmode alone, the action of adding a multiple of $\boldsymbol {g}_1$ to a steady solution, i.e. using (3.9), stretches/shrinks the interface according to the $\pm$ sign with no additional corrugations. Figure 9(a) shows the stable $\boldsymbol {A}_{1}$ and unstable $\boldsymbol {A}_{2}$ states with solid lines and the perturbation from the $\boldsymbol {A}_{1}$ state using (3.9) and $i=1$ with dotted lines. This figure demonstrates that we can continuously ‘stretch’ the nonlinear stable state by increasing the strength of the perturbation, $\rho$, and can eventually achieve an interface with an identical value of $Y$ to the unstable steady state and remarkably similar profile. This will have consequences, as discussed in the next section. We can also continuously deform the unstable branch in the same manner to match the interface of the stable branch. Therefore, we denote perturbations using the leading eigenmode only in (3.9) as ‘stretch’ perturbations.

Figure 9. Perturbations of the advancing $\boldsymbol {w}_{\star }=\boldsymbol {A}_1$ state. (a) ‘Stretch’ perturbations. The lower solid curve is the stable interface and the upper solid curve is the unstable interface. The dotted curves are different strength perturbations of the stable state according to (3.9) with $i=1$. The dashed curve corresponds to the perturbation which results in an identical value of $Y$ as the $\boldsymbol {A}_2$ state. (b) ‘Corrugation’ perturbations. The coloured curves are different strength perturbations of the stable state according to (3.10) for a fixed value of $\rho$ and increasing $N$ and the shaded region represents the nonlinear $\boldsymbol {A}_{1}$ steady solution. Parameter values are $Ca = 0.873,\chi =0.1,\lambda = 0.1$.

In a physical experiment, perturbations will naturally emerge from the presence of ‘noise’ in the system causing fluctuations of the interface and contact line. We will take advantage of the stability analysis and mimic this ‘noise’ by using the higher-order eigenmodes in the perturbation, i.e.

(3.10)\begin{equation} \boldsymbol{w} = \boldsymbol{w}_{{\star}} + \sum_{i=1}^{N}\rho_i\boldsymbol{g}_i, \end{equation}

where for the purpose of simple illustration we choose the amplitude coefficients, $\rho _i$, to be equal. Figure 9(b) shows the perturbed interface as the value of $N$ increases from 1 to 10, which demonstrates that by increasing the value of $N$ we are able to perturb the nonlinear steady interface with increasingly more corrugations or ‘noise’ in a systematic manner and hence be able to isolate specific geometric effects that act on the stability of the system (the first 10 modes are sufficient to analyse the response of the system as higher-order modes will decay rapidly). Strictly speaking, to model physical ‘noise’ would require performing a number of realisations with the values of $\rho _i$ chosen randomly from a probability distribution, which must be chosen in some reasonable manner, ideas which we do not pursue further here. Henceforth, perturbations of the form in (3.10), i.e. a combination of leading-order and higher-order modes, are called ‘corrugation’ perturbations.

The two forms of perturbation discussed here, in our nomenclature the ‘stretch’ and ‘corrugation’ perturbations, are used in the next section to understand the subsequent time-dependent behaviour of the system after systematically applying a perturbation to a steady state.

4.1. Perturbations from a steady-state IVP: $Ca< Ca_{{crit}}$

We now consider the first question and look at the dynamics of the system when $Ca< Ca_{{crit}}$. Our methodology is to start at either the stable or unstable state and perturb it using either a ‘stretch’ or ‘corrugation’ eigenmode expansion, which we will consider separately. We then run a series of IVPs to examine the transient behaviour and eventual dynamical outcome.

4.1.1. ‘Stretch’ perturbations

For the ‘stretch’ perturbations, we use the leading eigenmode only, and set the initial conditions to be

(4.3)\begin{equation} \boldsymbol{w}(t=0) = \boldsymbol{w}_{{\star}} \pm \rho \boldsymbol{g}_1. \end{equation}

Initially we concentrate on the advancing contact-line problem. Using the initial conditions stated in (4.3), small perturbations from the $\boldsymbol {A}_{1}$ stable state (figure 4a) decay and the system (unsurprisingly) relaxes back to its stable configuration. Figure 11 demonstrates this by tracking the value of $Y$ in time for two different perturbations, corresponding to $\Delta (0,\boldsymbol {A}_1)<0$ and $\Delta (0,\boldsymbol {A}_1)>0$, of the stable state near the fold (parameter values quoted in the caption). Furthermore, as seen in the insets on the right of figure 11, the decay rate excellently matches the value of the leading eigenvalue, $\sigma _1$, obtained from the linear stability analysis.

Figure 11. Time-dependent perturbations of the advancing contact line using (4.3). This figure shows perturbations from the solution labelled $\boldsymbol {A}_{1}$ in figure 4(a). The two different initial conditions are shown in the inset panels corresponding to the same colour lines in the main panel. Whether $\Delta (0,\boldsymbol {A}_1)>0$ or $\Delta (0,\boldsymbol {A}_1)<0$ the system relaxes back to the stable state. The inset panels on the right show the time signal of $\log (|\Delta (t,\boldsymbol {A}_1)|)$ as compared to the predicted growth rate $\sigma _1$. Parameter values are $Ca = 0.873,\chi =0.1,\lambda = 0.1$.

For the same parameter values, we can also perturb the $\boldsymbol {A}_{2}$ unstable state (figure 4a) by its leading eigenmode, which is unstable. Figure 12(a,b) shows the time signal of $Y$ as well as the phase-plane trajectories in the $(\overline {Ca},Y)$ projection. This case is more interesting, and we see that if $\Delta (0,\boldsymbol {A}_2)<0$ (i.e. the perturbation ‘contracts’ the $\boldsymbol {A}_{2}$ interface) then the system returns to the stable state (see figure 12a), whereas if $\Delta (0,\boldsymbol {A}_2) > 0$ (i.e. a ‘stretch’) then entrainment of the upper fluid occurs (see figure 12b).

Figure 12. Time-dependent perturbations of the unstable branch using (4.3). (a,b) The advancing contact-line problem. (c,d) The receding contact-line problem. This figure shows perturbations from the solution labelled $\boldsymbol {A}_{2}/\boldsymbol {R}_{2}$ in figure 4. For both the advancing and receding contact line, when $\Delta (0,\boldsymbol {A}_2)<0$ the perturbation relaxes to the corresponding stable branch (a,c). When $\Delta (0,\boldsymbol {A}_2)>0$ the perturbation grows and in the case of the advancing contact line, entrainment occurs; see inset panel labelled ‘Entrainment’ in (b). For the receding contact line a thin film develops; see inset panel labelled ‘Thin-film’ in (d). The other inset panels show the phase plane in the $(\overline {Ca},Y)$ projection. The initial conditions are marked as open circles, the system trajectory is marked with arrows and the steady bifurcation is shown without arrows. The unstable branch can be considered as the basin boundary of attraction of the stable steady state. Parameter values are $Ca = 0.873,\chi =0.1,\lambda = 0.1$ for the advancing problem and $Ca = 0.3,\chi =0,\lambda = 0.1$ for the receding problem.

Similar outcomes occur for the receding contact-line problem. Perturbations, using the leading eigenmode, of the $\boldsymbol {R}_{1}$ (figure 4b) stable state relax back to the stable equilibrium (results not shown). In contrast, figure 12(c,d) shows the time evolution of perturbations from the $\boldsymbol {R}_{2}$ (figure 4b) unstable state, where for $\Delta (0,\boldsymbol {R}_2) < 0$, the system relaxes back to the $\boldsymbol {R}_{1}$ state. But if $\Delta (0,\boldsymbol {R}_2)>0$ then, unlike the advancing contact line where entrainment occurs, a thin film develops that grows in size at a linear rate as $t\to \infty$ as shown in figure 12(d).

These IVP calculations show that if we consider the class of perturbations representing ‘stretches’, using the leading eigenmode only, then the indicator of whether the system returns to the stable state is that meniscus rise of the initial perturbation is smaller than the meniscus rise of the unstable state, i.e. the condition for stability is

(4.4)$$\begin{gather} \Delta(0,\boldsymbol{A}_2) < 0,\quad\text{advancing contact line}, \end{gather}$$

(4.5)$$\begin{gather}\Delta(0,\boldsymbol{R}_2) < 0,\quad\text{receding contact line}. \end{gather}$$

In the language of dynamical systems the unstable state represents the ‘boundary of the basin of attraction’ of the stable state, when considering simple stretches of the stable interface corresponding to perturbations using the leading eigenmode. The unstable state is therefore not just a trivial consequence of the steady bifurcation structure but also has an important role in dividing the phase plane into regions that have different transient dynamics.

Another important observation is that the calculation in figure 12(d), for the receding contact line, is consistent with the claim from Chan et al. (Reference Chan, Snoeijer and Eggers2012) that the solution curve represents the effective dynamics of the system ‘in which the state of the solution moves quasi-statically along the solution curve’. In the receding case, the trajectory of the system in the phase space $(\overline {Ca},Y)$ is qualitatively similar to the steady solution curve when plotted in the same diagram. The inset panels in figure 12(b,d) labelled ‘Phase plane’ show the steady solution curve and the trajectory of the system shown with an arrow. In the receding case, when $\Delta (0,\boldsymbol {R}_2)>0$ the trajectory closely follows the upper reaches of the unstable branch, whereas in the advancing problem this similarity does not occur. The comparison between the steady solution curve and the time-dependent trajectories is discussed in more detail in the next section. Finally, we note that the qualitative behaviour and conclusions that we have described are identical to the full model, as we show in Appendix C.

4.1.2. ‘Corrugation’ perturbations

We now concentrate on the ‘corrugation’ perturbations which arise from taking more eigenmodes in the perturbation expansion (3.10). As discussed in the previous section this class of perturbations includes higher modes, which we might reasonably expect in a noisy physical system. We choose the number of eigenmodes as $N=10$ for computational efficiency. We set the initial condition of the system therefore as

(4.6)\begin{equation} \boldsymbol{w}(t=0) = \boldsymbol{w}_{{\star}} + \rho\sum_{i=1}^{N}\boldsymbol{g}_i + \text{c.c.}, \quad N = 10. \end{equation}

We concentrate solely on the advancing case and, as before, examine perturbations from the stable $\boldsymbol {A}_{1}$ and unstable $\boldsymbol {A}_{2}$ (figure 4a) states using the initial conditions prescribed in (4.6). Figure 13 shows the time-dependent results for incrementally increasing values of $\rho$ in (4.6). Figure 13(a) shows trajectories in the $(L,Y)$ phase-plane projection, where $L$ is the total arclength of the interface. To help orient the trajectories we have added artificial axes which are centred on the unstable $\boldsymbol {A}_{2}$ state.

Figure 13. Time-dependent perturbations using the ‘corrugation’ perturbations. Perturbations from the solutions labelled in figure 4(a) using (4.6). (a–c) Perturbations from $\boldsymbol {A}_{1}$. (d–f) Perturbations from $\boldsymbol {A}_{2}$. (a,d) Time-dependent trajectories in the in the $(L,Y)$ phase-plane projection. Initial perturbations (IC) are denoted by open circles and steady states by large solid circles. The dashed blue lines indicate the unstable manifold of the $\boldsymbol {A}_{2}$ state and the solid black lines are artificial axes centred on the unstable $\boldsymbol {A}_{2}$ state. (a) Each red curve indicates a time-dependent trajectory with a different value of increasing $\rho$, as indicated by the black arrow, the initial conditions with minimum and maximum $\rho$ have been labelled. The solid trajectories correspond to initial perturbations that return to the steady state and dotted trajectories correspond to trajectories that result in entrainment. The inset shows the linear response given in (3.6) for a particular initial condition marked with a filled circle in the main panel. (b) The initial interfaces of the trajectories shown in (a) are shown as solid lines and the $\boldsymbol {A}_{1}$ and $\boldsymbol {A}_{2}$ steady states are shown using the shaded area and a dashed line, respectively. (c) The corresponding time signals of $Y$ for the trajectories shown in (a). (d) Red/black curves indicate trajectories with initial conditions $\Delta (0,\boldsymbol {A}_2)$ positive/negative, respectively. (e, f) The initial conditions corresponding to the trajectories in (d) with the unstable state indicated by the shaded area. Parameter values are $Ca = 0.873,\chi =0.1,\lambda = 0.1$.

After perturbing the $\boldsymbol {A}_{1}$ state, initially the system ‘smoothes’ the corrugations, which is indicated by the decreasing value of $L$ in the trajectories. Once these corrugations disappear, either the system becomes transient, and air entrainment occurs, or the system returns to the stable state. The dashed line between the $\boldsymbol {A}_{1}$ and $\boldsymbol {A}_{2}$ state is the response of the system if only the leading eigenmode is retained in the initial perturbation, as considered in the previous section on ‘stretch’ modes, which can be interpreted as the approximate unstable manifold of the $\boldsymbol {A}_{2}$ state.

The nonlinear trajectories, i.e. the time-dependent solution to the hybrid system given in (2.3)–(2.14) and (2.17a,b) with initial conditions (3.10), all eventually collapse on the unstable manifold of the unstable $\boldsymbol {A}_2$ state once the higher-order modes have sufficiently decayed, but, significantly, the combination of stable higher-order modes in the initial perturbation can cause transient behaviour, i.e. fluid entrainment, in the system, despite the corresponding eigenvalues being highly stable. This is in contrast to the linear response, i.e.

(4.7)\begin{equation} \boldsymbol{w}_{{lin}}(t) = \boldsymbol{w}_{{\star}} + \rho\sum_{i=1}^{N}\boldsymbol{g}_i\mathrm{e}^{\sigma_it} + \text{c.c.}, \quad N = 10, \end{equation}

which will always decay back to the stable state (as the eigenvalues are all negative), as demonstrated in the inset diagram for a particular initial condition.

The unstable $\boldsymbol {A}_{2}$ state, as indicated by a large circular symbol, plays a crucial role in the partition of behaviours. By examining the trajectories that explore the vicinity of the unstable branch in figure 13(a–c), it is not unreasonable to hypothesise that by continually refining the value of $\rho$ in the initial perturbation the system would be able to stay in the vicinity of the unstable state for an arbitrary time period, where the dynamics of the system is dominated by the stable eigenmodes of the unstable branch. This behaviour is typically indicative of interpreting the unstable branch as an ‘edge’ state, commonly used to describe weakly unstable states in the transition to turbulence and other fluid dynamics problems (see Kerswell, Pringle & Willis (Reference Kerswell, Pringle and Willis2014) for a description of an edge state). In these scenarios the weakly unstable ‘edge’ state acts as an ‘edge’ between two dynamical outcomes; in our problem it separates the system returning to the stable state or becoming transient.

The role of the unstable state is further emphasised in figure 13(d–f) for perturbations from the unstable $\boldsymbol {A}_{2}$ branch. In figure 13(d) different colour trajectories correspond to whether a positive or negative $\rho$ is chosen in the initial conditions given by (4.6), and figure 13(e, f) shows the initial conditions compared to the $\boldsymbol {A}_{2}$ steady state. Again, the system outcomes are partitioned by the presence of the unstable state which ‘deflects’ the system to either relax back to the stable state or become transient so that entrainment occurs. We note that in this case if $\Delta (0,\boldsymbol {A}_2)<0$ (initial conditions in the lower half-plane of figure 13d) then the system will become transient which is in direct contradiction to the ‘stretch’ mode perturbations. These calculations demonstrate that $L$ could be a more robust indicator of dynamical outcome; in this figure, all trajectories in the right-side $(L,Y)$ plane (centred on the unstable state) eventually become transient as entrainment occurs.

To determine the exact regions in phase space where initial conditions becomes transient or return to the stable state requires knowledge of the stable manifold of the $\boldsymbol {A}_2$ state (see figure 2b). In our high-dimensional system of PDEs, the calculation of the stable manifold of the $\boldsymbol {A}_2$ state is a non-trivial task (see e.g. Krauskopf et al. Reference Krauskopf, Osinga, Doedel, Henderson, Guckenheimer, Vladimirsky, Dellnitz and Junge2005) which we do not pursue here.

4.1.3. Physical significance of the perturbations

We now relate these results to the physical system. Firstly, we make the key observation that the system is able to experience instability for $Ca< Ca_{{crit}}$ due to finite-amplitude perturbations of the stable state. Furthermore, these results suggest that the system is susceptible to instability caused by perturbations that increase the total arclength; corrugations, representing ‘noise’, can be just as dangerous as perturbations that increase the meniscus rise. Given that noise is a ubiquitous phenomenon in physical systems we would expect to see this realised in a system with $Ca< Ca_{{crit}}$. Secondly, the extent to which we are able to perturb the stable $\boldsymbol {A}_{1}/\boldsymbol {R}_{1}$ state so that the system remains stable is dictated by the information encoded in the unstable $\boldsymbol {A}_{2}/\boldsymbol {R}_{2}$ steady state, in particular whether the perturbation causes the interface length or meniscus rise to increase beyond that of the unstable steady state. A consequence of this is that the closer $Ca$ is to $Ca_{{crit}}$ (from below), the smaller is the finite amplitude required to cause the system to become transient, as the stable and unstable branches are increasingly approaching each other.

We can be more precise and approximate the ‘size’ of the maximal stretch perturbation that maintains a steady contact line as a function of $Ca-Ca_{{crit}}$. Locally to the fold, the bifurcation curve will take the form

(4.8)\begin{equation} Ca - Ca_{{crit}} = a_1(Y - Y_{{crit}})^2 + a_2(Y - Y_{{crit}})^3 + \cdots, \end{equation}

where $Y_{{crit}}$ is the meniscus rise at the fold and $a_i$ are constants that can be determined by, for example, normal-form reduction methods (Kuznetsov Reference Kuznetsov1998). Therefore, close to the fold, the ‘length’ of the basin boundary is approximately equal to the vertical distance between the stable and unstable branches and is approximately

(4.9)\begin{equation} Y_{2} - Y_{1} = \sqrt{\frac{2}{a_1}}|Ca - Ca_{{crit}}|^{1/2}. \end{equation}

As a conclusion, finite disturbances from the steady stable state potentially lead to instability before that interface might become transient otherwise (due to a lack of steady state). Minimising ambient disturbances and minimising fluctuations in $Ca$ would help maintain a stable interface for higher $Ca$. The physical mechanisms underlying the observations are not straightforward to address due to the highly nonlinear nature of the problem. Moreover, three-dimensional perturbations will likely be important in practice and introduce additional physical mechanisms. Such three-dimensional calculations are exceedingly challenging and are currently the focus of ongoing research, so we defer a detailed study of physical mechanisms for future work.

4.2. Dynamics ‘beyond’ the fold: $Ca>Ca_{{crit}}$

We now turn our attention to starting the system from rest with a flat interface beyond the fold, i.e. $Ca>Ca_{{crit}}$, so that we are able to answer the second question stated at the start of § 4. The initial condition is

(4.10)\begin{equation} \boldsymbol{w}(t=0) = [\boldsymbol{u}(0),p(0),p_1(0),\boldsymbol{r}(0)]^{\rm T} = [\boldsymbol{0},\boldsymbol{0},\boldsymbol{0},(x_s,0)]^{\rm T}. \end{equation}

For values of $Ca$ exceeding the critical value there are no (known) steady states which can influence the system, unlike in the previous section. For the advancing case, the system becomes transient with entrainment occurring, as shown in figure 14 where the time signal of $Y$ is measured along with time snapshots of the interface at the times indicated for $Ca = 0.95,1.05$. We emphasise that the system is not attracted to a different steady state of the system that may exist and we are unable to find any additional steady states beyond the fold bifurcation. There are two distinct length scales present in the interface profile for sufficiently large times. The width of the fluid entrainment, $\hat {h}$, and the width of the thin film, $h_{{film}}$, appear to be weakly dependent on $Ca$ as shown by the inset of figure 14(b). These results are consistent with recent experimental work (He & Nagel Reference He and Nagel2019), but we leave a detailed analysis of the thickness of these entrainment regions, and their dependence on system parameters, to a future study; our focus is on the broad dynamical outcome, rather than the fine details of the solution. We do, however, remark that the hybrid model shows excellent promise in being able to predict the width of these films.

Figure 14. Time-dependent evolution for an advancing contact line when $Ca>Ca_{{crit}}$ and the system is at rest initially for $\lambda = 0.1,\chi = 0.1$. The different colour curves indicate different values of $Ca = 0.95,1.05$. (a) The time signal of $Y$ with interface profiles at the indicated times in the inset. (b) The trajectories (curves with arrows) in the $(\overline {Ca},Y)$ phase-plane projection compared with the steady bifurcation curve (curve without arrows). The inset compares the interface profiles for each $Ca$ at $t = 200$.

For the receding case, a thin film develops. Figure 15(a) shows the time signals of $Y$ in the receding case when $Ca = 0.4,0.5 > Ca_{{crit}}$. The system, although not approaching a steady state, forms a coherent structure whose meniscus rise grows at a constant velocity, independent of $Ca$, similar to the structure seen in figure 12(c). Unlike the advancing case, the time-dependent trajectories and the steady solution curve for the receding problem, when plotted in $(\overline {Ca},Y)$, are qualitatively similar. Figure 15(b) also compares the actual interface profiles for the time-dependent calculation and the corresponding steady solution for the same value of $Y$ (insets). It is striking how close the two corresponding interface profiles are, although we note that in the time-dependent case the horizontal height of the inflection point is smaller than that of the steady solution in all of the cases. This again shows compelling evidence that in the receding contact-line problem the unsteady solution branch represents the effective dynamics of the system.

Figure 15. Time-dependent evolution for a receding contact line when $Ca>Ca_{{crit}}$ and the system is at rest initially for $\lambda = 0.1,\chi =0$. The different colour curves indicate different values of $Ca = 0.4,0.5$. (a) The time signal of $Y$ with the final interface profiles in the inset. (b) The trajectories (curves with arrows) in the $(\overline {Ca},Y)$ phase-plane projection compared with the steady bifurcation curve (curve without arrows). The inset panels compare the interface profiles for the time-dependent evolution (red) and the steady solution branches (blue) at the same values of $Y = 1,2,3$ for the $Ca = 0.4$ case.

In both the advancing and receding contact-line problems, the system trajectory approaches a fixed value of $\overline {Ca}$ as $t\to \infty$. For the advancing contact line, the limiting meniscus rise velocity is dependent on $Ca$ (see slopes in figure 14a) so that the trajectories approach a limiting value of $\overline {Ca}$ which depends on $Ca$ (see figure 14b). However, for the receding contact line, the limiting value of $\overline {Ca}$ appears to be almost independent of $Ca$ (see figure 15b), indicating the contact-line region and the thin-film region are essentially de-coupled by this point. This would mean that if, after a thin film has developed, we slow the plate speed so that $Ca< Ca_{{crit}}$, the contact line would continue to move upwards and only the thin-film thickness, which is dependent on $Ca$ according to a Landau–Levich-type law (Landau & Levich Reference Landau and Levich1988), would change. Figure 16 shows this situation, where the trajectory in the $(\overline {Ca},Y)$ plane of a receding contact-line IVP initially starting at rest when $Ca>Ca_{{crit}}$. Once a film is developed we instantaneously change the value of $Ca< Ca_{{crit}}$, once $Y$ is large enough (we choose $Y = 6$ for convenience) and observe the contact line continuing to rise, whilst the thin film becomes smaller, adjusting to the changed value of $Ca$.

Figure 16. Time-dependent evolution for the receding contact line when $Ca>Ca_{{crit}}$ and the system is at rest initially. Once a film has formed and $Y = 6$, we set $Ca< Ca_{{crit}}$, as indicated in the figure by different colour trajectories and interface profiles. Other parameter values are $\lambda = 0.1,\chi = 0.0$.

We do not have a concrete physical explanation for why the trajectories and steady bifurcation curve closely match one another for the receding contact line. From a dynamical systems perspective it is possible that the stable manifolds of the unstable branch persist as $Ca > Ca_{{crit}}$ and that the system is constantly being ‘attracted’ by the remnants of the stable manifolds. This is, of course, speculation, and we leave this aspect as a challenge for future work.