Dip-coating flow in the presence of two immiscible liquids

Dip-coating is a common technique used to cover a solid surface with a thin liquid film, the thickness of which was successfully predicted by the theory developed by Landau&Levich and Derjaguin in the 1940's. In this work, we present an extension of their theory to the case where the dipping bath contains two immiscible liquids, one lighter than the other, resulting in the entrainment of two thin films on the substrate. We report how the thicknesses of the coated films depend on the capillary number, on the ratios of the properties of the two liquids and on the relative thickness of the upper fluid layer in the bath. We also show that the liquid/liquid and liquid/gas interfaces evolve independently from each other as if only one liquid was coated, except for a very small region where their separation falls quickly to its asymptotic value and the shear stresses at the two interfaces peak. Interestingly, we find that the final coated thicknesses are determined by the values of these maximum shear stresses.


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 1988), with applications ranging from anti-corrosion treatments and optically modified glasses, to surface functionalisation of bio-implants. The hydrodynamics of dip-coating have been known for nearly 80 years. The original description of this process dates back to the now classical work of Landau & Levich (1942), followed by the contribution of Derjaguin (1943), 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 (1982).
Nurtured by practical applications, experimental and theoretical aspects of dip-coating have kept attracting scientific attention ever since, as highlighted by the recent review by Rio & Boulogne (2017). New effects and regimes have been found to arise from the use of partially-wetting (Snoeijer et al. 2008;Tewes et al. 2019) or textured (Seiwert et al. 2011) solid substrates, non-Newtonian liquids (Hurez & Tanguy 1990;Maillard et al. 2016), surface active molecules (Park 1991;Scheid et al. 2010;Champougny et al. 2015) or jammed hydrophobic micro-particles (Dixit & Homsy 2013;Ouriemi & Homsy 2013) 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 et al. 2004), while surfactant monolayers at liquid/gas interfaces are known to exhibit analogs to two-dimensional liquid or gaseous phases (Vollhardt & Fainerman 2010).
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 1982;Liss & Duce 2005) or even more so in the event of an oil spill (Fingas 2015), hence the relevance of processes such as bubble bursting (Feng et al. 2014;Stewart et al. 2015) or bouncing (Feng et al. 2016) 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 in 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. (2017). 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 nanotubes 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. (2011) 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. 2019). 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. 2017;Kreder et al. 2018).
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 section 2, we develop the problem formulation and comment on the various dimensionless control parameters. Our results are presented and discussed in sections 3 and 4. In section 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 section 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 physicallymeaningful second film is limited to specific, finite regions in the parameter space. In section 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, section 6 is devoted to the conclusions.

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 ρ 1 and viscosity µ 1 ), covered by a layer of liquid 2 (with density ρ 2 < ρ 1 and viscosity µ 2 ), which has a thickness ∆H far from the plate. The interface between liquid 1 and 2, denoted by (I), has an interfacial tension σ 12 , while the interface between liquid 2 and the surrounding air, denoted by (II), has a surface tension σ 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 = σ 1a − (σ 2a + σ 12 ), where σ 1a denotes the surface tension of liquid 1, can be either positive or negative for the moment (see e.g. De Gennes et al. (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.
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 δh(z) = h 2 (z) − h 1 (z). The velocity fields are u 1 e z + v 1 e x and u 2 e z + v 2 e x in liquids 1 and 2, respectively.
Following the approach by Landau & Levich (1942) and Derjaguin (1943), 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: The x-independent pressures p i (i = 1, 2) are related to the atmospheric pressure p a and to the interfacial curvatures κ i through where, for i = 1, 2, 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 ∂Q i ∂z = 0 for i = 1, 2. (2.7) 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.

Non-dimensionalization
A convenient choice for the length scale to non-dimensionalise the problem is the capillary length c based on the properties of liquid 1 in the absence of the buoyancy effect created by liquid 2, defined as c = σ 12 ρ 1 g . (2.8) Note that c neither represents the scale of the static meniscus of liquid 1, nor the one of liquid 2. The "physical" dimensional capillary lengths associated to the menisci of liquids 1 and 2 are respectively c,1 = σ 12 /(ρ 1 − ρ 2 )g and c,2 = σ 2a /ρ 2 g, whose dimensionless versions indeed appear in the static configuration presented in § 2.4 (equations (2.34) and (2.35), respectively). The independent and dependent variables of the problem are made dimensionless as follows: The thickness of the layer of liquid 2 is also scaled by c , still keeping the same notation: ∆H −→ c ∆H. The capillary number is introduced as Ca = µ 1 U/σ 12 . Also, the following dimensionless parameters are defined: Based on these definitions, the dimensionless momentum equations deduced from equations (2.1) and (2.2) read where the dimensionless pressure gradients are (2.17) For 0 < z ∆H, equation (2.16) must be replaced by Notice that the pressure gradient Π 1 and its derivative are continuous at z = ∆H, as Π 2 smoothly approaches 0 as z → 0 + . In the above expressions, κ 1 and κ 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: Velocity continuity at x = h 1 : Stress continuity at x = h 1 : No slip at x = 0 : 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: Replacing the flow rates Q i given by (2.23) in the continuity equations (2.7) -which remain unchanged upon non-dimensionalization -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}, ∂h i /∂z = ∂ 2 h i /∂z 2 = 0 when z → +∞. Towards the reservoir, the static menisci connect to the flat surfaces of liquids 1 and 2. We therefore impose that κ 1 = 0 and ∂h 1 /∂z → ∞ at z = 0 for liquid 1; and κ 2 = 0 and ∂h 2 /∂z → ∞ at z = ∆H for liquid 2.

Scaling and simplification
The problem consisting of equations (2.7), along with equations (2.23) and the abovementioned 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 (1982) 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 ε as h i = εĥ i and 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 ε 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 κ i are of order unity (so dimensional curvatures are of order −1 c ) when the films approach the static menisci, we find α = 1/2 (Wilson 1982). Expanding equations (2.7) and (2.23) -(2.27) for ε 1 and retaining terms up to order ε 1/2 yields whereF ij are the functions defined in equations (2.24) -(2.27) but in terms ofĥ i instead of h i . One assumption that we make here is that the dynamic menisci region is located above z = ∆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 ε.
Analogously to the formulation developed in Wilson (1982) for the one-liquid case, equation (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 equation (2.28), viscous stresses will be of the same order as capillary ones if ε = Ca 2/3 , which is the same scaling obtained in the one-liquid case by Landau & Levich (1942). Finally, at leading order, the equations satisfied by the thicknesses are ∂ ∂ẑ The boundary conditions far from the liquid bath remain unchanged, as compared to the full formulation (developed in the last paragraph of 2.2):ĥ i =ĥ i = 0 forẑ → +∞, where the prime denotes the derivative with respect toẑ. 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 linearized and gravity was neglected.

Static configuration
For the purpose of asymptotic matching, we briefly study the static configuration depicted in figure 2. Setting Ca = 0 in equations (2.14) -(2.15) (which also implies Q i = 0) leads to Π 1 = Π 2 = 0, or equivalently: As done in the one-liquid case (see for instance Landau & Lifshitz (1987), Chapter 7, page 243), these equations can be integrated to find the exact solutions for the shape of Δz cl 0 air Figure 2. 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 ∆z cl = z cl,2 − z cl,1 .
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. Because h 1 (z cl,1 ) = h 2 (z cl,2 ) = 0 (by definition) and h 1 (z cl,1 ) = h 2 (z cl,2 ) = 0 (perfect wetting conditions on the plate), the Taylor expansions of the menisci profiles in the vicinity of the contact lines read, for z z cl,2 : (2.33) By solving (2.30) and (2.31), using the definition (2.6), we obtain the expressions for the curvatures at the contact lines, for interfaces (I) and (II) respectively.

Asymptotic matching
In the one-liquid case, thanks to the invariance of the problem in theẑ-direction, imposing the value of the curvature in the limitẑ → −∞ is enough to perform the asymptotic matching (Landau & Levich 1942). In our two-liquid case, however, we benefit from the invariance alongẑ 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 ∆H). Consequently, the boundary conditions in the limitẑ → −∞ 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 fulfill these two conditions, the film thicknessesĥ i are required to follow the parabolic approximations of the static menisci (Equations (2.32) and (2.33)), expressed in terms of the rescaled variablesĥ i = Ca −2/3 h i andẑ = Ca −1/3 z: in the limitẑ → −∞. We have introducedẑ 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 c Ca 1/3 from the contact lines of the hydrostatic solutions. For ε 1/2 = Ca 1/3 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.

Dimensionless control parameters
The problem formulated above depends on a large number of dimensionless control parameters: Ca, R, Σ, M , and ∆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 ρ 1 = 1100 kg/m 3 for the glycerol solution (Takamura et al. 2012) and ρ 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 σ 12 = 30 mN/m and σ 2a = 20 mN/m, respectively, yielding Σ = 0.667. Finally, the viscosity of the glycerol solution is µ 1 = 3.6 mPa · s (Takamura et al. 2012) and silicone oils can have a viscosity µ 2 ranging from a fraction of mPa · s to hundreds of Pa · s, 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 summarized in Table 1 and discussed in the following.
Capillary number Ca and floating layer thickness ∆H Importantly, after rescaling the problem as described in § 2.3, the capillary number Ca disappears from the formulation (see equations (2.29)), except in the matching conditions (2.36) -(2.37). In the static menisci, with which the matching is performed in the limit z → −∞, the dimensionless thickness ∆H of the floating layer remains related to the rescaled distance ∆ẑ cl between the contact lines of liquids 1 and 2 (see § 2.4) through Thus, for given fluid properties, changing the capillary number Ca has the same effect as changing the thickness ∆H. For this reason, throughout this work, we will vary only ∆H while keeping Ca = 10 −3 constant. This specific value was chosen in order to meet the condition Ca 1/3 1, ensuring that gravity effects are negligible in the dynamic meniscus region, while being relevant in many applications (Rio & Boulogne 2017). The values of ∆H explored in this work correspond to the range where solutions with two coated films were found to exist, i.e. between ∆H = 2.57 and 5.20 for our values of Ca, Σ and R.

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.

Surface tension ratio Σ
Four values of the surface tension ratio will be explored: Σ = 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.

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.

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 = σ 1a −(σ 2a +σ 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 1992). This is the case for the 40 % glycerol aqueous solution / silicone oil system described above, for which σ 1a = 70 mN/m (Takamura et al. 2012), leading to S = 20 mN/m > 0.
In the case where S < 0 however, the floating layer will be stable only above a certain critical thickness ∆H c (De Gennes et al. 2013). This critical thickness, made dimensionless with the capillary length c , can be written with our notations where S = S/σ 12 < 0 is the dimensionless spreading coefficient. For ∆H < ∆H c , the floating layer is metastable and, if perturbed, will dewet (Brochard-Wyart et al. 1993).

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 ε 1/2 1, that is to say Ca 1/3 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é 1998).

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 to viscous and capillary terms read where W e is the Weber number (Rio & Boulogne 2017), Re = ρ 1 U c /µ 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 W e ≈ 4 × 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 to gravitational ones in both liquid films requires Re Ca 2/3 1, that is to say W e Ca 1/3 , which is fulfilled for our parameters. Note that the softer condition Re Ca 2/3 ∼ 1 (i.e. W e ∼ Ca 1/3 ) is sufficient in the dynamic menisci, as gravity effects are already small there.

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 section 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).

