Dynamics of Thin Film Under a Volatile Solvent Source Driven by a Constant Pressure Gradient Flow

: The evolution of a thin liquid ﬁlm subject to a volatile solvent source and an air-blow effect which modiﬁes locally the surface tension and leads to Marangoni-induced ﬂow is shown to be governed by a degenerate fourth order nonlinear parabolic h -evolution equation of the type given by ∂ t h = − div x (cid:0) M 1 ( h ) ∂ 3 x h + M 2 ( h ) ∂ x h + M 3 ( h ) (cid:1) , where the mobility terms M 1 ( h ) and M 2 ( h ) result from the presence of the source and M 3 ( h ) results from the air-blow effect. Various authors assume M 2 ( h ) ≈ 0 and exclude the air-blow effect into M 3 ( h ) . In this paper, the authors show that such assumption is not necessarily correct, and the inclusion of such effect does disturb the dynamics of the thin ﬁlm. These emphasize the importance of the full deﬁnition (cid:126) t · grad ( γ ) = grad x ( γ ) + ∂ x h grad y ( γ ) of the surface tension gradient at the free surface in contrast to the truncated expression (cid:126) t · grad ( γ ) ≈ grad x ( γ ) employed by those authors and the effect of the air-blow ﬂowing over the surface.


Introduction
In this work, we study the influence of a volatile solvent on the evolution of a thin liquid film on a solid surface. As the volatile solvent diffuses through the atmosphere, a non-uniform solvent concentration arises which induces surface tension gradient and Marangoni-driven flow [1][2][3]. This effect is well-known and was described in the pioneering works of Marangoni [4] and Thomson [5].
This physical phenomena is exploited in the industrial process, and is known as Marangoni drying whereby a jet of ethanol vapor (or other solvent) is blown onto a wet surface. As the ethanol interacts with the liquid free surface, a Marangoni flow is generated which helps drying the surface. This process has motivated a number of studies [6][7][8][9]. On experimental side, Marangoni drying was first studied by Leenaars et al. [6]. They coined the name Marangoni drying to this process. Leenaars et al. [6] have shown that this drying technique is more effective than spin drying which is usually carried out by centrifugation: Marangoni drying yields to extremely clean surface whereas spin drying can contaminate the surface. On the theoretical side, Matar and Craster [7] have given a mathematical description based on this drying process. The mathematical framework of Matar and Craster [7] permits the quantitative evaluation of the thinning process due to Marangoni drying. The model also makes it possible through parametric study to examine the response of the liquid film to optimize the thinning process. O'Brien [8] has considered the dynamics of a liquid film subjected to the action of alcohol-vapor-induced Marangoni effects using asymptotic methods. Taking into account only terms of leading orders into the generic equations, O'Brien reduces the mathematical description to a single non-linear partial differential equation whose characteristics are similar to the so-called Korteweg-de Vries and Burgers equations.
Preluding these studies, Carles and Cazabat had already demonstrated that the surrounding atmosphere plays an important role in the dynamics of a droplet spreading [10,11]. More specifically, these authors observe that the usual laws of spreading no longer hold. For some droplets, the spreading was strongly accelerated and accompanied by instabilities at the contact line, whereas for others, the process was somewhat reversed: retraction of droplets after a fast initial spreading. In particular, when a non-volatile droplet was investigated under an inhomogeneous surrounding of volatile solvent, they observe both the Marangoni effect and instabilities at its edge (see Reference [10]). In the second paper of Carles and Cazabat (see Reference [11]), the spreading dynamics of droplets under an inhomogeneous vapor phase of its volatile solvent were further investigated. The authors have investigated the spreading of an oil droplet under an atmosphere contaminated with a volatile solvent. Still, the observed behaviors were completely different from the normal ones. For, in some situations the surface tension gradient counter-accelerates the spreading mechanism while in others it accelerates instead, both along its directions. After a short time, when the interfacial effect was no longer operative, it was observed that the droplet either resumes spreading or retracts, again depending on the surface tension of the solvent (see Reference [11]).
Much research has been devoted to the understanding of Marangoni-driven thin liquid films. An exhaustive review of such studies is beyond the scope of this introduction but the interested reader is referred to the reviews by Oron et al. [12] and Craster and Matar [13]. Most of these study have considered a leading order expression for the surface tension gradient at the free surface. One aim of this paper is to explore the effect of these terms which have typically been disregarded by others.
Another important aspect in the study of Marangoni drying is the two way-coupling between the liquid film and the surrounding atmosphere. Few papers appear to have investigated this effect. One notable exception is the work of Sultan et al. [14] who studied the stability of a two-dimensional bi-layered liquid-vapor system over a solid substrate. In this paper, the authors generalize the one-sided study of Burelbach et al. [15] on the evaporation of thin film by including the diffusion of the vapor phase region. Their results are described in terms of both interfacial and mass transport phenomena. In contrast to earlier work, we consider here diffusion driven by an air-blow effect in the surrounding atmosphere and a two-way coupling whereby the air-blow effect in the vapor phase affects the liquid film. This paper is organized as follows: In Section 2.1, the mathematical descriptions is exposed in dimensional form; the spatial description of a volatile solvent and the character of the interfacial force it induces are argued. Thereafter, Sections 2.2 and 2.3 are devoted to perturbation theory so as to render the mathematical descriptions of earlier sections dimensionless and at leading order. Section 2.4 concerns the numerical side, yet some leading features are exposed in Appendix A. Finally, the results are discussed through a series of case studies under Section 3 and concluded on this basis under Section 4.

