1. Introduction
The process of dip coating aims to deposit a thin liquid layer on the surface of an object by withdrawing the latter from a liquid bath in which it is initially immersed. This simple technique is employed in many industrial processes in order to modify the properties of a solid surface (Scriven Reference Scriven1988), with applications ranging from anti-corrosion treatments and optically modified glasses, to surface functionalisation of bio-implants. The hydrodynamics of dip coating has been known for nearly 80 years. The original description of this process dates back to the now classical work of Landau & Levich (Reference Landau and Levich1942), followed by the contribution of Derjaguin (Reference Derjaguin1943), who included the effect of gravity. We will refer to their description of the one-liquid dip-coating problem as the LLD theory. Later on, when the relevant mathematical tools became available, the rigorous asymptotic theory underlying these initial developments was proposed by Wilson (Reference Wilson1982).
Nurtured by practical applications, experimental and theoretical aspects of dip coating have continued to attract scientific attention ever since, as highlighted by the recent review by Rio & Boulogne (Reference Rio and Boulogne2017). New effects and regimes have been found to arise from the use of partially wetting (Snoeijer et al. Reference Snoeijer, Ziegler, Andreotti, Fermigier and Eggers2008; Tewes et al. Reference Tewes, Wilczek, Gurevich and Thiele2019) or textured (Seiwert, Clanet & Quére Reference Seiwert, Clanet and Quére2011) solid substrates, non-Newtonian liquids (Hurez & Tanguy Reference Hurez and Tanguy1990; Maillard et al. Reference Maillard, Bleyer, Andrieux, Boujlel and Coussot2016), surface active molecules (Park Reference Park1991; Scheid et al. Reference Scheid, Delacotte, Dollet, Rio, Restagno, Van Nierop, Cantat, Langevin and Stone2010; Champougny et al. Reference Champougny, Scheid, Restagno, Vermant and Rio2015) or jammed hydrophobic micro-particles (Dixit & Homsy Reference Dixit and Homsy2013; Ouriemi & Homsy Reference Ouriemi and Homsy2013) adsorbed at the liquid/gas interface. The two latter cases can be regarded as stepping stones to multi-phase dip-coating situations. Indeed, interfacial particle rafts were shown to behave as two-dimensional elastic sheets floating on top of the liquid bath (Vella, Aussillous & Mahadevan Reference Vella, Aussillous and Mahadevan2004), while surfactant monolayers at liquid/gas interfaces are known to exhibit analogues to two-dimensional liquid or gaseous phases (Vollhardt & Fainerman Reference Vollhardt and Fainerman2010).
Beyond these analogies, the dip-coating configuration in which the bath contains two immiscible liquids, one floating on top of the other, has never been explored. The objective of the present paper is to extend the LLD dip-coating theory to describe liquid entrainment at such a gas/liquid/liquid compound interface. Processes occurring at compound interfaces, for example made of a layer of oil floating on water, raise interest in the context of environmental science or semiconductor electronics. The surface of the oceans can be seen as a compound interface, due to the presence of the sea surface microlayer (Hardy Reference Hardy1982; Liss & Duce Reference Liss and Duce2005) and even more so in the event of an oil spill (Fingas Reference Fingas2015), hence the relevance of processes such as bubble bursting (Feng et al. Reference Feng, Roché, Vigolo, Arnaudov, Stoyanov, Gurkov, Tsutsumanova and Stone2014; Stewart et al. Reference Stewart, Feng, Kimpton, Griffiths and Stone2015) or bouncing (Feng et al. Reference Feng, Muradoglu, Kim, Ault and Stone2016) at gas/liquid/liquid interfaces. These are examples of elementary processes that occur at the millimetric or sub-millimetric scale in the ocean and whose physics is related to the one we consider in this work. In a different context, dip coating through a compound interface, consisting of lifting a substrate through a layer of carbon-nanotube-laden ink floating on top of a water bath, was experimentally investigated by Jinkins et al. (Reference Jinkins, Chan, Brady, Gronski, Gopalan, Evensen, Berson and Arnold2017). They showed the potential of this method, known as floating evaporative self-assembly in the context of semiconductor electronics, for the deposition of well-aligned carbon nanotube arrays.
In the present study, we will restrict ourselves to the situation in which both liquid phases are dragged, giving birth to a superposition of two liquid films on the substrate. It is worth mentioning the connection of this configuration to the model introduced by Seiwert et al. (Reference Seiwert, Clanet and Quére2011) to describe the dip coating of a textured solid with a single liquid phase. In that work, the effect of the texture was modelled as a secondary, uniform layer made of a fluid with a viscosity higher than that of the coating liquid. Finally, in a different geometry, the entrainment of a thin liquid film at a liquid/liquid interface has been described in the context of lubricant-infused surfaces (Li et al. Reference Li, Ueda, Paulssen and Levkin2019). When a water droplet is sliding on an oil-imbibed surface, thin oil films observed behind, beneath or around the droplet have been shown to obey similar scaling laws as in the LLD theory (Daniel et al. Reference Daniel, Timonen, Li, Velling and Aizenberg2017; Kreder et al. Reference Kreder, Daniel, Tetreault, Cao, Lemaire, Timonen and Aizenberg2018).
Our theoretical investigation of dip coating through a gas/liquid/liquid compound interface, leading to the deposition of a double liquid coating, is organised as follows. In § 2, we develop the problem formulation and comment on the various dimensionless control parameters. Our results are presented and discussed in §§ 3 and 4. In § 3, we focus on the typical interfacial shapes and flow structures obtained from the numerical solutions of the model. These observations allow us to rationalise quantitatively the asymptotic film thicknesses coated on the substrate, highlighting the universality of the entrainment mechanism (viscous stresses vs capillary suction). In § 4, we examine the dependence of the film thicknesses on a number of dimensionless control parameters of the problem. Importantly, we show that the existence of a physically meaningful second film is limited to specific, finite regions in the parameter space. In § 5, we comment on the scope and limitations of our model in relation to disjoining pressure effects and partial wetting conditions between the two liquids. Finally, § 6 is devoted to the conclusions.
2. Description of the model
2.1. Dimensional flow equations
We present here the dimensional formulation of the two-liquid dip-coating problem sketched in figure 1. A solid plate is vertically lifted at a constant speed 
$U$ through a stratified bath made of two immiscible liquid layers. This bath consists of a pool of liquid 1 (with density 
$\rho _1$ and viscosity 
$\mu _1$), covered by a layer of liquid 2 (with density 
$\rho _2 < \rho _1$ and viscosity 
$\mu _2$), which has a thickness 
${\rm \Delta} H$ far from the plate. The interface between liquid 1 and 2, denoted by (I), has an interfacial tension 
$\sigma _{12}$, while the interface between liquid 2 and the surrounding air, denoted by (II), has a surface tension 
$\sigma _{2a}$. Both liquids 1 and 2 are supposed to perfectly wet the plate. However, at this stage, we do not impose perfect wetting conditions between liquids 1 and 2. In other words, the spreading factor 
$S = \sigma _{1a} - (\sigma _{2a} + \sigma _{12})$, where 
$\sigma _{1a}$ denotes the surface tension of liquid 1, can be either positive or negative for the moment (see e.g. De Gennes, Brochard-Wyart & Quéré Reference De Gennes, Brochard-Wyart and Quéré2013). We will discuss the effect of the wetting properties of liquid 2 on liquid 1 in § 2.6 and later on in § 5.1.
Sketch of the flow configuration described in this work: a solid plate is pulled vertically at constant speed through a compound bath made of a lighter liquid (2) on top of a denser one (1). The liquid 1/liquid 2 and liquid 2/air interfaces are denoted by (I) and (II), respectively.
We consider the problem to be steady and two-dimensional in space, described by the horizontal coordinate 
$x$ and the vertical coordinate 
$z$ (see figure 1). As the plate crosses the compound bath, we assume it entrains two thin films: a lower film of liquid 1 with thickness 
$h_1(z)$ and an upper film of liquid 2 with thickness 
$\delta h (z) = h_2(z) - h_1(z)$. The velocity fields are 
$u_1 \boldsymbol {e_z} + v_1 \boldsymbol {e_x}$ and 
$u_2 \boldsymbol {e_z} + v_2 \boldsymbol {e_x}$ in liquids 1 and 2, respectively.
Following the approach by Landau & Levich (Reference Landau and Levich1942) and Derjaguin (Reference Derjaguin1943), we describe the flow of each phase in the transition region, referred to as the dynamic meniscus, that connects the films of constant thickness downstream to the static meniscus upstream. For each liquid, the vertical extent of the static meniscus is assumed much larger than the corresponding film thickness, allowing us to apply the lubrication theory to the flow. In this steady lubrication approach, the momentum conservation equations in the 
$z$-direction read
$$\begin{gather} 0 ={-}\frac{\partial p_1}{\partial z} + \mu_1 \frac{\partial^2 u_{1}}{\partial x^2} - \rho_1 g,\quad (0 < x < h_1(z),\ z > 0) \end{gather}$$
$$\begin{gather}0 ={-}\frac{\partial p_2}{\partial z} + \mu_2 \frac{\partial^2 u_{2}}{\partial x^2} - \rho_2 g, \quad (h_1(z) < x < h_2(z),\ z > {\rm \Delta} H). \end{gather}$$
The 
$x$-independent pressures 
$p_i$ (
$i = 1, 2$) are related to the atmospheric pressure 
$p_a$ and to the interfacial curvatures 
$\kappa _i$ through
$$\begin{gather} p_1 = p_2 - \sigma_{12} \kappa_1,\quad (0 < x < h_1(z),\ z > 0) \end{gather}$$
$$\begin{gather}p_2 = p_a - \sigma_{2a} \kappa_2,\quad (h_1(z) < x < h_2(z),\ z > {\rm \Delta} H), \end{gather}$$
$$\begin{gather}p_2 = p_a + \rho_2 g ({\rm \Delta} H - z),\quad (x > h_1(z),\ 0 < z \le {\rm \Delta} H), \end{gather}$$
where, for 
$i = 1, 2$,
\begin{equation} \kappa_i = \frac{\partial^2h_i}{\partial z^2} \left[1 + \left(\frac{\partial h_i}{\partial z}\right)^2\right]^{{-}3/2}. \end{equation}
Denoting 
$Q_1$ the flow rate of liquid 1 (i.e. under interface (I)) and 
$Q_2$ the total flow rate of liquids 1 and 2 (i.e. under interface (II)), the quasi-steady thickness-averaged continuity equations read
\begin{equation} \frac{\partial Q_i}{\partial z} = 0 \ \mathrm{for} \ i = 1, 2. \end{equation}
The expressions relating the flow rates 
$Q_i$ with the flow velocities, 
$u_i$ and the locations of the interfaces 
$h_i$ will be provided below, once the dimensionless formulation is introduced.
2.2. Non-dimensionalisation
A convenient choice for the length scale to non-dimensionalise the problem is the capillary length 
$\ell _c$ based on the properties of liquid 1 in the absence of the buoyancy effect created by liquid 2, defined as
\begin{equation} \ell_c = \sqrt{\frac{\sigma_{12}}{\rho_1 g}}. \end{equation}
Note that 
$\ell _c$ neither represents the scale of the static meniscus of liquid 1, nor the one of liquid 2. The ‘physical’ dimensional capillary lengths associated with the menisci of liquids 1 and 2 are respectively 
$\ell _{c,1} = \sqrt {\sigma _{12}/(\rho _1-\rho _2)g}$ and 
$\ell _{c,2} = \sqrt {\sigma _{2a}/\rho _2 g}$, whose dimensionless versions indeed appear in the static configuration presented in § 2.4 ((2.34) and (2.35), respectively). The independent and dependent variables of the problem are made dimensionless as follows:
$$\begin{gather} (x, z) \quad\longrightarrow\quad \ell_c\ (x, z) \end{gather}$$
$$\begin{gather}(h_1,h_2) \quad\longrightarrow\quad \ell_c\ (h_1, h_2), \end{gather}$$
$$\begin{gather}(p_1,p_2) \quad\longrightarrow\quad \sigma_{12}/\ell_c\ (p_1, p_2), \end{gather}$$
$$\begin{gather}(u_1,u_2,v_1, v_2) \quad\longrightarrow\quad U\ (u_1, u_2,v_1, v_2). \end{gather}$$
The thickness of the layer of liquid 2 is also scaled by 
$\ell _c$, still keeping the same notation: 
${\rm \Delta} H \longrightarrow \ell _c {\rm \Delta} H$. The capillary number is introduced as 
$Ca = \mu _1 U / \sigma _{12}$. Also, the following dimensionless parameters are defined:
\begin{equation} M = \frac{\mu_2}{\mu_1},\quad \varSigma = \frac{\sigma_{2a}}{\sigma_{12}},\quad R = \frac{\rho_2}{\rho_1}. \end{equation}Based on these definitions, the dimensionless momentum equations deduced from (2.1) and (2.2) read
$$\begin{gather} 0 = \varPi_1 + Ca \, \frac{\partial^2 u_{1}}{\partial x^2},\quad (0 < x < h_1(z),\ z > 0) \end{gather}$$
$$\begin{gather}0 = \varPi_2 + M \, Ca \, \frac{\partial^2 u_{2}}{\partial x^2},\quad (h_1(z) < x < h_2(z),\ z > {\rm \Delta} H) , \end{gather}$$where the dimensionless pressure gradients are
$$\begin{gather} \varPi_1 = \frac{\partial \kappa_1}{\partial z} + \varSigma \frac{\partial \kappa_2}{\partial z} - 1,\quad (0 < x < h_1(z),\ z > {\rm \Delta} H) \end{gather}$$
$$\begin{gather}\varPi_2 = \varSigma \frac{\partial \kappa_2}{\partial z} - R,\quad (h_1(z) < x < h_2(z),\ z > {\rm \Delta} H). \end{gather}$$
For 
$0 < z \le {\rm \Delta} H$, (2.16) must be replaced by
\begin{equation} \varPi_1 = \frac{\partial \kappa_1}{\partial z} - 1 + R. \end{equation}
Notice that the pressure gradient 
$\varPi _1$ and its derivative are continuous at 
$z = {\rm \Delta} H$, as 
$\varPi _2$ smoothly approaches 0 as 
$z \rightarrow 0^+$. In the above expressions, 
$\kappa _1$ and 
$\kappa _2$ are the dimensionless curvatures, which are still given by (2.6) after non-dimensionalisation. Equations (2.14) and (2.15) can be integrated with respect to 
$x$ across the corresponding layers, using the following boundary conditions:
$$\begin{gather} \text{Stress-free at } x = h_2:\quad \frac{\partial u_{2}}{\partial x} = 0, \end{gather}$$
$$\begin{gather}\text{Velocity continuity at } x = h_1:\quad u_1 = u_2, \end{gather}$$
$$\begin{gather}\text{Stress continuity at } x = h_1:\quad \frac{\partial u_{1}}{\partial x} = M \frac{\partial u_{2}}{\partial x}, \end{gather}$$
$$\begin{gather}\text{No slip at } x = 0:\quad u_1 = 1. \end{gather}$$
After obtaining the vertical velocity fields 
$u_{1}$ and 
$u_{2}$ (whose expressions are given in Appendix A), we compute the dimensionless flow rate 
$Q_1$ in the lower layer and the total dimensionless flow rate 
$Q_2$ carried by the two layers
\begin{equation} Q_i = Ca \, \int_0^{h_i} u_i (x) \, \mathrm{d}x = Ca\,h_i + \varPi_1 F_{i1} + \varPi_2 F_{i2}, \end{equation}
with 
$i = \{1, 2\}$ and
$$\begin{gather} F_{11} = \frac{h_1^3}{3} \end{gather}$$
$$\begin{gather}F_{12} = \left(h_2-h_1\right)\frac{h_1^2}{2} \end{gather}$$
$$\begin{gather}F_{21} = \frac{h_1^3}{3} + \left(h_2-h_1\right)\frac{h_1^2}{2} \end{gather}$$
$$\begin{gather}F_{22} = \frac{\left(h_2-h_1\right)^3}{3M} + h_1\left(h_2-h_1\right)^2 + \left(h_2-h_1\right) \frac{h_1^2}{2}. \end{gather}$$
Replacing the flow rates Q i given by (2.23) in the continuity (2.7) – which remain unchanged upon non-dimensionalisation – we arrive at a system of two fourth-order ordinary differential equations. This system is closed by implementing eight boundary conditions as follows. Far up from the reservoir, we impose that the thicknesses of the coated films converge asymptotically to constants: for 
$i = \{1, 2\}$, 
$\partial h_i/\partial z = \partial ^2 h_i/\partial z^2 = 0$ when 
$z \rightarrow +\infty$. Towards the reservoir, the static menisci connect to the flat surfaces of liquids 1 and 2. We therefore impose that 
$\kappa _1=0$ and 
$\partial h_1/\partial z \rightarrow \infty$ at 
$z = 0$ for liquid 1; and 
$\kappa _2 = 0$ and 
$\partial h_2/\partial z \rightarrow \infty$ at 
$z = {\rm \Delta} H$ for liquid 2.
2.3. Scaling and simplification
The problem consisting of (2.7), along with (2.23) and the above-mentioned boundary conditions, could be solved numerically for 
$h_1$ and 
$h_2$. However, in order to gain physical insight into the scalings governing the different regions of the flow, we further simplify the problem using the matched asymptotic expansion treatment proposed by Wilson (Reference Wilson1982) in the case of the one-liquid LLD flow.
The film thicknesses 
$h_i$ and the vertical coordinate 
$z$ are rescaled with a dimensionless parameter 
$\varepsilon$ as 
$h_i = \varepsilon \hat {h}_i$ and 
$z = \varepsilon ^{\alpha } \hat {z}$. This technique will allow us to find the spatial scales at which viscous stress and capillary pressure gradient are equally important or, in other words, to find the scale of the dynamic menisci. The small parameter 
$\varepsilon$ can be interpreted as the ratio between the length scale (in the direction perpendicular to the plate) at which viscous effect are important, and the capillary length. Requiring that the curvatures 
$\kappa _i$ are of order unity (so dimensional curvatures are of order 
$\ell _c^{-1}$) when the films approach the static menisci, we find 
$\alpha =1/2$ (Wilson Reference Wilson1982). Expanding (2.7) and (2.23)–(2.27) for 
$\varepsilon \ll 1$ and retaining terms up to order 
$\varepsilon ^{1/2}$ yields
\begin{equation} \frac{\partial}{\partial \hat{z}} \left[\frac{Ca}{\varepsilon^{3/2}} \hat{h}_i+ \left(\frac{\partial^3\hat{h}_1}{\partial\hat{z}^3} + \varSigma \frac{\partial^3\hat{h}_2}{\partial\hat{z}^3}-\varepsilon^{1/2}\right)\hat{F}_{i1} + \left(\varSigma \frac{\partial^3\hat{h}_2}{\partial\hat{z}^3} - \varepsilon^{1/2}R\right)\hat{F}_{i2}\right] = 0, \end{equation}
where 
$\hat {F}_{ij}$ are the functions defined in (2.24)–(2.27) but in terms of 
$\hat {h}_i$ instead of 
$h_i$. One assumption that we make here is that the dynamic meniscus region is located above 
$z = {\rm \Delta} H$. This means that the two interfaces 
$h_1(z)$ and 
$h_2(z)$ are defined there. Note that the curvatures are reduced to the second derivatives of the thicknesses, as the next term in their expansion is of order 
$\varepsilon$.
Analogously to the formulation developed in Wilson (Reference Wilson1982) for the one-liquid case, (2.28) encompasses the effects of both capillary- and gravity-driven drainage on the steady-state film thicknesses. In this work, we will focus on the case where capillary effects prevail over gravity. In (2.28), viscous stresses will be of the same order as capillary ones if 
$\varepsilon = Ca^{2/3}$, which is the same scaling obtained in the one-liquid case by Landau & Levich (Reference Landau and Levich1942). Finally, at leading order, the equations satisfied by the thicknesses are
\begin{equation} \frac{\partial}{\partial \hat{z}} \left[\hat{h}_i + \frac{\partial^3 \hat{h}_1}{\partial \hat{z}^3}\hat{F}_{i1} + \varSigma \frac{\partial^3 \hat{h}_2}{\partial \hat{z}^3}\left(\hat{F}_{i1} + \hat{F}_{i2}\right)\right] = 0.\end{equation}
The boundary conditions far from the liquid bath remain unchanged, as compared with the full formulation (developed in the last paragraph of 2.2): 
$\hat {h}_i^{\prime } = \hat {h}_i^{\prime \prime } = 0$ for 
$\hat {z} \rightarrow + \infty$, where the prime denotes the derivative with respect to 
$\hat {z}$. On the contrary, the boundary conditions near the reservoir have to be imposed through an asymptotic matching with the static menisci, as our ability to rigorously implement these boundary conditions at the liquid bath was lost when the curvature was linearised and gravity was neglected.
2.4. Asymptotic matching with the static menisci
2.4.1. Static configuration
For the purpose of asymptotic matching, we briefly study the static configuration depicted in figure 2. Setting 
$Ca = 0$ in (2.14) and (2.15) (which also implies 
$Q_i = 0$) leads to 
$\varPi _1 = \varPi _2 = 0$, or equivalently
$$\begin{gather} \frac{\partial\kappa_1}{\partial z} + R - 1 = 0, \end{gather}$$
$$\begin{gather}\varSigma \frac{\partial\kappa_2}{\partial z} - R = 0 . \end{gather}$$
As done in the one-liquid case (see for instance Landau & Lifshitz (Reference Landau and Lifshitz1987), chapter 7, page 243), these equations can be integrated to find the exact solutions for the shape of the static menisci 
$h_1(z)$ and 
$h_2(z)$. Since our purpose is to use these solutions to obtain asymptotic matching conditions for the dynamic thin film equations (2.29), we only need the approximated shape of the static interfaces (I) and (II) near the contact lines on the plate, i.e. close to positions 
$z_{cl,1}$ and 
$z_{cl,2}$, respectively.
Sketch of the static configuration: an immobile vertical plate is wetted by a compound bath at rest, made of a lighter liquid (2) on top of a denser one (1). The liquid 1/liquid 2 and liquid 2/air interfaces, denoted by (I) and (II), climb up to heights 
$z_{cl,1}$ and 
$z_{cl,2}$, respectively. The distance between the two corresponding contact lines on the plate is therefore 
${\rm \Delta} z_{cl} = z_{cl,2} - z_{cl,1}$.
Because 
$h_1 (z_{cl,1}) = h_2 (z_{cl,2}) = 0$ (by definition) and 
$h^{\prime }_1 (z_{cl,1}) = h^{\prime }_2 (z_{cl,2}) = 0$ (perfect wetting conditions on the plate), the Taylor expansions of the meniscus profiles in the vicinity of the contact lines read, for 
$z \le z_{cl,2}$:
$$\begin{gather} h_1 = \tfrac{1}{2}\kappa_1(z_{cl,1})\left(z - z_{cl,1}\right)^2, \end{gather}$$
$$\begin{gather}h_2 = \tfrac{1}{2}\kappa_2(z_{cl,2})\left(z - z_{cl,2}\right)^2. \end{gather}$$By solving (2.30) and (2.31), using the definition (2.6), we obtain the expressions for the curvatures at the contact lines,
$$\begin{gather} \kappa_1(z_{cl,1}) = \sqrt{2(1-R)} \quad \text{at}\ z_{cl,1} = \sqrt{2/(1-R)}, \end{gather}$$
$$\begin{gather}\kappa_2(z_{cl,2}) = \sqrt{2R/\varSigma} \quad \text{at}\ z_{cl,2} = {\rm \Delta} H + \sqrt{2\varSigma/R}, \end{gather}$$for interfaces (I) and (II) respectively.
2.4.2. Asymptotic matching
In the one-liquid case, thanks to the invariance of the problem in the 
$\hat {z}$-direction, imposing the value of the curvature in the limit 
$\hat {z} \rightarrow -\infty$ is enough to perform the asymptotic matching (Landau & Levich Reference Landau and Levich1942). In our two-liquid case, however, we benefit from the invariance along 
$\hat {z}$ only for one of the interfaces because the relative position of interfaces (I) and (II) far from the plate is fixed (they are separated by a given distance 
${\rm \Delta} H$). Consequently, the boundary conditions in the limit 
$\hat {z} \rightarrow -\infty$ not only result from imposing the curvatures of the two static menisci, but also require specifying how far apart the interfaces are. In order to fulfil these two conditions, the film thicknesses 
$\hat {h}_i$ are required to follow the parabolic approximations of the static menisci ((2.32) and (2.33)), expressed in terms of the rescaled variables 
$\hat {h}_i = Ca^{-2/3} h_i$ and 
$\hat {z} = Ca^{-1/3} z$
$$\begin{gather} \hat{h}_1 = \tfrac{1}{2}\sqrt{2\left(1-R\right)} \left(\hat{z} - \hat{z}_{cl,1}\right)^2, \end{gather}$$
$$\begin{gather}\hat{h}_2 = \tfrac{1}{2}\sqrt{2R/\varSigma}\left(\hat{z} - \hat{z}_{cl,2}\right)^2, \end{gather}$$
in the limit 
$\hat {z} \rightarrow -\infty$. We have introduced 
$\hat {z}_{cl,i} = z_{cl,i} \, Ca^{-1/3}$, where 
$z_{cl,i}$ denotes the dimensionless location where interfaces (I) and (II) (for 
$i=1$ and 
$2$ respectively) meet the plate in the static configuration and for perfect wetting conditions. Note that the matching between the dynamic and the static menisci will occur within a vertical distance of order 
$\ell _c \, Ca^{1/3}$ from the contact lines of the hydrostatic solutions. For 
$\varepsilon ^{1/2} = Ca^{1/3} \ll 1$ (condition to neglect gravity in the dynamic menisci, see §§ 2.3 and 2.6), this distance is much smaller than the capillary length, therefore justifying that the static menisci can be approximated by parabolas.
2.5. Dimensionless control parameters
The problem formulated above depends on a large number of dimensionless control parameters: 
$Ca$, 
$R$, 
$\varSigma$, 
$M$ and 
${\rm \Delta} H$. For the sake of conciseness, we will restrict our analysis to the parameters that are expected to have the largest impact on the flow structure and whose effect cannot be easily accounted for by some scaling argument.
We can estimate practically relevant orders of magnitude by considering a system made of a 
$40\,\%$ glycerol aqueous solution (liquid 1) with a floating layer of silicone oil (liquid 2). The corresponding densities are 
$\rho _1 = 1100$ kg m
$^{-3}$ for the glycerol solution (Takamura, Fischer & Morrow Reference Takamura, Fischer and Morrow2012) and 
$\rho _2 = 970$ kg m
$^{-3}$ for the oil, yielding a density ratio 
$R = 0.885$. The liquid 1/liquid 2 interfacial tension and liquid 2/air surface tension are 
$\sigma _{12} = 30$ and 
$\sigma _{2a} = 20$ mN m
$^{-1}$, respectively, yielding 
$\varSigma = 0.667$. Finally, the viscosity of the glycerol solution is 
$\mu _1 = 3.6$ mPa s (Takamura et al. Reference Takamura, Fischer and Morrow2012) and silicone oils can have a viscosity 
$\mu _2$ ranging from a fraction of a milliPascal second to hundreds of Pascal seconds, while keeping their other properties (density and surface tension) roughly constant. The values and ranges chosen in our computations for the dimensionless control parameters are summarised in table 1 and discussed in the following.
Main dimensionless control parameters of the problem and corresponding values or ranges explored in this work.
2.5.1. Capillary number 
$Ca$ and floating layer thickness ${\rm \Delta} H$
Importantly, after rescaling the problem as described in § 2.3, the capillary number 
$Ca$ disappears from the formulation (see (2.29)), except in the matching conditions (2.36) and (2.37). In the static menisci, with which the matching is performed in the limit 
$\hat {z} \rightarrow -\infty$, the dimensionless thickness 
${\rm \Delta} H$ of the floating layer remains related to the rescaled distance 
${\rm \Delta} \hat {z}_{cl}$ between the contact lines of liquids 1 and 2 (see § 2.4) through
\begin{equation} {\rm \Delta}\hat{z}_{cl}\,Ca^{1/3} = {\rm \Delta} H + \sqrt{\frac{2\varSigma}{R}} - \sqrt{\frac{2}{1-R}}. \end{equation}
Thus, for given fluid properties, changing the capillary number 
$Ca$ has the same effect as changing the thickness 
${\rm \Delta} H$. For this reason, throughout this work, we will vary only 
${\rm \Delta} H$ while keeping 
$Ca = 10^{-3}$ constant. This specific value was chosen in order to meet the condition 
$Ca^{1/3} \ll 1$, ensuring that gravity effects are negligible in the dynamic meniscus region, while being relevant in many applications (Rio & Boulogne Reference Rio and Boulogne2017). The values of 
${\rm \Delta} H$ explored in this work correspond to the range where solutions with two coated films were found to exist, i.e. between 
${\rm \Delta} H = 2.57$ and 
$5.20$ for our values of 
$Ca$, 
$\varSigma$ and 
$R$.
2.5.2. Density ratio $R$
The density ratio will be set to 
$R = 0.885$ (corresponding to the glycerol solution/silicone oil system described above) throughout this study. We choose to keep this parameter constant since, in practice, 
$R$ will most of the time be of this order with aqueous solution/oil systems.
2.5.3. Surface tension ratio $\varSigma$
Four values of the surface tension ratio will be explored: 
$\varSigma = 0.333$, 
$0.667$ (corresponding to the glycerol solution/silicone oil system described above), 
$1.0$ and 
$1.27$. These values are chosen around 
$1$, which is what can be achieved with most common fluids.
2.5.4. Viscosity ratio $M$
For glycerol solution/silicone oil systems, the viscosities of both phases can be tuned over several orders of magnitude by varying the glycerol concentration in the water phase and changing the average chain length in the oil phase. Consequently, the values of viscosity ratio explored will span several decades, from 
$M = 10^{-2}$ to about 
$10^1$, the exact upper value being limited by the existence of solutions with two coated films, as will be shown in § 4.1.
2.6. Validity conditions of the model
2.6.1. Stability of the hydrostatic configuration
For the static configuration of figure 2 to exist in practice, we need the stability of the floating layer to be warranted. So far, no assumption has been made on the (dimensional) spreading coefficient 
$S = \sigma _{1a} - (\sigma _{2a} + \sigma _{12})$ of liquid 2 on liquid 1. In the case where 
$S>0$, the hydrostatic configuration where liquid 2 forms a continuous layer in the bath will be stable regardless of the thickness of this layer, at least as long as long range molecular forces (e.g. van der Waals) can be neglected (Léger & Joanny Reference Léger and Joanny1992). This is the case for the 
$40\,\%$ glycerol aqueous solution/silicone oil system described above, for which 
$\sigma _{1a} = 70$ mN m
$^{-1}$ (Takamura et al. Reference Takamura, Fischer and Morrow2012), leading to 
$S = 20\ \textrm {mN}\ \textrm {m}^{-1} > 0$.
In the case where 
$S<0$, however, the floating layer will be stable only above a certain critical thickness 
${\rm \Delta} H_c$ (De Gennes et al. Reference De Gennes, Brochard-Wyart and Quéré2013). This critical thickness, made dimensionless with the capillary length 
$\ell _c$, can be written with our notations
\begin{equation} {\rm \Delta} H_c = \sqrt{\frac{-2\mathscr{S}}{R(1-R)}}, \end{equation}
where 
$\mathscr {S} = S/\sigma _{12} < 0$ is the dimensionless spreading coefficient. For 
${\rm \Delta} H < {\rm \Delta} H_c$, the floating layer is metastable and, if perturbed, will dewet (Brochard-Wyart, Martin & Redon Reference Brochard-Wyart, Martin and Redon1993).
2.6.2. Negligible gravity
We included gravity in our first, most general formulation developed in §§ 2.1 and 2.2. When performing the appropriate scaling in the dynamic menisci, retaining only leading-order terms (see § 2.3), we showed that the gravity term can be neglected provided 
$\varepsilon ^{1/2} \ll 1$, that is to say 
$Ca^{1/3} \ll 1$. Note that this condition, which is met for our value of 
$Ca = 10^{-3}$, is the same as in the one-liquid LLD problem (de Ryck & Quéré Reference de Ryck and Quéré1998).
2.6.3. Negligible inertia
Even at low capillary numbers, inertial effects may also arise when low-viscosity liquids are used. In our model, scaled as presented in § 2.3, the conditions for inertia to be negligible as compared with viscous and capillary terms read
$$\begin{gather} \text{in liquid 1}\quad We = Re \, Ca \ll 1, \end{gather}$$
$$\begin{gather}\text{in liquid 2}\quad We = Re \, Ca \ll \mathrm{min}(\varSigma / R; M/R), \end{gather}$$
where 
$We$ is the Weber number (Rio & Boulogne Reference Rio and Boulogne2017), 
$Re = \rho _1 U \ell _c/ \mu _1$ the Reynolds number (based on the properties of liquid 1) and 
$Ca$ the capillary number (as defined in § 2.2). Taking the above-mentioned 
$40\,\%$ glycerol aqueous solution as liquid 1 and 
$Ca = 10^{-3}$, we find 
$We \approx 4 \times 10^{-3}$, showing that conditions (2.40)–(2.41) are fulfilled in the whole range of parameters explored in our work (see table 1).
Additionally, in the general formulation where gravity is still present, neglecting inertial forces as compared with gravitational ones in both liquid films requires 
$Re \, Ca^{2/3} \ll 1$, that is to say 
$We \ll Ca^{1/3}$, which is fulfilled for our parameters. Note that the softer condition 
$Re \, Ca^{2/3} \sim 1$ (i.e. 
$We \sim Ca^{1/3}$) is sufficient in the dynamic menisci, as gravity effects are already small there.
3. Results: flow morphology and entrainment mechanism
In this section, we will present and discuss the typical shape of interfaces (I) and (II) (§ 3.1), as well as the flow structure in the dynamic menisci (§ 3.2), obtained by solving the model described in § 2 using the pseudo-transient numerical strategy described in Appendix B. Based on the observation of these results, we will propose a simplified, geometrical description (§ 3.3) allowing us to derive scaling laws for the asymptotic thickness of the coated liquid films (§ 3.4).
3.1. Shape of the interfaces
Figure 3(a) allows us to appreciate the typical shape of interfaces (I) and (II), as well as their matching to the static menisci (parameters 
$\varSigma = 0.667$, 
$R = 0.885$, 
$M = 1$, 
$Ca = 10^{-3}$ and 
${\rm \Delta} H = 3.404$). We focus on scales of the order of 
$Ca^{2/3}$ and 
$Ca^{1/3}$, on the horizontal and vertical axes, respectively. The regions displayed correspond to the dynamic menisci and thin liquid films dragged on the plate. Interfaces (I) (separating liquids 1 and 2) and (II) (between liquid 2 and air) are represented by blue and orange solid lines, respectively. These are the dynamic solutions obtained by solving (2.29), but plotted in their dimensionless non-rescaled form 
$h_i = \hat {h}_i Ca^{2/3}$ so that they may be represented on the same graph as the static menisci. The blue and orange dashed lines, obtained by solving the static equations (2.30) and (2.31) respectively, show the shapes of the static menisci expected when the plate is not moving and the liquids have zero contact angle with it (perfect wetting conditions). The location of the plate is outlined by the thick vertical grey line at 
$x = 0$. In addition to the static menisci themselves, their parabolic expansions, given by (2.32) and (2.33), are displayed with magenta and red dotted lines. Towards the bath, figure 3a shows that interfaces (I) and (II) (solid lines) approach these parabolic expansions, as required by the matching conditions (2.36) and (2.37). Using the common terminology in matched asymptotic expansions, this illustrates that the parabolas described by (2.32) and (2.33) are simultaneously the inner limits (i.e. close to the plate) of the outer solutions (static menisci) and the outer limits (approaching the bath) of the inner solutions (thin liquid films computed using lubrication theory).
Shape of the interfaces and matching to the static menisci for 
$\varSigma = 0.667$, 
$R = 0.885$, 
$M = 1$, 
$Ca = 10^{-3}$ and 
${\rm \Delta} H = 3.403$. (a) In the dynamic meniscus region, the interfaces (solid lines) depart from the static solutions (dashed lines) to connect to two thin films of uniform thicknesses (that can barely be distinguished at this scale). The dotted lines show the parabolic approximations of the static menisci near the plate, used in the matching conditions (2.36) and (2.37). (b) The horizontal distance 
$\delta h = h_2 - h_1$ between interfaces (I) and (II) exhibits a strong and localised decrease towards its asymptotic value 
$\delta h_{\infty }$ (black dashed line). (c) Similarly, the vertical distance 
${\rm \Delta} z$ between interfaces (I) and (II) also plummets, shortly before the two interfaces reach their asymptotic positions in the 
$x$-direction.
Further downstream from the menisci, interfaces (I) and (II) get very close to each other (until they can no longer be told apart at the scale of the plot) and flatten out to converge towards their asymptotic positions. Figures 3(b) and 3(c) allow us to better appreciate this convergence by looking at, on the one hand, the evolution of the horizontal distance 
$\delta h = h_2 - h_1$ between interfaces (I) and (II) with the vertical coordinate 
$z$ (panel b) and, on the other hand, the evolution of the vertical distance 
${\rm \Delta} z$ between these interfaces with the horizontal coordinate 
$x$ (panel c). Both panels reveal a very abrupt decay of the distance between interfaces (I) and (II) as they approach their asymptotic positions 
$x=h_{1,\infty }$ and 
$x=h_{2,\infty }$, respectively. Asymptotically, the uniform and steady thicknesses bounded by these interfaces are 
$h_{1,\infty }$ for the film of liquid 1 (referred to as lower film) and 
$\delta h_{\infty } = h_{2,\infty } - h_{1,\infty }$ for the film of liquid 2 (referred to as the upper film).
Remarkably, the thickness of the upper film is much smaller than that of the lower one: in the example of showed in figure 3, 
$\delta h_{\infty } = 1.97 \times 10^{-4}$ while 
$h_{1,\infty } = 1.33 \times 10^{-2}$ or, in terms of rescaled quantities, 
$\delta \hat {h}_{\infty } =\delta h_{\infty } \, Ca^{-2/3} = 1.97 \times 10^{-2}$ while 
$\hat {h}_{1,\infty } = h_{1,\infty } \, Ca^{-2/3} = 1.33$. Note that, with this normalisation, the rescaled thickness for a one-liquid Landau–Levich–Derjaguin film would simply be 
$\hat {h}_{\mathrm {LLD}} = 0.9458$, as the corresponding dimensional thickness is equal to 
$0.9458 \, \ell _c \, Ca^{2/3}$ (Landau & Levich Reference Landau and Levich1942). As we will see further on in § 4, this feature (
$\delta \hat {h}_{\infty } \ll \hat {h}_{1,\infty }$) is observed in all the parameter ranges explored in this work. We can already notice that differences in the properties of liquids 1 and 2 are not sufficient to account for such a discrepancy because the dimensional capillary lengths 
$\ell _{c,1}$ and 
$\ell _{c,2}$, corresponding to interfaces (I) and (II) respectively, remain of the same order of magnitude: 
$\ell _{c,2}/\ell _{c,1} = \sqrt {\varSigma (1/R - 1)} = 0.208 - 0.406$ for all parameters considered in our study (see table 1). As we will develop in next paragraphs, the difference in the order of magnitude of 
$h_{1,\infty }$ and 
$\delta h_{\infty }$ is connected to the very different flow patterns conveying liquid into the lower and upper films.
3.2. Flow structure
In this paragraph, we turn our attention to the structure of the flow in the dynamic meniscus region. As an example, figure 4 displays various relevant flow magnitudes, obtained for parameters 
$\varSigma = 0.667$, 
$R = 0.885$, 
$M = 1$, 
$Ca = 10^{-3}$ and 
${\rm \Delta} H = 3.404$, as functions of the vertical coordinate 
$\hat {z}$. In figure 4(a), we plot the streamlines in the dynamic menisci for liquid phases 1 (in blue) and 2 (in orange). The corresponding analytical expressions for the velocity fields are given in Appendix A by (A1) and (A2) for liquid 1 and (A3) and (A4) for liquid 2, respectively. While, as expected, the plate drags the lower liquid up, streamlines reveal that the flow in the vicinity of interface (I) (liquid 1/liquid 2 interface, thick blue line) is actually going downwards upstream the stagnation point (blue dot). The consequences of this flow pattern on the coated liquid films can be better understood by looking at the main physical effects at play: viscous entrainment by shear stresses at the dragging interfaces (plate/liquid 1 for the lower phase, liquid 1/liquid 2 for the upper phase) and capillary suction generated by corresponding interfacial curvatures. The former is presented in figure 4(b), where the dimensionless rescaled shear stress at different interfaces of interest is plotted as a function of the vertical coordinate 
$\hat {z}$, while the latter is quantified in figure 4(c), displaying the dimensionless rescaled pressure gradient as a function of 
$\hat {z}$.
Flow structure for dimensionless parameters 
$\varSigma = 0.667$, 
$R = 0.885$, 
$M = 1$, 
$Ca = 10^{-3}$ and 
${\rm \Delta} H = 3.403$ as a function of the vertical coordinate 
$\hat {z}$. (a) Streamlines in liquid 1 (blue) and liquid 2 (orange). (b) Shear stresses at the plate/liquid 1 interface, 
$\hat {\tau }_{01}$ (black solid line), and at the liquid 1/liquid 2 interface, 
$\hat {\tau }_{12}$ (blue solid line). (c) Pressure gradients in liquid 1, 
$\hat {\Pi }_1$ (solid blue line), and in liquid 2, 
$\hat {\Pi }_2$ (solid orange line). For comparison, the red dashed lines represent the corresponding magnitudes in the one-liquid LLD theory: shear stress at the plate/liquid interface, 
$\hat {\tau }_{LLD}$ (b), and pressure gradient in the liquid, 
$\hat {\Pi }_{LLD}$ (c).
Let us first focus on the viscous stresses, which promote liquid film entrainment. The solid black line in figure 4(b) represents the shear stress 
$\hat {\tau }_{01}$, defined by (A9), at the interface between the plate and liquid 1 in our two-liquid configuration. For comparison, the dashed red line shows the shear stress 
$\hat {\tau }_{LLD}$ at the plate/liquid interface in the one-liquid film coating problem; the corresponding expression is given by (A11a,b). Both curves exhibit qualitatively the same shape: the shear stress at the plate increases progressively as the height above the surface of the bath increases, goes through a maximum and then decays back to zero as 
$\hat {z}$ keeps increasing. This is in sharp contrast with the trend followed by the shear stress 
$\hat {\tau }_{12}$ at interface (I), separating liquids 1 and 2 (solid blue line, defined by (A10)). The shear stress at interface (I) is found to be essentially zero everywhere, except for a small peak in a narrow region (around 
$\hat {z} \approx 24$ in this example), which coincides with the zone where interfaces (I) and (II) get very close to each other (see figure 3 and § 3.1). This observation is consistent with the streamlines of figure 4(a) that show that liquid 2 is essentially recirculating on top of liquid 1, as no shear is transmitted along most of interface (I). Note that the shear stress at interface (II) (separating liquid 2 and the atmosphere) is identically zero, as set by the boundary condition (2.19).
We now turn to the capillary pressure gradients, which impede thin liquid film entrainment. Figure 4(c) displays the dimensionless rescaled pressure gradient 
$\hat {\Pi }_1$ (blue) and 
$\hat {\Pi }_2$ (orange) in liquid phases 1 and 2, defined by (A5) and (A6), respectively. Again, the two liquid phases exhibit very different behaviours: the lower liquid is subjected to a mild pressure gradient, spread over distances of several units in 
$\hat {z}$. On the contrary, there is virtually no pressure gradient in the upper liquid, except for a high peak around 
$\hat {z} \approx 24$, which also corresponds to the region of highest shear stress (in absolute value) at the liquid 1/liquid 2 interface.
3.3. Geometrical approach: virtual contact point
Let us focus on the region where we observe the peak in the shear stress at the liquid/liquid interface, 
$\hat {\tau }_{12}$, and in the pressure gradient in the upper film, 
$\hat {\Pi }_2$. Figures 4(b) and 4(c) show that this region is very narrow in the vertical direction, as compared with the overall extension of the dynamic menisci. Moreover, figure 3(b) reveals that, concomitantly, the distance between interfaces (I) and (II) decays quickly down to the asymptotic one, 
$\delta h_{\infty }$. For these reasons, we find meaningful to model this region as a point, called the virtual contact point, where interfaces (I) and (II) virtually meet. Note that, in general, this point does not coincide with a stagnation point at interface (I), as can be seen in the example of figure 4. The vertical coordinate of the virtual contact point is denoted by 
$z_{\ast }$. In the following, all quantities evaluated at this point will be marked by an asterisk.
Introducing the virtual contact point allows us to develop a geometrical, asymptotic (‘zoomed-out’) description of the flow, as sketched in figure 5. By construction, we have 
$\delta h \ll h_{1}$, and asymptotically 
$\delta h_{\infty } \ll h_{1,\infty }$, downstream of the virtual contact point. We therefore simplify the geometry in this region, assuming that interfaces (I) and (II) are merged into a single interface (III) with effective dimensionless surface tension 
$1 + \varSigma$. In this framework, the different interfaces pictured in figure 5 can be described as follows.
(i) Interface (I): in the region below the virtual contact point, the very small shear stress acting on interface (I) (see figure 4b) implies that this surface can be regarded as shear free. As a consequence, the lower liquid 1 bounded by this interface behaves as a one-liquid Landau–Levich film.
(ii) Interface (II): the upper liquid 2 is bounded by the approximately shear-free interface (I) and the rigorously shear-free interface (II). Consequently, the pressure gradient inside this liquid must be zero, which is consistent with figure 4(c).
(iii) Interface (III): in the region above the virtual contact point, by construction, the dynamics of the effective interface (III) also obeys the one-liquid Landau–Levich equation.
Sketch showing the virtual contact point, where interfaces (I) and (II) are assumed to come into contact. Above this point, located at height 
$z^{\ast }$, interfaces (I) and (II) merge into a single interface, denoted by (III), with an effective surface tension equal to the sum of that of interfaces (I) and (II).
This simplified description will be useful to explain some of the most salient features of the two-liquid flow using well-known properties of its one-liquid counterpart. More specifically, we will show that the asymptotic film thicknesses can be rationalised, and even quantitatively predicted to some extent, looking at the flow variables at the virtual contact point.
3.4. Scaling laws for the film thicknesses
The results presented in figure 4 show that the two main competing forces – viscous stresses and capillary pressure gradient – reach their extreme values in the vicinity or at the virtual contact point. This suggests that the amount of liquid entrained in each film may be rationalised by considering solely the local shear stresses around that location. In what follows, we use scaling arguments – in the spirit of the ones developed, for instance, in Champougny et al. (Reference Champougny, Scheid, Restagno, Vermant and Rio2015) – to relate the (dimensionless rescaled) shear stresses at the virtual contact point to the steady-state (dimensionless rescaled) thicknesses of the coated liquid films, namely 
$\hat {h}_{1,\infty }$ for liquid 1 and 
$\delta \hat {h}_{\infty }$ for liquid 2.
3.4.1. Lower film thickness $\hat {h}_{1,\infty }$
In the case of the lower film, the shear stress responsible for liquid entrainment is the one at the plate/liquid 1 interface (black solid line in figure 4b). In dimensional terms, the maximum shear stress must be of the order of the viscosity times the plate velocity divided by the minimum thickness, i.e. that of the film
\begin{equation} -\tau_{01, \ast} \sim \frac{\mu_1 U}{h_{1,\ast}}. \end{equation}
In this expression, the star indicates that the quantities are evaluated at the virtual contact point. The dimensional thickness and stress are related to their dimensionless rescaled counterparts through 
$h_{1,\ast } = \hat {h}_{1,\ast } \, \ell _c \, Ca^{2/3}$ and 
$\tau _{01, \ast } = \hat {\tau }_{01,\ast }\ (\sigma _{12}/\ell _c) \, Ca^{1/3}$, respectively, yielding 
$-\hat {\tau }_{01,\ast } \hat {h}_{1,\ast } \sim 1$. Making the approximation that the film thickness has already reached its asymptotic value at the virtual contact point, namely that 
$\hat {h}_{1,\infty } \approx \hat {h}_{1,\ast }$, we arrive at
\begin{equation} \hat{h}_{1,\infty} \sim{-}\hat{\tau}_{01,\ast}^{{-}1}. \end{equation}
In figure 6(a), we plot the lower film thickness 
$\hat {h}_{1,\infty }$ as a function of the inverse of the shear stress at the virtual contact point, 
$-\hat {\tau }_{01,\ast }^{-1}$, for a variety of 
$\varSigma$, 
$M$ and 
${\rm \Delta} H$. All numerical data are found to collapse on a single master curve, which is convincingly adjusted by a linear fit (solid black line) with a prefactor equal to 
$0.67$, therefore supporting the scaling law (3.2). Note that the presence of a second lighter fluid always causes the lower film thickness 
$\hat {h}_{1,\infty }$ to be larger than that of a one-liquid Landau–Levich film 
$\hat {h}_{LLD} = 0.9458$.
Comparison of the asymptotic thicknesses 
$\hat {h}_{1,\infty }$ and 
$\delta \hat {h}_{\infty }$ obtained numerically with the predictions of scaling laws derived from the virtual contact point model for (a) the lower film (3.2) and (b) the upper film (3.6). The solid lines are best linear fits, yielding prefactors 
$0.67$ and 
$1.29$ for the lower and upper films, respectively. The data presented correspond to various values of 
$\varSigma$ (symbol shape), 
$M$ (greyscale) and 
${\rm \Delta} H$ (not marked).
3.4.2. Upper film thickness $\delta \hat {h}_{\infty }$
In the case of the upper film, the shear stress responsible for liquid entrainment is the one at interface (I), i.e. at the liquid 1/liquid 2 interface (blue solid line in figure 4b). Unlike the one at the plate, the velocity at interface (I) is not given a priori but depends in a non-trivial way on the parameters of the problem. For this reason, we cannot start from the definition of the shear stress, as we did in the case of the lower film. Instead, we write that, around the virtual contact point, the total shear force at interface (I) must balance the capillary suction exerted by interface (II). Denoting by 
$\ell$ the streamwise extension of the virtual contact point region (typically the width of the peaks in 
$\hat {\tau }_{12}$ and 
$\hat {\varPi }_2$ in figure 4), this force balance reads in dimensional terms
\begin{equation} -\tau_{12}^{{\ast}} \ell \sim \sigma_{2a} \frac{\delta h_{{\ast}}}{\ell^2} \delta h_{{\ast}}. \end{equation}Additionally, matching the curvatures of the dynamic and static menisci of liquid 2 imposes (still using dimensional quantities)
\begin{equation} \frac{\delta h_{{\ast}}}{\ell^2} \sim \frac{1}{\ell_{c,2}}, \end{equation}
where we recall that 
$\ell _{c,2} = \sqrt {\sigma _{2a}/\rho _2 g}$ is the capillary length related to interface (II). This condition is analogous to the asymptotic matching condition with the static meniscus introduced by Landau & Levich (Reference Landau and Levich1942) in the one-liquid dip-coating problem. Combining (3.3) and (3.4), we arrive at the dimensional scaling expression
\begin{equation} \delta h_{\infty} \sim \ell_{c,2} \left(-\tau_{12,\ast} \frac{ \ell_{c,2}}{\sigma_{2a}} \right)^2, \end{equation}
where we approximated 
$\delta h_{\infty } \approx \delta h_{\ast }$. Using the definitions 
$\delta h_{\infty } = \delta \hat {h}_{\infty } \ell _c Ca^{2/3}$ and 
$\tau _{12, \ast } = \hat {\tau }_{12,\ast } (\sigma _{12}/\ell _c) \, Ca^{1/3}$, we finally express (3.5) in terms of the rescaled dimensionless variables
\begin{equation} \delta\hat{h}_{\infty} \sim \left(-\hat{\tau}_{12,\ast} \right)^2 R^{{-}3/2} \varSigma^{{-}1/2}. \end{equation} In figure 6(b), we plot the upper film thickness 
$\delta \hat {h}_{\infty }$ as a function of the right-hand term of (3.6) for a variety of 
$\varSigma$, 
$M$ and 
${\rm \Delta} H$. When plotted in this way, the numerical data are found to collapse on a reasonably linear master curve. The linear fit (solid black line) yields a prefactor equal to 
$1.29$, showing that (3.6) can be used to estimate the upper film thickness 
$\delta \hat {h}_{\infty }$.
To conclude, perhaps the most important lesson learnt from scalings (3.2) and (3.6) is the universality of the entrainment mechanism. Given the values of the shear stresses at the virtual contact point (
$\hat {\tau }_{01,\ast }$ for liquid 1, 
$\hat {\tau }_{12,\ast }$ for liquid 2), the amount of fluid dragged in each film can be readily estimated using the same ideas exposed in the original work of Landau & Levich (Reference Landau and Levich1942).
4. Results: parametric dependence of film thicknesses
In this section, we turn our attention to the parametric dependence of the asymptotic thicknesses of the coated liquid films, 
$\hat {h}_{1,\infty }$ for the lower film and 
$\delta \hat {h}_{\infty }$ for the upper film. The main control parameters varied in our study are the dimensionless floating layer thickness 
${\rm \Delta} H$, the viscosity ratio 
$M$ and the surface tension ratio 
$\varSigma$ (see § 2.5). Importantly, our results reveal that double coating solutions only exist in finite areas of the parameter space, which we dub existence islands (§ 4.1). We will propose arguments explaining some trends and boundaries observed for the film thicknesses as a function of 
${\rm \Delta} H$ (§§ 4.2 and 4.3).
4.1. Thickness maps in the $M - {\rm \Delta} H$ parameter space
In figure 7, the asymptotic film thicknesses are presented as colour maps in the 
$M - {\rm \Delta} H$ parameter space for four different values of the surface tension ratio 
$\varSigma$. Panels (a,c,e,g) show the lower film thickness 
$\hat {h}_{1,\infty }$ while (b,d,f,h) displays the upper film thickness 
$\delta \hat {h}_{\infty }$. For an easier quantitative reading, cuts along the 
$M$ and 
${\rm \Delta} H$ directions for the surface tension ratio 
$\varSigma = 1.27$ are presented in figure 8. Panels (a,c) show the variation of the lower film thickness 
$\hat {h}_{1,\infty }$ with 
$M$ for fixed values of 
${\rm \Delta} H$ (figure 8a) and with 
${\rm \Delta} H$ for fixed values of 
$M$ (figure 8c). Similarly, (b,d) display the variation of the upper film thickness 
$\delta \hat {h}_{\infty }$ with the same parameters 
$M$ (figure 8b) and 
${\rm \Delta} H$ (figure 8d).
Asymptotic thicknesses 
$\hat {h}_{1,\infty }$ (lower film, a,c,e,g) and 
$\delta \hat {h}_{\infty }$ (upper film, b,d,f,h), shown as colour contour maps in the (
$M$, 
${\rm \Delta} H$) parameter space. Each row corresponds to a different surface tension ratio 
$\varSigma$. The red dashed lines enclose the ‘existence islands’ of a double coating, namely the areas in the parameter space where 
$\delta \hat {h}_{\infty } \ge 5 \times 10^{-4}$. The values of 
$\hat {h}_{1,\infty }$ obtained outside these islands are shaded.
Cuts of the thickness maps presented in figures 7(g) and 7(h) (
$\varSigma = 1.27$), within the existence island enclosed in the dashed red contour. The top panel shows the variation of (a) the lower film thickness 
$\hat {h}_{1,\infty }$ and (b) the upper film thicknesses 
$\delta \hat {h}_{\infty }$ with the viscosity ratio 
$M$, for constant values of the floating layer thickness 
${\rm \Delta} H$. The bottom panel shows the variation of the same quantities (on (c) and (d), respectively) with the floating layer thickness 
${\rm \Delta} H$, for constant values of the viscosity ratio 
$M$.
Together, figures 7 and 8 reveal dramatic qualitative and quantitative differences between the lower and upper coated films. As already observed and rationalised in § 3, the lower film exhibits final thicknesses 
$\hat {h}_{1,\infty }$ of order unity (as expected from the Landau–Levich scaling, see § 3.4), while the upper film reaches steady-state thicknesses 
$\delta \hat {h}_{\infty }$ that are about 
$10^{-3}$ to 
$10^{-2}$ times smaller. Not only are the lower and upper film thicknesses very disparate, but they also depend very differently on the control parameters. Figures 7(a), 7(c), 7(e) and 7(g) reveal that the thickness of the lower film weakly depends on the viscosity ratio 
$M$ and grows monotonically with the depth of the floating layer 
${\rm \Delta} H$. On the contrary, figures 7(b), 7(d), 7(f) and 7(h) show that the thickness 
$\delta \hat {h}_{\infty }$ of the upper film depends non-monotonically on both 
$M$ and 
${\rm \Delta} H$, exhibiting a maximum whose exact position in the (
$M, {\rm \Delta} H$) parameter space depends on the surface tension ratio 
$\varSigma$.
Around this maximum, the upper film thickness decreases in all directions, eventually going down to values reaching our numerical accuracy (
$\delta \hat {h}_{\infty } \approx 5 \times 10^{-4}$, dashed red lines in figure 7). We observed that increasing the resolution of our numerical scheme barely affects the position of the red contours, leading us to conclude that there is no solution with two entrained films beyond these limits (black areas in figure 7b,d,f,h). The zones of the (
$M, {\rm \Delta} H$) parameter space enclosed within the red dashed contours will therefore be referred to as existence islands for the double-layer coating. Note that values of the lower film thickness 
$\hat {h}_{1,\infty }$ can still be obtained outside these islands (see shaded areas in figure 7a,c,e,g). However, since the theory used to obtain them postulates the presence of two entrained liquid films, only the data enclosed in the existence islands – where a double coating solution exists – can be discussed.
4.2. Minimum and maximum ${\rm \Delta} H$ for the existence of two films
The very different scales on the 
$M$ and 
${\rm \Delta} H$ axes in figure 7 highlight that the existence islands extend only along a finite range of values of 
${\rm \Delta} H$, getting narrower and narrower as 
$\varSigma$ increases, while they span several orders of magnitude in 
$M$ values. In this paragraph, we aim at providing some physical arguments to rationalise this observation.
The existence islands exhibit a relatively well-defined lower boundary in 
${\rm \Delta} H$, which becomes independent of the viscosity ratio 
$M$ for 
$M \ll 1$. This limit can be understood from hydrostatic considerations. A necessary condition to have two entrained films is that the static menisci of the two liquids touch the plate, as depicted in figure 2. In other words, the apparent contact line of the upper interface should lie above that of the lower one. For this to occur, the distance between the apparent contact lines of the lower and the upper menisci, 
${\rm \Delta} z_{cl}$, must be positive. Combining (2.34) and (2.35) this condition translates into
\begin{equation} {\rm \Delta} H \ge {\rm \Delta} H_{\mathrm{min}} = \sqrt{\frac{2}{1-R}} - \sqrt{\frac{2\varSigma}{R}}. \end{equation}
The values of 
${\rm \Delta} H_{\mathrm {min}}$ predicted by (4.1) for surface tension ratios 
$\varSigma =$ 
$0.33$, 
$0.67$, 
$1.00$ and 
$1.27$ are presented in table 2 and displayed in figures 7(b), 7(d), 7(f) and 7(h) as dashed white lines. Comparison to the numerical data reveals a good agreement with the lower limit of the existence islands for 
$M \ll 1$.
Comparison between some features of the results presented in figure 7 and the corresponding values predicted using simplified approaches (§ 5), for different values of the surface tension ratio 
$\varSigma$. The quantities extracted from figure 7 are evaluated in the limit 
$M \ll 1$ and within the existence islands (red dashed contours), where a double coating solution exists. These quantities are (a) the lower limit in 
${\rm \Delta} H$ of the existence islands and (b) the minimum values of the asymptotic lower film thickness 
$\hat {h}_{1,\infty }$.
Regarding the maximum 
${\rm \Delta} H$ for the existence of a double-film configuration, 
${\rm \Delta} H_{{max}}$, the physical origin is different. For the upper film to be dragged, shear must be transmitted from the plate to the interface between the two liquids. However, the region where the plate exerts shear on the lower film is limited to an extension of order 
${\rm \Delta} \hat {z} \sim 10$ in the streamwise direction, as can be seen in figure 4(b). This means that if the virtual contact point – where the two interfaces come close to each other – is outside this limited region in 
$\hat {z}$, no shear can be transferred to the liquid/liquid interface, and thus the upper film cannot be dragged. In terms of 
${\rm \Delta} H$, (2.38) allows us to translate the extension 
${\rm \Delta} \hat {z} \sim 10$ into
\begin{equation} {\rm \Delta} H_{{max}} - {\rm \Delta} H_{{min}} \sim 10 \, Ca^{1/3}. \end{equation}
The condition 
$Ca^{1/3} \ll 1$ (see §§ 2.3 and 2.6) explains why we expect the existence islands to span only along a reduced extension in 
${\rm \Delta} H$. In our particular case where 
$Ca = 10^{-3}$, (4.2) predicts island extensions of the order of unity in 
${\rm \Delta} H$, which is compatible with the data displayed in figure 7.
4.3. Effect of ${\rm \Delta} H$ on the lower film thickness $\hat {h}_{1,\infty }$
In this paragraph we use the equivalent description presented in § 3.3 to explain the effect of the floating layer thickness 
${\rm \Delta} H$ on the asymptotic thickness 
$\hat {h}_{1,\infty } \approx \hat {h}_{III,\infty }$ of the lower film (see figure 5). To do so, we first examine the curvatures of the different interfaces represented in figure 5 and make the following observations.
(i) Interface (I): since interface (I) behaves as in the one-liquid LLD theory, its curvature
$\hat {h}''_{I}$ decays monotonically with $\hat {z}$ (see Appendix C). The maximum curvature, $\sqrt {2(1-R)}$, is found in the limit $\hat {z} \rightarrow -\infty$, corresponding to the lower static meniscus.(ii) Interface (II): since the pressure gradient in liquid 2 is approximately zero, interface (II) has a constant curvature, given by that of the upper static meniscus:
$\hat {h}''_{II} =\sqrt {2R/\varSigma }$.(iii) Interface (III): for the same reason as interface (I), the curvature of interface (III) decreases monotonically with
$\hat {z}$, reaching its maximum value $\hat {h}''_{III, \ast }$ at its lowest point, i.e. the virtual contact point.
Taking again advantage of the fact that interface (III) can be described by the one-liquid LLD theory, we approximate (see (C5))
\begin{equation} \hat{h}_{III,\infty} \approx \frac{1.336}{(1+\varSigma)^{2/3} \hat{h}''_{III, \ast}}. \end{equation}
Moreover, the condition that the pressure inside liquid 1 must be continuous across the virtual contact point determines the value of the curvature of interface (III) at this location. On the one hand, the pressure slightly downstream the virtual contact point is given by the capillary pressure jump across interface (III), 
$(1 + \varSigma )\ \hat {h}''_{III,*}$. On the other hand, just upstream that point, the pressure is equal to the sum of the capillary pressure jumps across interfaces (II) and (I), which are respectively 
$\varSigma \hat {h}''_{II,*}$ and 
$\hat {h}''_{I,*}$. Requiring pressure continuity at the virtual contact point leads to
\begin{equation} \hat{h}''_{III, \ast} = \frac{\hat{h}''_{I, \ast} + \varSigma \hat{h}''_{II}}{1+\varSigma}. \end{equation}Combining (4.3) and (4.4), we arrive at the following expression, relating the asymptotic lower film thickness to the curvature of interface (I) at the virtual contact point
\begin{equation} \hat{h}_{1,\infty} \approx \hat{h}_{III,\infty} \approx \frac{1.336\ (1+\varSigma)^{1/3}}{\hat{h}''_{I, \ast} + \sqrt{2\varSigma R}}. \end{equation}4.3.1. Monotonic behaviour of $\hat {h}_{1,\infty }$ with ${\rm \Delta} H$
This simplified view allows us to explain why 
$\hat {h}_{1,\infty }$ grows monotonically with the thickness of the floating layer, 
${\rm \Delta} H$. As this parameter grows, the virtual contact point, where interfaces (I) and (II) meet, displaces up downstream or, in other words, 
$\hat {z}^{\ast }$ increases. Since 
$\hat {h}''_{I}$ is a decreasing function of 
$\hat {z}$, the larger 
$\hat {z}^{\ast }$, the lower the corresponding curvature 
$\hat {h}''_{I, \ast }$ of interface (I). We therefore deduce that 
$\hat {h}''_{I, \ast }$ decays with 
${\rm \Delta} H$. Finally, (4.5) allows us to conclude that 
$\hat {h}_{1,\infty }$ must be a monotonically increasing function of 
${\rm \Delta} H$, as observed in our numerical results in figures 7 and 8.
4.3.2. Lower bound for $\hat {h}_{1,\infty }$
Elaborating on these ideas, we can also provide a lower bound for the asymptotic thickness of the lower film. If the virtual contact point is displaced far upstream (
$\hat {z} \rightarrow -\infty$, i.e. closer to the liquid bath), it will eventually reach the region where the curvature of interface (I) has its asymptotic (maximum) value: 
$\hat {h}''_{I, \ast } = \sqrt {2(1-R)}$. Substituting this value in (4.5) yields a lower bound for the lower film thickness
\begin{equation} \mathrm{inf}\ (\hat{h}_{1,\infty}) \approx 0.9458 \frac{\left(1+\varSigma\right)^{1/3}}{\sqrt{1-R} + \sqrt{\varSigma R}}. \end{equation}
We can compare the predictions of (4.6) to the values of 
$\hat {h}_{1,\infty }$ observed in figure 7. To do so, we should restrict ourselves to the area of the parameter space enclosed in the dashed red line (where solutions for two entrained films exist) and to the limit 
$M \ll 1$, in which the simplified geometrical model is valid. As shown in table 2, the minimum values 
$\mathrm {min}\ (\hat {h}_{1,\infty })$ obtained from the simulations are not only compatible with the lower bounds provided by (4.6) but also follow a similar trend with 
$\varSigma$.
5. Discussion
In this last section, we discuss some limitations of our hydrodynamic model in relation to practically relevant effects, such as partial wetting conditions between the two liquids and long range intermolecular forces across the upper film, when it is sufficiently thin.
5.1. Partial wetting conditions
As mentioned in § 2.6, the stability of a floating layer of liquid 2 on liquid 1 depends on the dimensionless spreading parameter 
$\mathscr {S} = S/\sigma _{12} = \varSigma ^{\prime } - (\varSigma + 1)$, where we have introduced 
$\varSigma ^{\prime } = \sigma _{1a}/\sigma _{12}$. For 
$\mathscr {S} > 0$, the floating layer is stable regardless of its thickness while, for 
$\mathscr {S} < 0$, only floating layers of dimensionless thickness 
${\rm \Delta} H > {\rm \Delta} H_c$ given by (2.39) are stable (De Gennes et al. Reference De Gennes, Brochard-Wyart and Quéré2013). In § 4.2, we estimated the minimum floating layer thickness 
${\rm \Delta} H_{\mathrm {min}}$ needed for a double liquid film to be entrained (4.1). In partial wetting conditions, our hydrodynamic description of the configuration of figure 1 is therefore warranted as long as 
${\rm \Delta} H_{\mathrm {min}} > {\rm \Delta} H_c$. This translates into the following condition on the dimensionless spreading parameter
\begin{equation} -\left(\sqrt{R} - \sqrt{\varSigma (1-R)} \right)^2 < \mathscr{S} < 0, \end{equation}
or, equivalently, on the dimensionless liquid 1/air surface tension 
$\varSigma ^{\prime }$
\begin{equation} 1 + \varSigma - \left( \sqrt{R} - \sqrt{\varSigma (1-R)} \right)^2 < \varSigma^{\prime} < 1 + \varSigma. \end{equation}5.2. Thinness of the upper film and disjoining pressure effects
As can be seen in figure 7, the dimensionless rescaled thickness of the upper film, 
$\delta \hat {h}_{\infty }$, is of the order of 
$10^{-3}$ to 
$10^{-2}$ in the range of parameters explored. The corresponding dimensional thickness, 
$\delta h_{\infty } = \ell _c Ca^{2/3} \times \delta \hat {h}_{\infty }$, is therefore expected to be in the range 20–200 nm for the typical values 
$Ca = 10^{-3}$ and 
$\ell _c \approx 2$ mm. Albeit small, these thicknesses are withing reach of techniques such as reflection interference contrast microscopy, as exemplified by the recent measurements of Kreder et al. (Reference Kreder, Daniel, Tetreault, Cao, Lemaire, Timonen and Aizenberg2018) on the oil layer wrapping water droplets advancing on lubricant-infused surfaces.
Given the thinness of the upper liquid film, one may wonder in which conditions intermolecular long-range forces would be expected to have a noticeable effect on the upper film dynamics. These interactions in the film are usually quantified by the disjoining pressure isotherm (Derjaguin & Churaev Reference Derjaguin and Churaev1978), which measures the relative force acting between its two interfaces as a function of their separation. For asymmetric films made of a pure liquid, such as our upper film of liquid 2, the disjoining pressure isotherm only contains van der Waals interactions, whose nature (attractive or repulsive) depends on the sign of the Hamaker constant 
$A$ (Israelachvili Reference Israelachvili2011). Following the approach of Champougny et al. (Reference Champougny, Rio, Restagno and Scheid2017), the effect of van der Waals interactions can be included in our formulation by adding a disjoining pressure gradient 
$\hat {\varPi }^{d}$ to the capillary pressure gradients 
$\hat {\varPi }_1$ and 
$\hat {\varPi }_2$ defined in Appendix A. Considering only non-retarded van der Waals forces (Léger & Joanny Reference Léger and Joanny1992) and using the rescaling exposed in §§ 2.2 and 2.3, the disjoining pressure gradient reads
\begin{equation} \hat{\varPi}^{d} = Ca^{{-}2} \mathcal{A} \frac{\hat{h}'_2 - \hat{h}'_1}{(\hat{h}_2-\hat{h}_1)^4}, \end{equation}
where 
$\mathcal {A} = A / 2{\rm \pi} \sigma _{12} \ell _c^2$ is a dimensionless Hamaker constant. Using the values of 
$Ca$, 
$\sigma _{12}$ and 
$\ell _c$ reported in § 2.5, as well as the typical dimensional Hamaker constant 
$|A| \approx 10^{-20}$ J for water–oil–air systems (Israelachvili Reference Israelachvili2011), we estimate 
$\mathcal {A} \approx 1.6 \times 10^{-14}$.
The upper film thickness 
$\delta \hat {h}_{\infty }^{crit}$ at which the disjoining pressure gradient becomes of the order of the capillary one can be estimated by requiring
\begin{equation} \hat{\varPi}_2 ={-}\varSigma \hat{h}'''_2 \sim \hat{\varPi}^{d} = Ca^{{-}2} \mathcal{A} \frac{\hat{h}'_2 - \hat{h}'_1}{(\hat{h}_2-\hat{h}_1)^4}. \end{equation}
In the rescaled variables used here, 
$\hat {h}_{2,\infty } \sim \hat {h}_{1,\infty }$ is of order unity, and so is the characteristic distance 
${\rm \Delta} \hat {z}$ over which the thicknesses vary. Thus, we have 
$\hat {h}'''_2 \sim 1$ and 
$\hat {h}'_2-\hat {h}'_1 \sim \hat {h}_2-\hat {h}_1 = \delta \hat {h}_{\infty }$, yielding
\begin{equation} \delta\hat{h}_{\infty}^{crit} \sim \left(|\mathcal{A}| \, Ca^{{-}2} \, \varSigma^{{-}1}\right)^{1/3}. \end{equation}
With the values of 
$\mathcal {A}$ and 
$Ca$ given above, we arrive at thresholds 
$\delta \hat {h}_{\infty }^{crit} \sim 3.6 \times 10^{-3}$, 
$2.9 \times 10^{-3}$, 
$2.5 \times 10^{-3}$ and 
$2.3 \times 10^{-3}$ for 
$\varSigma = 0.333$, 
$0.667$, 
$1.00$ and 
$1.27$, respectively. For each 
$\varSigma$, this means that the parts of the existence islands where 
$\delta \hat {h}_{\infty } \lesssim \delta \hat {h}_{\infty }^{crit}$ (essentially the edges, see figure 7) may be affected by disjoining pressure. Conversely, disjoining pressure effects are not expected to play a significant role in the interior of the existence islands. Note that the above reasoning allows us to estimate when van der Waals interactions are expected to influence the value of the asymptotic thickness of the upper film. A detailed discussion of the effect of long-range intermolecular forces on the stability of such a film (Fisher & Golovin Reference Fisher and Golovin2005; Bandyopadhyay & Sharma Reference Bandyopadhyay and Sharma2006) lies beyond the scope of the present work.
6. Conclusions
In this work, we investigated the dip-coating flow generated by a plate lifted at constant speed through a stratified bath made of two immiscible liquids, focusing on the case when both of them are entrained on the plate. We presented first a general formulation within the framework of the lubrication approximation that includes the full (nonlinear) expression for the curvature and the effect of gravity. Following a similar scheme as Wilson (Reference Wilson1982), we then simplified the problem for small capillary numbers 
$Ca^{1/3} \ll 1$. Using this asymptotic formulation, we explored the effect of some of the control parameters on the asymptotic thicknesses of the resulting thin liquid films. We mostly focused our attention on the viscosity ratio 
$M$ and the thickness of the floating layer 
${\rm \Delta} H$, as these parameters could be varied more easily in experiments. For completeness we also showed results for a limited number of values of the surface tension ratio 
$\varSigma$, although this parameter would be harder to vary experimentally using common liquids.
We could rationalise the numerical results by examining the flow in the narrow zone where the liquid/liquid and liquid/air interfaces get very close, which we dubbed virtual contact point. We showed that the shear stresses at the plate and at the liquid/liquid interface at that point alone are sufficient to predict the thicknesses of the two films, through simple scaling laws derived from the original ideas of Landau & Levich (Reference Landau and Levich1942). Regarding the influence of the different control parameters of the problem, we found that the thickness 
$\hat {h}_{1,\infty }$ of the lowermost film is barely affected by the viscosity ratio 
$M$. This is a consequence of the shear stress at the interface separating the two liquids being negligible as compared with the one at the plate/lower liquid interface. On the contrary, the thickness of the floating layer, 
${\rm \Delta} H$, has a strong impact on 
$\hat {h}_{1,\infty }$, which grows monotonically with 
${\rm \Delta} H$. In this two-liquid configuration, the thickness of the lower film is always larger than the corresponding thickness for a one-liquid Landau–Levich flow. The thickness 
$\delta \hat {h}_{\infty }$ of the uppermost film exhibits a comparatively more complex behaviour, showing a non-monotonic trend with both 
$M$ and 
${\rm \Delta} H$, with a maximum for a given pair of these parameters. More importantly, there is a finite range of values in the (
$M$, 
${\rm \Delta} H$) parameter space where 
$\delta \hat {h}_{\infty }$ takes physically realisable values, which amounts to say, where a solution with two entrained films exist.
In summary, we provided evidence that a dip-coating configuration with two entrained films is feasible using existing liquids and we developed the framework to understand and predict the corresponding entrained thicknesses. The present theory of plate coating could also be extended to fibre coating provided the curvatures of the static menisci, near the fibre and in perfect wetting conditions, are modified to account for the second principal radius of curvature contributing primarily to the capillary suction mechanism (Quéré Reference Quéré1999). In the limit of a fibre radius 
$b \ll \ell _c$, the adaptation is straightforward and consists in replacing 
$\sqrt {2}/ \ell _c$ by 
$1/b$ in the expressions for the two static menisci.
From the point of view of applications, two improvements for industrial dip-coating processes could be envisioned based on the present findings:
(i) Since adding a floating liquid layer has been shown to increase the thickness of the lower coated layer, this double-layer configuration could be used to reduce the number of passes in multi-pass dip-coating processes (Li, Haertling & Howng Reference Li, Haertling and Howng1993; Petropoulos et al. Reference Petropoulos, Foley, Hedrick and Nealey1997), without having to add undesired additives such as surfactants or nanoparticles.
(ii) Since the upper film is up to
$10^{3}$ times thinner than the lower one, this configuration could be used to deposit thin layers of a very viscous fluid (e.g. polymers), which would be impossible otherwise, at least at speeds compatible with industrial production. Naturally, this would require the removal of the lower layer, for example by taking advantage of a porous substrate (Aradian, Raphael & De Gennes Reference Aradian, Raphael and De Gennes2000). This could be of interest in the fabrication of enhanced textile (Hu Reference Hu2016) or porous membranes (Jesswein et al. Reference Jesswein, Uebele, Dieterich, Keller, Hirth and Schiestel2018).
Acknowledgements
All authors are thankful to the four referees for their very relevant comments and suggestions, which helped improving substantially the quality of the manuscript.
Funding
L.C. acknowledges funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 882429. B.S. thanks the F.R.S-FNRS for financial support. A.A.K. thanks Universidad Carlos III de Madrid and Banco Santander for the financial support of his Chair of Excellence. J.R.-R. acknowledges the support of the Spanish Ministry of Economy and Competitiveness through grants DPI2017-88201-C3-3-R and DPI2018-102829-REDT, partly funded with European funds.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Velocity field and shear stress
Here, we provide expressions for the velocity field, pressure gradients and shear stresses used to compute the data shown in figures 4 and 6. The vertical velocities, 
$\hat {u}_i$, arise from solving (2.14) and (2.15) with boundary conditions (2.19)–(2.22), while the horizontal ones, 
$\hat {v}_i$, come from solving the continuity equation. In the following expressions, a hat is used to denote rescaled variables, as introduced in § 2.3.
A.1. Velocity field in fluid 1
\begin{align} \hat{u}_1 &= u_1 = \hat{\Pi}_1\left(\frac{\hat{x}^2}{2} - \hat{h}_1 \hat{x}\right) - \hat{\Pi}_2 \hat{x} \left(\hat{h}_2-\hat{h}_1\right) + 1 ,\end{align}
\begin{align} \hat{v}_1 &= v_1\,Ca^{{-}1/3} = \hat{\Pi}'_1\left(-\frac{\hat{x}^3}{6} + \hat{h}_1\frac{\hat{x}^2}{2}\right) + \hat{\Pi}'_2\frac{\hat{x}^2}{2}\left(\hat{h}_2-\hat{h}_1\right) + \hat{\Pi}_1\frac{\hat{x}^2}{2}\hat{h}_1 \nonumber\\ &\quad +\hat{\Pi}_2\frac{\hat{x}^2}{2}\left(\hat{h}'_2-\hat{h}'_1\right) . \end{align}A.2. Velocity field in fluid 2
\begin{align} \hat{u}_2 &= u_2 ={-}\hat{\Pi}_1\frac{\hat{h}_1^2}{2} + \hat{\Pi}_2\left[\frac{1}{M}\left(\frac{\hat{x}^2}{2} - \hat{h}_2\hat{x} - \frac{\hat{h}_1^2}{2} + \hat{h}_1\hat{h}_2\right) - \hat{h}_1\left(\hat{h}_2-\hat{h}_1\right)\right]+1, \end{align}
\begin{align} \hat{v}_2 &= v_2\,Ca^{{-}1/3} = \hat{\Pi}'_1 \hat{x}\frac{\hat{h}_1^2}{2} - \frac{\hat{\Pi}'_2}{M}\left[\frac{\hat{x}^3}{6} - \frac{\hat{h}_2\hat{x}^2}{2} + \hat{x}\left(-\frac{\hat{h}_1^2}{2} + \hat{h}_1\hat{h}_2 - M\hat{h}_1\left(\hat{h}_2-\hat{h}_1\right)\right)\right] \nonumber\\ &\quad +\hat{\Pi}_1 \hat{x} \hat{h}_1 \hat{h}'_1 - \frac{\hat{\Pi}_2}{M}\left[-\hat{h}'_2\frac{\hat{x}^2}{2} + \hat{x}\left(-\hat{h}_1\hat{h}'_1 + \left(\hat{h}_1\hat{h}_2\right)' - M\left(\hat{h}_1\left(\hat{h}_2-\hat{h}_1\right)\right)'\right)\right] \nonumber\\ &\quad + C(\hat{z}) . \end{align}
The function 
$C(\hat {z})$ is the one that guarantees 
$\hat {v}_1(\hat {z}, \hat {h}_1) = \hat {v}_2(\hat {z}, \hat {h}_1)$. In these equations, the prime denotes the derivative with respect to 
$\hat z$. Moreover, the linearised rescaled pressure gradients 
$\hat {\Pi }_1$ and 
$\hat {\Pi }_2$ and their derivatives with respect to 
$\hat {z}$, 
$\hat {\Pi }'_1$ and 
$\hat {\Pi }'_2$, are given by
$$\begin{gather} \hat{\Pi}_1 = \lim_{Ca \rightarrow 0} Ca^{1/3} \varPi_1 ={-}\hat{h}_1''' - \varSigma \hat{h}_2''' \end{gather}$$
$$\begin{gather}\hat{\Pi}_2 = \lim_{Ca \rightarrow 0} Ca^{1/3} \varPi_2 ={-}\varSigma \hat{h}_2''' \end{gather}$$
$$\begin{gather}\hat{\Pi}'_1 ={-}\hat{h}_1'''' - \varSigma \hat{h}_2'''' \end{gather}$$
$$\begin{gather}\hat{\Pi}'_2 ={-}\varSigma \hat{h}_2'''' . \end{gather}$$
Note that taking the limit in expressions (A5) and (A6) amounts to simultaneously neglecting the gravity term in the definitions of the pressure gradients 
$\varPi _{i}$ ((2.16) and (2.17)) and to linearising the curvature, i.e. to setting 
$\kappa _{i} \approx h_{i}''$.
Finally, the shear stresses at the plate/liquid 1 interface, 
$\hat {\tau }_{01}$, and at the liquid 1/liquid 2 interface, 
$\hat {\tau }_{12}$, are
$$\begin{gather} \hat{\tau}_{01} = Ca^{{-}1/3} \tau_{xz}(x=0) ={-}\hat{\Pi}_1\,\hat{h}_1 - \hat{\Pi}_2\left(\hat{h}_2-\hat{h}_1\right), \end{gather}$$
$$\begin{gather}\hat{\tau}_{12} = Ca^{{-}1/3} \tau_{xz}(x=h_1) ={-}\hat{\Pi}_2\left(\hat{h}_2-\hat{h}_1\right). \end{gather}$$In figure 4(b,c), we also show the shear stress at the plate and the derivative of the pressure gradient, corresponding to the one-liquid Landau–Levich–Derjaguin problem. Their respective expressions are
\begin{equation} \hat{\tau}_{LLD} = \hat{h}_1 \hat{h}_1''' \quad \text{and}\quad \hat{\Pi}_{LLD} ={-} \hat{h}_1'''. \end{equation}Appendix B. Time-dependent formulation and numerical method
Although we ultimately seek for steady-state solutions of the problem, as described in § 2, we develop here a quasi-steady formulation (§ B.1), where the unsteady terms are kept in the mass conservation equations. This formulation allowed us to find the steady numerical solutions for 
$h_1$ and 
$h_2$ by time marching (§ B.2). The fact that we reach the steady solution by time marching ensures that the steady state is stable.
B.1. Time-dependent formulation
In the framework of our quasi-steady formulation, the dimensional thickness-averaged continuity equations (2.7) are replaced by
\begin{equation} \frac{\partial h_i}{\partial t} + \frac{\partial Q_i}{\partial z} = 0,\quad \mathrm{for}\ i = 1, 2, \end{equation}while the momentum conservation equations (2.1) and (2.2) remain unchanged. Non-dimensionalisation is performed as explained in § 2.2, making the time dimensionless with the capillary time, namely
\begin{equation} t \quad\longrightarrow\quad t \, \ell_c\mu_1/\sigma_{12}. \end{equation}
Equations (B1) are unaltered upon non-dimensionalisation and the flow rates 
$Q_i$ are still given by (2.23). Using the rescaling and simplifications developed in § 2.3, the quasi-steady equations satisfied by the film thicknesses at leading order are
\begin{equation} \frac{\partial \hat{h}_i}{\partial t} + \frac{\partial}{\partial \hat{z}} \left[\hat{h}_i + \frac{\partial^3 \hat{h}_1}{\partial \hat{z}^3}\hat{F}_{i1} + \varSigma \frac{\partial^3 \hat{h}_2}{\partial \hat{z}^3}\left(\hat{F}_{i1} + \hat{F}_{i2}\right)\right] = 0,\end{equation}
where the 
$\hat {F}_{ij}$ are the same as in the steady formulation (2.24)–(2.27).
B.2. Numerical method
The system of (B3) with the boundary conditions described in §§ 2.3 and 2.4 for 
$\hat {z} \rightarrow \pm \infty$ are discretised in space using first-order one-dimensional finite volumes in a staggered grid. The thicknesses are defined at nodes placed at the centre of the elements, while the flow rates,
\begin{equation} \hat{q}_i = \hat{h}_i + \frac{\partial^3 \hat{h}_1}{\partial \hat{z}^3}\hat{F}_{i1} + \varSigma \frac{\partial^3 \hat{h}_2}{\partial \hat{z}^3}\left(\hat{F}_{i1} + \hat{F}_{i2}\right), \end{equation}
are defined at nodes located at their boundaries. The resulting set of ordinary differential equations for the discretised thicknesses, 
$\hat {h}_i$, is then time marched with the routine odeint implemented in the scientific package SciPy of Python.
To impose the upstream boundary conditions as we approach the static meniscus, which in the scaled variables is equivalent to 
$\hat {z} \rightarrow -\infty$, we set 
$\hat {h}_1$ to a large value, say 
$\hat {h}_{1,{b}} = 100$, and then we vary 
$\hat {h}_2$ at that boundary, 
$\hat {h}_{2,{b}}$, along a range of values larger than 
$\hat {h}_{1,{b}}$. Notice that, numerically, we can impose these boundary conditions at 
$\hat z=0$ without loss of generality, thanks to the translational invariance of the problem. For every pair 
$(\hat {h}_{1,{b}}, \hat {h}_{2,{b}})$ we can then compute 
${\rm \Delta} \hat {z}_{cl} = \hat {z}_{cl,2} - \hat {z}_{cl,1}$ using (2.36) and (2.37) and, finally, relate this parameter to 
${\rm \Delta} H$ using (2.34) and (2.35).
Besides imposing 
$\hat {h}_1$ and 
$\hat {h}_2$ at the upstream boundary of the numerical domain, we also need to impose the second derivatives (taken from the static meniscus solution) there, for which purpose we use a non-centred finite-difference scheme. The same approach is used to impose the boundary conditions at the downstream boundary (corresponding physically to 
$\hat {z} \rightarrow \infty$), where we set 
$\partial \hat {h}/\partial \hat {z} = \partial ^2\hat {h}/\partial \hat {z}^2 = 0$. In our numerical method, this boundary condition is applied at 
$\hat {z} = 60$, a value sufficiently large so that the results do not depend on it.
As for the initial conditions, although the formulation introduced in this paper is able to describe transient phenomena, here, we are only interested in the long-term, steady solution. For this reason we start the time-marching procedure from initial conditions that do not represent any actual physical configuration but that satisfy the boundary conditions described above. In particular, we use a linear function of decaying exponentials, which ensures a smooth start-up of the time-marching procedure, namely,
$$\begin{gather} \hat{h}_1(\hat{z}, 0) = \sqrt{2(1-R)} \exp{({-}L\hat{z})}/L^2, \end{gather}$$
$$\begin{gather}\hat{h}_2(\hat{z}, 0) = \sqrt{\frac{2R}{\varSigma}} \exp{({-}L\hat{z})}/L^2, \end{gather}$$
where the constant 
$L$ has the value 
$L=1$ for all the simulations reported here.
Appendix C. One-liquid Landau–Levich vade mecum
In this appendix, we derive some properties of the classical one-liquid dip-coating flow (Landau & Levich Reference Landau and Levich1942) that are relevant for the discussion of its two-liquid counterpart. For a single-phase LLD flow, the film thickness 
$\hat {h}$ obeys the differential equation
\begin{equation} \hat{h}_{\infty} = \hat{h} + s \hat{h}'''\frac{\hat{h}^3}{3}, \end{equation}
where 
$s$ represents a dimensionless surface tension. For instance, if we wish to apply this equation to interface (I), 
$s = 1$, while for interfaces (II) and (III) (see figure 5) it would be 
$s = \varSigma$ and 
$s = 1 + \varSigma$, respectively. To write this equation, the thickness 
$\hat {h}$ and streamwise coordinate 
$\hat {z}$ have been made dimensionless with a capillary length times 
$Ca^{2/3}$ and 
$Ca^{1/3}$, respectively. Using this notation, 
$\hat {h}_{\infty }$ represents the dimensionless flat film thickness, far above the bath, and is equivalent also to the dimensionless flow rate transported by the film.
We start by proving that the approximate film curvature, 
$\hat {h}''$, decays monotonically with 
$\hat {z}$, i.e. moving up along the plate. In a first step, we introduce the change of variables 
$\eta = \hat {h}/\hat {h}_{\infty }$ and 
$\xi = \hat {z}/(\hat {h}_{\infty } s^{1/3})$, which yields
\begin{equation} 1 = \eta + \eta'''\frac{\eta^3}{3}, \end{equation}
where the prime denotes now the derivative with respect to 
$\xi$. In a second step, as suggested in Landau & Levich (Reference Landau and Levich1942), we take advantage of the autonomous character of the equation to reduce its order through the substitution 
$\eta ' = -F^{1/2}$. In terms of this function 
$F$, the interface curvature becomes 
$\eta '' = \frac {1}{2}\mathrm {d}F/\mathrm {d}\eta$ and (C2) turns into
\begin{equation} \frac{\mathrm{d}^2F}{\mathrm{d}\eta^2} = \frac{6\left(\eta-1\right)}{\eta^3 F^{1/2}}. \end{equation}
Since, by definition, 
$\eta > 1$, this equation reveals that 
$\mathrm {d}^2F/\mathrm {d}\eta ^2 > 0$ and therefore that 
$\mathrm {d}F/\mathrm {d}\eta = 2\eta ''$ grows monotonically with the film thickness 
$\eta$. Likewise, since the film thickness decreases monotonically with the height 
$\hat {z}$, as 
$\eta ' = -F^{1/2} < 0$, we can conclude that the interfacial curvature 
$\eta ''$ also decreases monotonically as we move up along the plate.
Regarding the flat film thickness, 
$\hat {h}_{\infty }$, it is possible to relate it to the curvature in the limit 
$\eta \gg 1$, that is, where the film connects with the static meniscus. Integrating (C3) with the initial conditions 
$F = \mathrm {d}F/\mathrm {d}\eta = 0$ at 
$\eta = 1$, as suggested in Landau & Levich (Reference Landau and Levich1942), we get 
$\mathrm {d}F/\mathrm {d}\eta (\eta \rightarrow \infty ) = 2 \eta ''(\eta \rightarrow \infty ) = 2.673$.
In terms of the original variables 
$\hat {h}$ and 
$\hat {z}$,
\begin{equation} K = \frac{\mathrm{d}^2\hat{h}}{\mathrm{d}\hat{z}^2}(\hat{z} \rightarrow -\infty) = \frac{\eta''(\eta \rightarrow \infty)}{s^{2/3} \hat{h}_{\infty}}, \end{equation}
where we denoted the curvature of the static meniscus near the wall by 
$K$. In perfect wetting conditions, 
$K = \sqrt {2}$. Finally,
\begin{equation} \hat{h}_{\infty} = \frac{1.336}{s^{2/3} K}. \end{equation}
This result shows that, the higher the curvature far away from the flat film region, the thinner the film thickness. This has been used in § 4.3 to understand the variation of lower film thickness 
$\hat {h}_{1, \infty }$ with 
${\rm \Delta} H$ and to estimate a lower bound for this quantity.
Open access