Shape of the interfaces
Figure 3a allows us to appreciate the typical shape of interfaces (I) and (II), as well as their matching to the static menisci (parameters Σ = 0.667, R = 0.885, M = 1, Ca = 10 −3 and ∆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 equations (2.29), but plotted in their dimensionless non-rescaled form h i =ĥ 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).
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 3b and c allow us to better appreciate this convergence by looking at, on the one hand, the evolution of the horizontal distance δ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 ∆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,∞ and x = h 2,∞ , respectively. Asymptotically, the uniform and steady thicknesses bounded by these interfaces are h 1,∞ for the film of liquid 1 (referred to as lower film) and δh ∞ = h 2,∞ − h 1,∞ for the film of liquid 2 (referred to as 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, δh ∞ = 1.97 × 10 −4 while h 1,∞ = 1.33 × 10 −2 or, in terms of rescaled quantities, δĥ ∞ = δh ∞ Ca −2/3 = 1.97 × 10 −2 whileĥ 1,∞ = h 1,∞ Ca −2/3 = 1.33. Note that, with this normalization, the rescaled thickness for a one-liquid Landau-Levich-Derjaguin film would simply beĥ LLD = 0.9458, as the corresponding dimensional thickness is equal to 0.9458 c Ca 2/3 (Landau & Levich 1942). As we will see further down in section 4, this feature (δĥ ∞ ĥ 1,∞ ) 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 c,1 and c,2 , corresponding to interfaces (I) and (II) respectively, remain of the same order of magnitude: c,2 / c,1 = Σ (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,∞ and δh ∞ is connected to the very different flow patterns conveying liquid into the lower and upper films.

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 Σ = 0.667, R = 0.885, M = 1, Ca = 10 −3 and ∆H = 3.404, as functions of the vertical coordinateẑ. In figure 4a, 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 equations (A 1-A 2) for liquid 1 and (A 3-A 4) 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 entrainement 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 4b, where the dimensionless rescaled shear stress at different interfaces of interest is plotted as a function of the vertical coordinateẑ, while the latter is quantified in figure 4c, displaying the dimensionless rescaled pressure gradient as a function ofẑ.
Let us first focus on the viscous stresses, which promote liquid film entrainment. The solid black line in figure 4b represents the shear stressτ 01 , defined by equation (A 9), at the interface between the plate and liquid 1 in our two-liquid configuration. For comparison, the dashed red line shows the shear stressτ LLD at the plate/liquid interface in the one-liquid film coating problem; the corresponding expression is given by equation (A 11). 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ẑ keeps increasing. This is in sharp contrast with the trend followed by the shear stressτ 12 at interface (I), separating liquids 1 and 2 (solid blue line, defined by equation (A 10)). The shear stress at interface (I) is found to be essentially zero everywhere, except for a small peak in a narrow region (aroundẑ ≈ 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 4a 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 4c displays the dimensionless rescaled pressure gradientΠ 1 (blue) and Π 2 (orange) in liquid phases 1 and 2, defined by equations (A 5) and (A 6), respectively. Again, the two liquid phases exhibit very different behaviors: the lower liquid is subjected to a mild pressure gradient, spread over distances of several units inẑ. On the contrary, there is virtually no pressure gradient in the upper liquid, except for a high peak around z ≈ 24, which also corresponds to the region of highest shear stress (in absolute value) at the liquid 1/liquid 2 interface.

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,τ 12 , and in the pressure gradient in the upper film,Π 2 . Figures 4b  and 4c show that this region is very narrow in the vertical direction, as compared to the overall extension of the dynamic menisci. Moreover, figure 3b reveals that, concomitantly, the distance between interfaces (I) and (II) decays quickly down to the asymptotic one, δh ∞ . For these reasons, we find meaningful to model this region as a point, called 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 * . In the following, all quantities evaluated at this point will be marked by an * .
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 δh h 1 , and asymptotically δh ∞ h 1,∞ , downstream 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 + Σ. In this framework, the different interfaces pictured in figure 5 can be described as follows.
• 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.
• 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 4c.
• 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.
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 rationalized, and even quantitatively predicted to some extent, looking at the flow variables at the virtual contact point.

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 rationalized 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. (2015) -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ĥ 1,∞ for liquid 1 and δĥ ∞ for liquid 2.

Lower film thicknessĥ 1,∞
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: − τ 01, * ∼ µ 1 U h 1, * . (3.1) 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, * =ĥ 1, * c Ca 2/3 and τ 01, * =τ 01, * (σ 12 / c ) Ca 1/3 , respectively, yielding −τ 01, * ĥ1, * ∼ 1. Making the approximation that the film thickness has already reached its asymptotic value at the virtual contact point, namely that h 1,∞ ≈ĥ 1, * , we arrive atĥ 1,∞ ∼ −τ −1 01, * . (3.2) In figure 6a, we plot the lower film thicknessĥ 1,∞ as a function of the inverse of the shear stress at the virtual contact point, −τ −1 01, * , for a variety of Σ, M and ∆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ĥ 1,∞ to be larger than that of a one-liquid Landau-Levich filmĥ LLD = 0.9458.

Upper film thickness δĥ ∞
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 the streamwise extension of the virtual contact point region (typically the width of the peaks inτ 12 and Π 2 in figure 4), this force balance reads in dimensional terms: Additionally, matching the curvatures of the dynamic and static menisci of liquid 2 imposes (still using dimensional quantities) where we recall that c,2 = σ 2a /ρ 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 (1942) in the one-liquid dip-coating problem. Combining equations (3.3) and (3.4), we arrive at the dimensional scaling expression where we approximated δh ∞ ≈ δh * . Using the definitions δh ∞ = δĥ ∞ c Ca 2/3 and τ 12, * =τ 12, * (σ 12 / c ) Ca 1/3 , we finally express equation (3.5) in terms of the rescaled dimensionless variables: In figure 6b, we plot the upper film thickness δĥ ∞ as a function of the right-hand term of equation (3.6) for a variety of Σ, M and ∆H. When plotted in that way, the numerical data is found to collapse on a reasonably linear master curve. The linear fit (solid black line) yields a prefactor equal to 1.29, showing that equation (3.6) can be used to estimate the upper film thickness δĥ ∞ . 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 (τ 01, * for liquid 1,τ 12, * 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 (1942).

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,ĥ 1,∞ for the lower film and δĥ ∞ for the upper film. The main control parameters varied in our study are the dimensionless floating layer thickness ∆H, the viscosity ratio M and the surface tension ratio Σ (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 ∆H ( § 4.2 and § 4.3).

Thickness maps in the M − ∆H parameter space
In figure 7, the asymptotic film thicknesses are presented as color maps in the M −∆H parameter space for four different values of the surface tension ratio Σ. The left column (panels a, c, e, g) shows the lower film thicknessĥ 1,∞ while the right column (panels b, d, f, h) displays the upper film thickness δĥ ∞ . For an easier quantitative reading, cuts along  Together, figures 7 and 8 reveal dramatic qualitative and quantitative differences between the lower and upper coated films. As already observed and rationalized in section 3, the lower film exhibits final thicknessesĥ 1,∞ of order unity (as expected from the Landau-Levich scaling, see § 3.4), while the upper film reaches steady-state thicknesses δĥ ∞ that are about 10 −3 to 10 −2 times smaller. Not only the lower and upper film thicknesses are very disparate, but they also depend very differently on the control parameters. Figures 7a, c, e, 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 ∆H. On the contrary, figures 7b, d, f, h show that the thickness δĥ ∞ of the upper film depends non-monotonically on both M and ∆H, exhibiting a maximum whose exact position in the (M, ∆H) parameter space depends on the surface tension ratio Σ.
Around this maximum, the upper film thickness decreases in all directions, eventually going down to values reaching our numerical accuracy (δĥ ∞ ≈ 5 × 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 figures 7b, d, f, h). The zones of the (M, ∆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ĥ 1,∞ can still be obtained outside these islands (see shaded areas in figures 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.

Minimum and maximum ∆H for the existence of two films
The very different scales on the M and ∆H axes in figure 7 highlight that the existence islands extend only along a finite range of values of ∆H, getting narrower and narrower as Σ increases, while they span several orders of magnitude in M values. In this paragraph, we aim at providing some physical arguments to rationalize this observation.
The existence islands exhibit a relatively well-defined lower boundary in ∆H, which becomes independent of the viscosity ratio M for M 1. This limit can be understood from hydrostatics 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, that the apparent contact line of the upper interface lays above that of the lower one. For this to occur, the distance between the apparent contact lines of the lower and the upper menisci, ∆z cl , must be positive. Combining equations (2.34) and (2.35) this condition translates into The values of ∆H min predicted by equation (4.1) for surface tension ratios Σ = 0.33, 0.67, 1.00 and 1.27 are presented in table 2 and displayed in figures 7b, d, f, h as dashed white lines. Comparison with the numerical data reveals a good agreement with the lower limit of the existence islands for M 1. Regarding the maximum ∆H for the existence of a double-film configuration, ∆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 ∆ẑ ∼ 10 in the streamwise direction, as can be seen in figure 4b. This means that if the virtual contact point -where the two interfaces come close to each other -is outside this limited region inẑ, no shear can be transferred to the liquid / liquid interface, and thus the upper film cannot be dragged. In terms of ∆H, equation (2.38) allows us to translate the extension ∆ẑ ∼ 10 into ∆H max − ∆H min ∼ 10 Ca 1/3 . (4.2) The condition Ca 1/3 1 (see § 2.3 and § 2.6) explains why we expect the existence islands to span only along a reduced extension in ∆H. In our particular case where Ca = 10 −3 , equation (4.2) predicts island extensions of the order of unity in ∆H, which is compatible with the data displayed in figure 7.

Effect of ∆H on the lower film thicknessĥ 1,∞
In this paragraph we use the equivalent description presented in § 3.3 to explain the effect of the floating layer thickness ∆H on the asymptotic thicknessĥ 1,∞ ≈ĥ III,∞ 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.
• Interface (I) -Since interface (I) behaves as in the one-liquid LLD theory, its curvatureĥ I decays monotonically withẑ (see appendix C). The maximum curvature, 2 (1 − R), is found in the limitẑ → −∞, corresponding to the lower static meniscus.  figure 7 and the corresponding values predicted using simplified approaches (section 5), for different values of the surface tension ratio Σ. The quantities extracted from figure 7 are evaluated in the limit M 1 and within the existence islands (red dashed contours), where a double coating solution exists. These quantities are (a) the lower limit in ∆H of the existence islands and (b) the minimum values of the asymptotic lower film thicknessĥ1,∞.
• 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: h II = 2R/Σ.
• Interface (III) -For the same reason as interface (I), the curvature of interface (III) decreases monotonically withẑ, reaching its maximum valueĥ III, * at its lowest point, i.e. the virtual contact point.

(4.3)
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 + Σ)ĥ 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 Σĥ II, * andĥ I, * . Requiring pressure continuity at the virtual contact point leads tô Combining equations (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:ĥ (4.5) Monotonous behaviour ofĥ 1,∞ with ∆H This simplified view allows us to explain whyĥ 1,∞ grows monotonically with the thickness of the floating layer, ∆H. As this parameter grows, the virtual contact point, where interfaces (I) and (II) meet, displaces up downstream or, in other words,ẑ * increases. Sinceĥ I is a decreasing function ofẑ, the largerẑ * , the lower the corresponding curvatureĥ I, * of interface (I). We therefore deduce thatĥ I, * decays with ∆H. Finally, equation (4.5) allows us to conclude thatĥ 1,∞ must be a monotonously increasing function of ∆H, as observed in our numerical results in figures 7 and 8.
Lower bound forĥ 1,∞ 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 (ẑ → −∞, i.e. closer to the liquid bath), it will eventually reach the region where the curvature of interface (I) has its asymptotic (maximum) value:ĥ I, * = 2 (1 − R). Substituting this value in equation (4.5) yields a lower bound for the lower film thickness: (4.6) We can compare the predictions of equation (4.6) to the values ofĥ 1,∞ 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 1, in which the simplified geometrical model is valid. As shown in table 2, the minimum values min (ĥ 1,∞ ) obtained from the simulations are not only compatible with the lower bounds provided by equation (4.6) but also follow a similar trend with Σ.

Discussion
In this last section, we discuss some limitations of our hydrodynamic model in relation to practically relevant effects, such as partial wetting conditions betweeen the two liquids and long range intermolecular forces across the upper film, when it is sufficiently thin.

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 S = S/σ 12 = Σ − (Σ + 1), where we have introduced Σ = σ 1a /σ 12 . For S > 0, the floating layer is stable regardless of its thickness while, for S < 0, only floating layers of dimensionless thickness ∆H > ∆H c given by equation (2.39) are stable (De Gennes et al. 2013). In § 4.2, we estimated the minimum floating layer thickness ∆H min needed for a double liquid film to be entrained (equation (4.1)). In partial wetting conditions, our hydrodynamic description of the configuration of figure 1 is therefore warranted as long as ∆H min > ∆H c . This translates into the following condition on the dimensionless spreading parameter (5.1) or, equivalently, on the dimensionless liquid 1 / air surface tension Σ

Thinness of the upper film and disjoining pressure effects
As can be seen in figure 7, the dimensionless rescaled thickness of the upper film, δĥ ∞ , is of the order of 10 −3 to 10 −2 in the range of parameters explored. The corresponding dimensional thickness, δh ∞ = c Ca 2/3 × δĥ ∞ , is therefore expected to be in the range 20 − 200 nm for the typical values Ca = 10 −3 and c ≈ 2 mm. Albeit small, these thicknesses are withing reach of techniques such as reflection interference contrast microscopy (RICM), as examplified by the recent measurements of Kreder et al. (2018) 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 1978), 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 2011). Following the approach of Champougny et al. (2017), the effect of van der Waals interactions can be included in our formulation by adding a disjoining pressure gradientΠ d to the capillary pressure gradientsΠ 1 andΠ 2 defined in Appendix A. Considering only non-retarded van der Waals forces (Léger & Joanny 1992) and using the rescaling exposed in § 2.2 and § 2.3, the disjoining pressure gradient readŝ c is a dimensionless Hamaker constant. Using the values of Ca, σ 12 and c reported in § 2.5, as well as the typical dimensional Hamaker constant |A| ≈ 10 −20 J for water-oil-air systems (Israelachvili 2011), we estimate A ≈ 1.6 × 10 −14 . The upper film thickness δĥ crit ∞ at which the disjoining pressure gradient becomes of the order of the capillary one can be estimated by requiring: In the rescaled variables used here,ĥ 2,∞ ∼ĥ 1,∞ is of order unity, and so is the characteristic distance ∆ẑ over which the thicknesses vary. Thus, we haveĥ 2 ∼ 1 and h 2 −ĥ 1 ∼ĥ 2 −ĥ 1 = δĥ ∞ , yielding With the values of A and Ca given above, we arrive to thresholds δĥ crit ∞ ∼ 3.6 × 10 −3 , 2.9 × 10 −3 , 2.5 × 10 −3 and 2.3 × 10 −3 for Σ = 0.333, 0.667, 1.00 and 1.27, respectively. For each Σ, this means that the parts of the existence islands where δĥ ∞ δĥ crit ∞ (essentially the edges, see figure 7) may be affected by disjoining pressure. Conversly, 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 2005;Bandyopadhyay & Sharma 2006) lies beyond the scope of the present work.

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 (1982), we then simplified the problem for small capillary numbers Ca 1/3 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 ∆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 Σ, although this parameter would be harder to vary experimentally using common liquids.
We could rationalize 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 (1942). Regarding the influence of the different control parameters of the problem, we found that the thicknessĥ 1,∞ 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 to the one at the plate / lower liquid interface. On the contrary, the thickness of the floating layer, ∆H, has a strong impact onĥ 1,∞ , which grows monotonically with ∆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 δĥ ∞ of the uppermost film exhibits a comparatively more complex behaviour, showing a non-monotonous trend with both M and ∆H, with a maximum for a given pair of these parameters. More importantly, there is a finite range of values in the (M , ∆H) parameter space where δĥ ∞ takes physically-realizable values, which amounts to say, where a solution with two entrained films exists.
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 fiber coating provided the curvatures of the static menisci, near the fiber 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é 1999). In the limit of a fiber radius b c , the adaptation is straightforward and consists in replacing √ 2/ 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 et al. 1993;Petropoulos et al. 1997), without having to add undesired additives such as surfactants or nanoparticles.
(ii) Since the upper film is up to 10 3 times thinner than 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 et al. 2000). This could be of interest in the fabrication of enhanced textile (Hu 2016)  (2.16) -(2.17)) and to linearize the curvature, i.e. to set κ i ≈ h i .
Finally, the shear stresses at the plate / liquid 1 interface,τ 01 , and at the liquid 1 / liquid 2 interface,τ 12 , are: In figure 4b and 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 arê