Modeling the Inhomogeneous Vapor Phase Region
Following the works of Carles et al. (see References [10,11]), we describe the presence of a volatile solvent in the atmosphere, region Ω v ⊂ R 2 , by a volatile source function S ( x) of the form In Equation (1), x 0 = (x s , y s ) ∈ Ω v is the point at which the source is located. Further, the scalar quantities M s , µ s , λ (v) > 0 define the rate of increase in moles of solvent per unit length, the strength of the source and the diffusion constant of the volatile solvent vapor. The time scale T > 0 serves to establish a diffusion length l T = 2 √ λ (v) T. In the following, we will consider the transport of volatile solvent from a point source given by Equation (1).

On the Surface Tension and Its Derivative
For some real t ∞ < ∞, let the set I t = (0, t ∞ ) denote the period of time over which the study is carried out. Suppose the scalar field function c (l) (t, x) be the spatio-temporal distribution of the chemical concentration in the liquid domain Ω l , where (t, x) ∈ I t × Ω l . Then, according to the experimental work of Carles et al. [10], the closure relationship relating the scalar field function c (l) (t, x) to the surface tension γ (t, x) at every (t, x) ∈ I t × Γ 4 is linear and hence, where κ 1 > 0 and κ 2 ≤ 0 are property-dependant constants. To this end, the directional derivative of the scalar field function γ (t, x) along the unit vector e x yields to To the best of our knowledge, the effect of the term ∂ x h∂ y c (l) was disregarded in the past. In the sequel, we will show that it has a non-negligible effect. Before moving on to the next section, we deduce the constants κ 1 , κ 2 of the above expression. Let γ 0,l > 0 be the surface tensions of the liquid in its uncontaminated state. Then, it is evident that one of these constants describes the surface tension of the liquid film when the liquid film is chemically stable, and the other equates the surface tension gradient of the liquid whenever that stable condition is being disturbed chemically. Thus, these equalities follow: respectively, where D η (·) denotes the ordinary derivative operator with respect to η; that η = c (l) in the above description is understood but when η = t, D η (·) defines itself as the Lagrangian derivative operator. Consequently, it results by substitution that,

Constitutive Equations
Consider a thin liquid film initially sitting on a solid substrate Γ 2 , having a thickness h 0 > 0 and a length L > 0 (with ord h 0 /L ∼ 10 −3 , see References [7,8]). For mathematical treatment this liquid film is designated by Ω l ⊂ R 2 and bounded by ∂Ω l = 4 i=1 Γ i . Henceforth, subscripts v and l indicate the vapor and liquid phases, respectively. The vapor phase surrounding is defined by Ω v ⊂ R 2 and bounded by ∂Ω v = 7 i=4 Γ i . The volatile solvent source is placed at x 0 = (x s , y s ) ∈ Ω v ; Figure 1 supports such geometrical descriptions. The presence of the source S ( x) gives rise to advection-convection phenomena. Therefore, if c (k) (t, x) designates the chemical concentration in Ω k , where (t, x) ∈ I t × Ω k and k ∈ {l, v}, then the following advection-diffusion equations result in: where the bi-subscripts δ kv ∈ {0, 1} is the Kronecker delta (with k ∈ {l, v}). The operator D t (·) = ∂ t (·) + u · grad (·) designates the Lagrangian derivative operator. Evidently, λ (l) > 0 is the diffusion constant of the volatile solvent into the liquid film, and u = u (t, x) for all (t, x) ∈ I t × Ω k (with k ∈ {l, v}); u (t, x) is called the velocity vector field distribution within either Ω l or Ω v . For every k ∈ {l, v}, the flow field u (t, x) within Ω k is described by the incompressible Navier-Stokes equations of constant fluids properties in a gravitational field, which in vector invariant forms read: where the tensor field T (p, u), which stands for the total stress tensor, is given by In Equation (4), tensors I and E ( u), respectively, are the identity and strain-rate tensors, and fields p (t, x), u (t, x) = u x , u y T and f (t, x) = f x , f y T , respectively, denote the pressure, velocity, and body force (per unit mass). The property constants ρ > 0 and µ > 0, respectively, define the density and dynamic viscosity of the k th fluid (with k = l for liquid and k = v for vapor phases). The pressure difference sustained across the interface Γ 4 is modeled with the Young-Laplace equation. Thus, (We note here that grad p (l) = − grad γ ∇ · n if the atmosphere has constant pressure; in the sequel, an alternative form for grad p (v) is established such that the latter operates as an air-blow effect.) In addition, it is understood that by the vector field u (t, x) defined in Ω l , we implicitly mean u (l) (t, x), and so forth for other property constants and variables. The vector fields n (t, x), t (t, x) are the outer unit normal and tangent vectors of their boundary. In particular, on Γ 4 (see Figure 1), we obtain for all (t, x) ∈ I t × Γ 4 .

Initial and Boundary Conditions
At time t = 0, the domains Ω l and Ω v are initially motionless and free from or low in volatile substances. Thus, these initial conditions follow: In these conditions, c ≥ 0 are arbitrary constants which will be used in the nondimensionalization process so as to obtain a zero concentration in dimensionless forms. On boundaries Γ 5 , Γ 7 , we impose a constant pressure gradient in the atmosphere. If K (v) > 0 stands for such constant and · denotes the Euclidean norm in R 2 , then On the boundaries Γ 1 , Γ 3 of the domain Ω l , we prescribe a so-called normal stress-free boundary conditions. Thus, for every i ∈ {1, 3}, where the scalar field p (l) 0 > 0 is an arbitrary constant and n (t, x) is the outer unit normal vector of the boundary in question. For the chemical concentrations, we impose a flux continuity boundary condition on the interface Γ 4 , Further, we assume the mass transfer process to vanish on Γ 6 and to result in a zero-flux boundary conditions on boundaries Γ 1 to Γ 3 and Γ 5 to Γ 7 . Thus, On the surface Γ 2 , we consider a Navier slip-with-friction conditions, where β > 0 is the Navier-slip coefficient. For, the movement of the liquid film entails a singularity of stress at the contact line if the usual no-slip boundary condition, i.e., t · E ( u) n = 0 on Γ 2 , is imposed between the liquid film and the solid substrate (see Reference [16] and references therein). On the interface Γ 4 , we require that the substantial derivative of y − h (x, t) vanishes. Therefore, Finally, to ensure continuity upon tangential shear stress on Γ 4 , we consider for all (t, x) ∈ I t × Γ 4 .
A remark deserves attention at this point. It is very true that Figure 1 does not indicate the presence of a contact line and, therefore, the discussion of the Navier slip-with-friction condition seems odd. The authors agreed that the Navier-slip-with-friction conditions on Γ 2 should instead be the classical no penetration boundary conditions; that is, the normal and tangential velocity vector fields vanish for all (t, x) ∈ I t × Γ 2 , n · u = 0 and t · u = 0. However, we wish to derive a degenerate fourth order nonlinear parabolic h-evolution equation that will apply not only in the case considered here (Navier-slip-with-friction coefficient β > 0 set to zero) but also to the case in which a contact line (or trijunction) exists and a liquid on a solid substrate spreads and displaces the surrounding fluid (say, vapor phase region Ω v ).

Asymptotic Approximations
Set = h 0 /L 1 and, by hypothesis, assume ord (h 0 /L) ≈ 10 −3 ; obviously, denotes the asymptotic parameter. Geometrical properties are scaled based on the initial film thickness h 0 = L and physical properties are on the contrary scaled based on interfacial and viscous forces. On By q → αq is meant that the quantity q is dimensional if it precedes → and dimensionless if it follows →; the quantity α is obviously dimensional. The parameters Fr l , Pe l , Re l > 0 are, respectively, the Froude, Péclet, and Reynolds numbers in Ω l . By substitution of the above scalings into the dimensional equations, one obtains a system of dimensionless variables equations.
/L 2 because the spreading of the volatile solvent is by diffusion process. On the whole, the following scalings hold reasonable in the domain Ω v : Substituting the above scalings into Equation (3), there results in In a like manner, the dimensionless forms of Equation (4) respectively. The surface tension of the above dimensionless scaling set is termed the spreading pressure, [7], where γ 0,k > 0 is the surface tension of the k th fluid (with k = l, v), as already pointed out. The quantities c ∞ are prescribed values characterizing the initial and final values of the concentrations of the k th fluids (with k = l, v). The control parameters Fr v , Pe v , Re v > 0, respectively, are obviously the Froude, Péclet, and Reynolds numbers in Ω v .
It now remains to nondimensionalize the set of initial and boundary conditions. Beginning with the initial conditions given in Section 2.1.4 for the field functions u (k) (t, x), c (k) (t, x) (with k = l, v), respectively, we obtain, for every k ∈ {l, v}, where the property constants c Likewise, nondimensionalizing the normal stress-free boundary condition given in Section 2.1.4 yields to −p + ∂ x u x = −P (l) . Imposing on boundaries Γ 1 , Γ 3 , respectively, a zero viscous stress along with a Dirichlet condition on p ( x, t) gives To nondimensionalize the flux-continuity boundary condition, one should be careful. In fact, away 0 , the dimensionless form of the flux-continuity boundary condition writes In like manner, after nondimensionalizing the other boundary conditions given in Section 2.1.4, one deduces For the dimensionless form of the so-called Navier slip-with-friction conditions, we proceed in two steps. In the first place, we nondimensionalize the conditions given there. Then, we replace the Navier-slip coefficient β by (2 βL) −1 . This gives Likewise, after nondimensionalizing the kinematic boundary conditions given in Section 2.1.4, we obtain Note that the dimensionless unit normal vector and tangential vectors n (t, x) and t (t, x), respectively, defined in Section 2.1.3 are approximated as follows, ord ( n) ∼ (− ∂ x h, 1) T , ord t ∼ (1, ∂ x h) T , respectively. We now wish to combine Equation (2) with Equation (5), given in Section 2.1.3. To do so, define two dimensionless numbers, namely, the Capillary (Ca) and Marangoni (Ma) numbers as respectively. The construction of these dimensionless numbers Ca, Ma > 0 is based on the following line of thought. Since ord u