Appendix B. Time-dependent formulation and numerical method
Although we ultimately seek for steady-state solutions of the problem, as described in section 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 Equations (B 1) are unaltered upon non-dimensionalization and the flow rates Q i are still given by (2.23). Using the rescaling and simplifications developed in § 2.3, the quasisteady equations satisfied by the film thicknesses at leading order are where theF ij are the same as in the steady formulation (equations (2.24) -(2.27)).

B.2. Numerical method
The system of equations (B 3) with the boundary conditions described in § 2.3 and § 2.4 forẑ → ±∞ are discretized in space using first-order one-dimensional finite volumes in a staggered grid. The thicknesses are defined at nodes placed at the center of the elements, while the flow rates, are defined at nodes located at their boundaries. The resulting set of ordinary differential equations for the discretized thicknesses,ĥ 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ẑ → −∞, we setĥ 1 to a large value, saŷ h 1,b = 100, and then we varyĥ 2 at that boundary,ĥ 2,b , along a range of values larger thanĥ 1,b . Notice that, numerically, we can impose these boundary conditions atẑ = 0 without loss of generality, thanks to the translational invariance of the problem. For every pair (ĥ 1,b ,ĥ 2,b ) we can then compute ∆ẑ cl =ẑ cl,2 −ẑ cl,1 using equations (2.36) -(2.37) and, finally, relate this parameter to ∆H using equations (2.34) -(2.35).
Besides imposingĥ 1 andĥ 2 at the upstream boundary of the numerical domain, we also need to impose the second derivatives (taken from the static menisci solution) there, for which purpose we use a non-centered finite-difference scheme. The same approach is used to impose the boundary conditions at the downstream boundary (corresponding physically toẑ → ∞), where we set ∂ĥ/∂ẑ = ∂ 2ĥ /∂ẑ 2 = 0. In our numerical method, this boundary condition is applied atẑ = 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, what ensures a smooth start-up of the time-marching procedure, namely, 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 1942) that are relevant for the discussion of its two-liquid counterpart. For a single-phase LLD flow, the film thicknessĥ obeys the differential 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 = Σ and s = 1 + Σ, respectively. To write this equation, the thicknessĥ and streamwise coordinateẑ have been made dimensionless with a capillary length times Ca 2/3 and Ca 1/3 , respectively. Using this notation,ĥ ∞ 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,ĥ , decays monotonically withẑ, i.e. moving up along the plate. In a first step, we introduce the change of variables η =ĥ/ĥ ∞ and ξ =ẑ/(ĥ ∞ s 1/3 ), which yields where the prime denotes now the derivative with respect to ξ. In a second step, as suggested in Landau & Levich (1942), we take advantage of the autonomous character of the equation to reduce its order through the substitution η = −F 1/2 . In terms of this function F , the interface curvature becomes η = 1 2 dF/dη and equation (C 2) turns into d 2 F dη 2 = 6 (η − 1) η 3 F 1/2 . (C 3) Since, by definition, η > 1, this equation reveals that d 2 F/dη 2 > 0 and therefore that dF/dη = 2η grows monotonically with the film thickness η. Likewise, since the film thickness decreases monotonically with the heightẑ, as η = −F 1/2 < 0, we can conclude that the interfacial curvature η also decreases monotonically as we move up along the plate.
In terms of the original variablesĥ andẑ, where we denoted the curvature of the static meniscus near the wall by K. In perfect wetting conditions, K = √ 2. Finally,ĥ ∞ = 1.336 s 2/3 K .

(C 5)
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ĥ 1,∞ with ∆H and to estimate a lower bound for this quantity.