in order to establish the
Capillary number Ca, it suffices to substitute (γ 0,l − γ 0,v ) /µ and µλ (l) /L for u (l) 0 and γ, respectively, into Ca = µu (l) 0 /γ. On the one hand, the Marangoni number Ma being regarded as proportional to (thermal-)surface tension forces divided by viscous forces, i.e., Ma = −∂ T γ (T ∞ − T 0 ) L/µα, where T and α, respectively, designate the temperature and thermal diffusivity, so to adapt this dimensionless number with the problem under consideration, it suffices to replace the field T by c (l) , and the coefficient α by λ (l) . Thence, the dimensionless numbers so far defined. Using those dimensionless numbers into Equation (2), the resulting expression in dimensional form reads, Combining the above expression with Equation (5), the resulting dimensional expression writes Let us now investigate into p (v) (t, x). In the immediate vicinity of every point which, in other words, suggests to scale δp as such, δp → −1 (γ 0,l − γ 0,v ) δp/L; obviously, the differential quantity δp is dimensional if it precedes → and dimensionless otherwise. Hence, taking into account the characteristic scales given in Section 2.2.1, setting γ 0,l → γ 0,l (γ 0,l − γ 0,v ) (γ is dimensional if it precedes → and dimensionless otherwise) and then rescaling the dimensionless The dimensionless form of the (chemical)-surface tension is now given by the expression γ c (l) = γ 0,l + (Ma / Ca) c (l) .

Leading Order Model
To begin, some simple analysis is now in order.
(These numerical values are extracted from References [7,8].) On the other hand, we may assert that ord Pe −1 l ∼ 2 , since the spreading pressure π s is a predominant factor. At this point, we are now well-prepared for the search of a simplified model.
In the limit that the quantities 2 and Pe v , respectively, tend to zero and to unity, the leading order forms of the dimensionless equations for the fields Some analysis is now in order to support the assumption made upon Pe l . Suppose that the flow induced within the region Ω l is an −|k| -order of the Péclet number, i.e., ord (Pe l ) ∼ −|k| (with k = 0). Then, if we express Fr l , Pe l and Re l in terms of the quantity (γ 0,l − γ 0,v ) /µ it can be shown that And, consequently, where ν ∼ |k| g (with k = 0). The above approximation enables us to obtain some information as to the nature of the order of magnitube of one dimensionless number with respect to the other. Moreover, for the case of water, we found k ≈ 2, based on the data of O'Brien [8]; hence, the assertion made in so far. With this view in mind, taking k ≈ 2 there follows that Fr l , which simply asserts a small Bond number (Bo l ), since Bo l = Re l / Fr l 1. Consequently, the system of Equation (7) reduces to On the other hand, paying little attention to the effects arising from body and viscous forces in Ω l , i.e., Fr −1 v , Re −1 v 1, the system of Equation (8) results in

Alternative Form for the Pressure Field in the Vapor Phase
We postulate that the pressure field in Ω v varies according to where j ≥ 1. Consequently, the leading order of the flow equations in Ω v in vector form writes where j ≥ 1. Therefore, D t u (v) + K (v) e x = 0 at leading order; that is, we force the acceleration of every fluid parcel carrying the volatile solvent to depart from its respective trajectory; thence, a translational movement of the volatile substances towards the x-direction.

Method of Solution
Computing the system Equation (10) using the leading order Neumann and Dirichlet kinematic boundary conditions given in Section 2.2.1, one obtains To settle the spatial derivative ∂ x p appearing in u ( x, t), the pressure field model p (v) (t, x; ) is invoked into the pressure expression p (l) (t, x), giving where j ≥ 1. To augment the interfacial hydrodynamics by an −2 -order, we assume that ord (γ 0,l ) ∼ −2 . Consequently, this assertion suggests to rescale the surface tension of the liquid film by substituting −2 γ 0,l for γ 0,l . Thence, A similar approach was suggested by Oron et al. [12]. In the following section, the h-evolution equation is established.

Evolution Equation for the Interface
Using the general form of the Leibniz integral rule, one obtains When all these are considered, we deduce after integration the h-evolution equation, which writes In particular, it may well be pointed out here that when the hitherto neglected term ∂ x h∂ y c (l) , the induced pressure gradient ∂ x p (v) = K (v) , and the Navier-Slip coefficient β > 0 are disregarded into Equation (13), we deduce the well-known extended thin film model, written as in its general and dimensional form.

Preliminaries
Define three mobility terms, namely, M 1 (h), M 2 (h), and M 3 (h), respectively, as follows so that Equation (13) may be rewritten as such: Having settled the h-evolution equation in condensed form, we now set forth the latter into a weak differential formulation by re-defining the functions h, ∂ 2 x h as thus, h, ∂ 2 x h = h, h ; the idea of doing so is to help numerical implementation. This yields to the following h, h -system: which is now a degenerate second order nonlinear parabolic system. To adapt the situation, we construct at every point x of the surface Γ 2 and time t a depth average flow fieldˆ u = û x ,û y T with the velocity field u (t, x) as followŝ Substituting the vector field u ( x, t), established in Section 2.3.2, into the above integral yieldŝ Having laid enough foundations, we are now well-prepared to address the problem numerically; this is the duty of what follows.

Computing with the COMSOL Multiphysics Software
In summary, the flow of interest takes place in the vapor domain Ω v and the liquid domain Ω l with a strong coupling through the free surface Γ 4 among the fields c (l) (t, x), c (v) (t, x), and h (t, x). Let their respective systems of equations be designated by S c (v) , S c (v) , and (S h ), respectively. Then, for the evolution of the chemical concentration function c (l) (t, x) into Ω l , the corresponding system S c (l) writes S c (l) :D t c (l) − ∂ 2 y c (l) = 0 for all (t, x) ∈ I t × Ω l , where the hatted operatorD t (·) = ∂ t (·) +ˆ u · grad (·) designates the Lagrangian derivative operator subjected to the previously established depth averaged flow fieldˆ u = û x ,û y T . For the evolution of and, finally, for the evolution of the free-surface function h (t, x), the corresponding system writes where the scalar function F (h, h ), designating the flux functions of the thin film equation, is defined as Appendix A gives an account on Modeling Methodology using the COMSOL Multiphysics interfaces. With all these features settled conveniently into appropriate COMSOL Multiphysics interfaces, we ran several numerical simulations, whose solutions are discussed in the next section.

Results and Discussion
Herein, the results are discussed, emphasizing the hitherto neglected term ∂ x h∂ y γ and the hitherto overlooked term grad p (v) = K (v) (air-blow effect; K (v) , a prescribed constant).

Preliminaries
By an interfacial surface tension gradient is meant the vector gradient ∇ Γ 4 γ, where the operator grad Γ 4 (·) = (I − n ⊗ n) grad (·) is termed the Γ 4 -interface gradient operator on Γ 4 ; ⊗ is called the tensor product operator: for given vector fields v = v x , v y , w = w x , w y ∈ R 2 , their tensor product results in a 3 × 3 matrix with the entries v ⊗ w = v i w j (i,j)∈{x,y} 2 . Let grad Γ 4 be its characteristic scale. Then, grad Γ 4 (·) → grad Γ 4 grad Γ 4 (·). (Note: grad Γ 4 (·) is dimensional if it precedes → and dimensionless otherwise.) After nondimensionalizing the latter, this results in Further, asserting that ord grad Γ 4 ∼ L −1 -which is reasonable-yields at leading order. With these features beforehand together with the following definition and approximation, respectively, there follows that grad Γ 4 (γ) ≈ ∇ γ + ∇ ⊥ γ. Upon perusal of earlier works, one finds that the term grad Γ 4 (γ) was approximated as follows: grad Γ 4 (γ) ≈ ∇ γ. In other words, their authors considered thin film equations of the form For example, in the paper of O'Brien [8], the explicit form of H π; h, ∇ γ at leading order-which is the 0 -order treated in the present work-reads whereas, in the present work the following extended model is proposed, namely, where the property π stands as the set of control parameters, for instance, ours reads π def = x 0 , µ s , Ca, Ma ; and the differential δH goes to zero whenever ∇ ⊥ γ and ∇p (v) simultaneously go to zero. Thus, to distinguish our results with respect to others, it suffices to demonstrate through numerical experiments the following inequality: Graphically, the effects/intepretations of the terms ∇ γ = e x ∂ x γ, ∇ ⊥ γ = e x ∂ x h∂ y γ and ∇p (v) = K (v) e x , and δH can be discussed by the aid of two new concepts, namely, coarse, fine approximations.
Physically speaking, if an orthogonal projection is considered the vector quantity ∇ Γ 4 γ approximates to ∇ Γ 4 γ ≈ ∇ γ. This is the case of a coarse approximation of ∇ Γ 4 γ. The works of Matar et al. [7] and O'Brien [8] are, among others, two examples in which are found coarse approximation. Contrarily, if a rotational projection is adopted instead, this would result to ∇ Γ 4 γ ≈ ∇ γ + ∇ ⊥ γ. Thus, the h-evolution proposed here yields a better approximation, termed fine approximation. The left-hand side of Figure 2 illustrates those two concepts. In addition, when the free surface of the liquid film is subjected to an air-blow effect, ∇p (v) = K (v) (K (v) constant), the right-hand side of Figure 2 follows. Throughout what follows, by the statement steady-state regime/profile is meant that there exists a property t ∞ > 0 such that the property h (t, x) at every point x ∈ Γ 2 is invariant for all t ≥ t ∞ .  Figure 3.  Table 1.)

Case 1. The h-Evolution
x 0 Ca Ma 10 −3 10 −2 1.8 10 10 −2 (k, 5 ) L/4 0.5 500 Reference [7] Reference [16] Reference [8] It will be remarked in Figure 3 that in both cases the observed film thickness distribution is somewhat very similar, differing only in magnitudes. Consequently, this put forward the following fact: which demonstrates that the term ∇ ⊥ γ (t, x) does affect the solution of the thin film proplem. It is, therefore, of interest to examine in which situations the term ∇ ⊥ γ (t, x) is of negligible consequence, which is argued in what follows.

Of the Strength of the Source
What happens when we augment the strength µ s ≥ 0 of the source is here discussed. We carried out two parametric studies, one including and the other excluding the term ∇ ⊥ γ (t, x), with respect to µ s ranging from µ s = 10 to µ s = 41. In Figure 4 are exposed the graphical results which we shall now discuss.  [10,40] in the steady-state regime. These plottings are based on a zero pressure gradient basis, i.e., ∇p (v) = 0. (Note that the numerical experiments are carried out based on the data given in Table 1.) Glancing at Figure 4, one sees that the evolution of min h (µ s ) with respect to µ s is linear only when ∇ ⊥ γ ( x, t) = 0. So to speak, when ∇ ⊥ γ ( x, t) = 0 the quantity min h (t, x) decays nonlinearly. On this basis, if the nonlinear effect arises when ∇ ⊥ γ (t, x) is included in the mathematical model, obviously the rate at which the liquid film will dry or coat its substrate will be faster. Consequently, the higher be the effect of the volatile solvent, i.e., µ s 1, the more effective will be the effect of the term ∇ ⊥ γ (t, x); thus, this shows that its presence does affect the solution of the problem.

A Liquid Film of Infinitesimal Thickness
It be interesting to see how fast the quantity min h (t, x) goes to zero. In this section, upon the following substitution, 0.15 h 0 for h 0 , we consider an extreme case, a liquid film of infinitesimal thickness, to describe these situations. With these features beforehand, two simulations were ran, of which one takes account of the term ∇ ⊥ γ (t, x) and the other disregards it. The resulting superimposed curves are shown in Figure 5. Of these curves, the lower curve-the one such that min h (t, x) = 0 (i.e., h [ ∇ ⊥ γ =0] )-describes the case for which ∇ ⊥ γ (t, x) = 0; the upper curve-the other one such that min h (t, x) = 0 (i.e., h [ ∇ ⊥ γ=0] )-describes the case for which ∇ ⊥ γ (t, x) = 0. Physically speaking, it would appear that if the liquid film is of infinitesimal thickness, the higher the magnitude of the diffusive flux, the larger will be the deformation of the liquid film.
For that matter of comparison, cases for which ∇ ⊥ γ = 0 and ∇ ⊥ γ = 0 are exposed in Figure 6. As can be easily seen in Figure 6, the curve h ∇ ⊥ γ =0 goes to zero faster than h ∇ ⊥ γ=0 . Figure 6. Plots of the function min h (t, x) versus time when the numerical experiments includes (circular markers, •) and excludes (square markers, ) the term ∇ ⊥ γ; solid and dotted lines are approximations based upon the linear expression h (t, x = L/2) = −π 1 t + π 2 . (We note that this expression does not hold when t goes to zero, say, when t is within the range [0, 0.5), for min h (t = 0, x) = h 0 , as can be easily seen.) Numerically, we found in the linear regime that (only in the linear regime).
At this point, we thought it desirable to ascertain how long does min h (t, x) take to goes to zero. Based on the foregoing statements, it would also appear that the time at which min h [ ∇ ⊥ γ =0] ≈ 0 is t ≈ 2.5 and that at which min h [ ∇ ⊥ γ =0] ≈ 0 is t ≈ 1.75 -so to speak, a discrepancy of about 30 %. Recall the characteristic scales of t and h (t, x), namely, t ∼ −1 µL/ (γ 0,l − γ 0,v ) and h ∼ L, respectively, the dimensional form of h (t, x = L/2) = −π 1 · t + π 2 writes In particular, if τ ref > 0 designates a reference time at which min h (t, x) vanishes, it follows that thence, knowing π 1 (= min ∂ t h), π 2 (≈ αh 0 ; α ≤ 1, a correcting factor fitting the above linear law), respectively. Thus, an error of one second in the time at which min h (x, t) ≈ 0 will account for the term ∇ ⊥ γ. (Of these numerical values, the values of the physical property constants γ 0,l , γ 0,v , and µ are those found in the paper of [8].)

On Capillary and Marangoni Effects
The impact of Capillary (Ca) and Marangoni (Ma) numbers on the behavior of the thin film equation is now discussed. The resulting graphical results are shown in Figures 7-9, respectively. On inspecting those figures, it is easily seen that, on augmenting the number Ca the liquid film thins (Figure 8), while on the contrary, on incrementing the number Ma, it thickens ( Figure 9). Thus, out of these situations, it follows that, by careful variations of the control parameters Ca, Ma, respectively, one can regulate the thinning rate. Figure 7. Graphs of the functions h (t, x) with respect to the Marangoni number Ma in the steady-rate regime. The term ∇ ⊥ γ is set to zero in one case. (Note the numerical experiment is carried out based on the data given in Table 1.) As can been easily seen, Figure 7 shows that the evolution of the quantity min h (Ma) thins almost linearly with respect to Ma when ∇ ⊥ γ (t, x) = 0. But, since this is not the case when such term is taken into account, there obviously follows that the higher be the quantity Ma the more effective be consequence arising from ∇ ⊥ γ (t, x). Clearly, in the context of so-called Marangoni Drying, it acts as a additive effect, and, therefore, should not be left out of account.
Since the Capillary number Ca is a measure of the relative effect of viscous effects with respect to surface tension forces acting across the interface Γ 4 , an interesting way of interpreting this number is as follows: Let F µ , F γ denote the viscous and surface tension forces. Write Ca ref = F µ / F γ to designate a reference Capillary number. Then, variations of these forces by δ F µ , δ F γ results in the following strict inequalities The resulting curves are exposed in Figure 8 for On the other hand, the Marangoni number Ma may be regarded as that dimensionless quantity which describes the relative effect of chemical surface tension forces with respect to viscous forces acting across the interface Γ 4 . Thus, if the quantity F γ( c (l) ) designates the chemical surface tension force, and δ F γ( c (l) ) an increment of it, the following strict inequalities hold In the above condition, the quantity Ma ref = For example, an air blow over the liquid film is an instance which reflects this situation and is discussed here. In Figure 10 is displayed such a situation where ∇p (v) ≡ K (v) = 0.1. It shows how the shapes of the function h (t, x) describe themselves for subsequent instants of time when subjected to such a constant pressure gradient.  Table 1.) From these, it is observed that they all retain their sinusoidal forms and, furthermore, resolve themselves in a more or less regular manner into a deformed travelling waves, as can be easily seen in Figure 10. It is interesting to see how the interface h (t, x) evolves when we increase the quantity K (v) > 0 by 0.1, the resulting quantity by the same amount, and so forth. In other words, we wish to see what happens when we augment the strength of the air blow. Physically, the following statements do hold. Increasing the quantity ∇p (v) will both weaken and deviate the diffusive fluxes hitting the interface Γ 4 . Consequently, min |h (t, x) − h 0 | with a constant pressure gradient amounting to K (v) = 0.1 would be less than that which would amount to K (v) = 1.0, and thence the graphical results of Figure 10. (Note that the placement of the curved arrow appearing there is merely to aid theoretical understanding.)

Effect of an Air Blow with
We saw in earlier sections that including the term ∇p (v) > 0 results in translating the contaminated zone to the east, whereas including the term ∇ ⊥ γ results in thinning the liquid film faster. In this section, both terms are operative. Clearly, under this condition, the dynamics of the interface h (t, x) would be a combined translational movement accompanied by a fast thinning process. Indeed, glancing at Figure 11, one observes that this is clearly the case. However, while in Figure 10 the curved arrow is towards the north-east direction, that of Figure 11 is towards the south-east. We found a simple explanation to explain this mechanism. In fact, if for a prescribed value of ∇p (v) = K (v) = 1.0 we additionally take into account the term ∇ ⊥ γ = ∂ x h∂ y γ, the resulting effect would be, more diffusive fluxes will strike the interface Γ 4 while undergoing motion. Consequently, when time evolves, that region of the interface Γ 4 which is contaminated will move to the east while thinning the thickness of the liquid film to the south; thence, the curved arrow shown in Figure 11 is towards the south-east. If the gradient ∇p is disregarded in Equation (14), and the latter expressed in dimensionless form, a γ-dominated h-evolution equation is obtained, whose mathematical description writes This equation is similar to that leading order model proposed by O'Brien [8]; it is identical when When those two orders of magnitudes are proved negligible, one says that the solutions of the variable fields in question give agreement with respect to themselves. Consequently, in order to show such agreement, it only necessitates to display by means of numerical experiments graphical results elucidating the relative O'Brien (t, x) and of h (t, x) with respect to h O'Brien (t, x). In Figure 12 is portrayed the graph of (h − h O'Brien ) (x, t), and in Figures 13 and 14, respectively, are portrayed the superposition of the graphs C (l) (t, x) and C Exploring the discussion a step further, one deduces from Figure 12 the following inequalities: O'Brien (t,  Figure 15 illustrates the distribution of chemical concentration c (l) (t, x) into the liquid film Ω l when the ∇ ⊥ γ = 0 (left-hand side), and when ∇ ⊥ γ = 0 (right-hand side). The Marangoni flow induced in Ω l by the presence of the volatile solvent is illustrated in Figure 16 for the case ∇ ⊥ γ = 0. In the next section, we conclude the work.

Conclusions
We have studied the movement of a liquid film in an endeavor to examine in which situation does the hitherto neglected term ∇ ⊥ γ = ∂ x h∂ y γ have an effect, and a constant pressure-gradient-driven flow ∇p (v) = K (v) promotes the equilibrium thickness of the liquid film. On trial, the means of several numerical experiments gave graphs which do not have the same shapes. When ∇ ⊥ γ = 0, our results showed that the thinning process produced by the presence of the volatile source is much increased by the inclusion of the term-a phenomenon evidently occasioned in Marangoni drying. Thus, if a liquid film is used to realise a specific functional property relative to its contact surface, asserting that the term ∇ ⊥ γ is negligibly small could in some cases have a profound effect upon that property, and, therefore, can led to inaccurate predictions. When ∇p (v) = K (v) is considered, the liquid film thins less rapidly but undergoes translational movement in the direction of the pressure-gradient-driven flow so as to deviate the diffusive fluxes from hitting the free surface of the liquid film horizontally. Although there are applications in coating and drying processes when thin film flow instability resulting from a pressure-gradient-driven flow is undesirable, it is to be noted that the latter is desirable insofar that other transport phenomena come into play as, for instance, when heat transfer is concerned. Furthermore, if the effect of an air blow is too stiff and the thickness of the liquid film quite thin, instabilities result. On the whole, we have expounded conclusively at the following points:

•
The discrepancy between the several results is owing to the effect of the hitherto neglected term ∇ ⊥ γ on the dyamics of the liquid film.

•
The inclusion of a constant pressure-gradient-driven flow ∇p (v) = K (v) in the vapor phase domain might be of some advantage in supporting the thinning process.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Modeling Methodology with COMSOL
The modeling methodology is herein explained. For illustration, Figure A1 is considered.  Figure A1. An overview of numerical method. The COMSOL File includes two models, Mod-k (with k = 1, 2). Of these, Mod-1 is a 2D-Model containing two Physics Interfaces for the scalar fields c (j) (with j = 1, 2); Mod-2 is a 1D-Model comprising only one Physics Interface, that computes the h-evolution equation. Additionally, each models includes a coupling operator Opt-k (·) (with k = 1, 2) which passes data from one model.
Let those models be Mod-1, Mod-2, respectively. With this view in mind the methods of solutions are based on the following sequences of operations. In the first place, add two models to the COMSOL model builder, one two dimensional model (i.e., Mod-1) to solve problem statements S c (l) , S c (v) , respectively, and one one dimensional model (i.e., Mod-2) to solve problem statement (S h ). In the second place, add two Coefficient Form PDE interfaces to the two dimensional model and to the one dimensional model one Coefficient Form PDE interface. To construct the so-called 1D, 2D computational geometries, we proceed as follows. Under Mod-1 construct a large rectangular box whose length from (x, y) = (0, 0) spans L = 1 to the east and h 0 + L = 2 to the north. To describe the free-surface of the liquid film, a horizontal line defined by the equation y (x) = h 0 is inserted. (we note that the respective scalings L, h 0 are all dimensionless quantities, but, however, they are here employed to reflect the nature of the thin film.) In like manner, under Mod-2 construct a line segment whose length from x = 0 spans L = 1 to the east. We pointed out that the fields c (l) (t, x), c (v) (t, x) and h (t, x) are coupled. Consequently, to ensure these couplings from Mod-1 to Mod-2, we add a coupling operator Opt-1 to Mod-1 and another coupling operator Opt-2 to Mod-2. Logically, the inclusion of these operators Opt-1, Opt-2 is merely to pass data from a (source) domain to a destination domain. For the operator Opt-1, we let the interface Γ 2 be its source domain and the line segment defined under the model Mod-2 be its destination domain. For the operator Opt-2, we reverse the situation, letting the interface Γ 2 be its destination domain and the line segment defined under the model Mod-2 be its source domain. Thus, to invoke the functions c (l) (t, x, 0) of Mod-1 into Mod-2 and h (t, x) of Mod-2 into Mod-1, we shall use respectively. (More precisely, we used two linear extrusion operators, since, we wish to force out (that is, extrude) the source domain of definition of its functions onto an equivalent destination domain; see the COMSOL Reference Guide for further details.) To implement (S h ), we switch the coefficient form PDE interface of model Mod-2 to a two-dimensional system h, h -system so to speak we set Ξ = h, h T . Note that for all the problem statements S c (l) , S c (v) , and (S h ), the discretization method for computing the solution is quadratic. In Table 1 is exposed all the scalar quantities used, of which, some are partly extracted from the papers of O'Brien, Matar et al. and Buckingham et al. [7,8,16], while others are based on the underlying philosophy of the problem at hand. These are entered in the parameters table, found under the COMSOL global definitions interface. We made no mention of the imposed initial and boundary conditions. Their implementation are immediate, for, it suffices to consider the default and added nodes to implement them.