Marangoni stability of a thin liquid film falling down above or below an inclined thick wall with slip

The linear and nonlinear instability of a thin liquid film flowing down above or below (Rayleigh-Taylor instability) an inclined thick wall with finite thermal conductivity are investigated in the presence of slip at the wall-liquid interface. A nonlinear evolution equation for the free surface deformation is obtained under the lubrication approximation. The curves of linear growth rate, maximum growth rate and critical Marangoni number are calculated. When the film flows below the wall it will be subjected to destabilizing and stabilizing Marangoni numbers. It is found that from the point of view of the linear growth rate the flow destabilizes with slip in a wavenumber range. However slip stabilizes for larger wavenumbers up to the critical (cutoff) wavenumber. From the point of view of the maximum growth rate flow slip may stabilize or destabilize increasing the slip parameter depending on the magnitude of the Marangoni and Galilei numbers. Explicit formulas were derived for the intersections (the wavenumber for the growth rate and the Marangoni number for the maximum growth rate) where slip changes its stabilizing and destabilizing properties. From the numerical solution of the nonlinear evolution equation of the free surface profiles, it is found that slip may suppress or stimulate the appearance of subharmonics depending on the magnitudes of the selected parameters. In the same way, it is found that slip may increase or decrease the nonlinear amplitude of the free surface deformation. The effect of the thickness and finite thermal conductivity of the wall is also investigated.

Coating of surfaces is a very important activity with many industrial applications.The cooling of surfaces has also relevant applications in systems with heat generation.When the thin liquid films are falling down walls applications are found when coating or cooling surfaces where a temperature gradient exists across the film.The research on wall cooling has relevance in microchips located inside modern computers.Along with the advances of technology the microchips generate more and more heat that should be dissipated efficiently.Due to the usefulness of this subject extensive reviews were written some years ago by Oron et al. [1] and Dávalos-Orozco [2].More recently, new advances in the problem of non-isothermal thin film flows were reviewed in Dávalos-Orozco [3].Thin liquid films falling down walls have been investigated since many years ago under different boundary conditions.The stability of these flows is of interest because, when the film is not able to wet the wall, rivulets [4] may form and for some applications the film becomes useless.Usually it is difficult to control atmospheric temperature variations.Therefore, it is of interest to impose thermal boundary conditions to the wall-film system to solve the thermocapillary problem.The Marangoni instability has been investigated since many years ago.Pearson [5] studied the case of a flat free surface.Free surface deformation was first taken into account by Scriven and Sternling [6].The effect of gravity was included by Takashima [7,8].This instability was considered by Dávalos-Orozco and You [9] and Moctezuma-Sánchez and Dávalos-Orozco [10] in a cylindrical wall.In [11] it is of interest to apply a horizontal temperature gradient.This is also relevant in annular pools [12][13][14].
The no-slip condition at the wall surface is commonly used.However, it is possible to find surfaces where this condition is not satisfied [31,32].In some cases, the characteristics of the surface can be modified by changing its chemical composition [33].When the wall has very small topography, the liquid trapped inside plays the role of lubricant when the liquid is in motion reducing the effective friction and presenting an apparent slip at the wall [34].
Thermocapillary effects were taken into account along with slip at the wall in Kowal et al. [35].In the case of thin films falling down walls slip effect was considered in a number of works [36][37][38].Slip was also contemplated in the presence of topography by Usha and Anjalaiah [39] and Zakaria and Selim [40].The stability of a viscoelastic film was examined in Pal and Samanta [41] and Chattopadhyay et al. [42].
The finite thickness and thermal conductivity of the wall were not taken into account in the papers with slip and thermocapillary effects [35,36,38].Therefore, in this paper the linear and nonlinear thermocapillary instability of a liquid layer falling down a thick wall with finite thermal conductivity is investigated.The liquid layer may be located above or below (Rayleigh-Taylor instability [43,44]) the thick wall.It important to point out here that Karapetsas and Mitsoulis [45] have shown that the slip condition is necessary to have good agreement with experiments in newtonian and non-newtonian non-isothermal creeping flows.
In a previous paper [46], the Marangoni stability of a thin liquid layer coating a horizontal thick wall in the presence of gravity was investigated when slip is present at the interface between the wall and the liquid.There, it was found the possibility of a stabilizing effect of slip under different conditions.In some of the papers mentioned above stabilizing effects of slip have been mentioned.However, no definite analytical results have been given to determine the limits between stabilizing and destabilizing behavior of slip.Therefore, in this paper it is of interest to find out analytically and numerically the limits of the parameters in which slip may stabilize or destabilized when a main flow is present in a thin film flowing down above or below a hot or cold thick wall.
The structure of the paper is as follows.In the next section, the basic equations are presented along with the steps followed to obtain the nonlinear evolution equation of the free surface deformation under the lubrication approximation.Then, the linear results are presented in section three.The solutions of the nonlinear evolution equation are given in section four and section five are the conclusions.

Equations of motion and evolution equation
The sketch of the system under investigation is shown in Fig. 1.Notice that gravity is in the x-direction and that the figure sketches the particular case of a vertical wall.However, in this paper also the case of a wall inclined at 45 • degrees with respect to the hori- zontal is considered.Only when the wall is inclined 45 • degrees gravity may have positive or negative (Rayleigh-Taylor) direction.The mean liquid layer thickness is d f and the wall thickness is d W .The free surface height h(x, y, t), the unknown of the problem, is represented by a solid line.The origin is located at the interface between the liquid and the slippery thick wall.The atmospheres below the wall and above the free surface have different temperatures The equations of motion, energy and continuity of the liquid, the energy equation of the wall and their corresponding boundary conditions are made non dimensional as follows.The distance in the z-direction is Fig. 1 Sketch of the system.The liquid may slip on the wall falling down in the x-direction.The dashed line located at z = d f shows the mean thickness of the liquid layer.The wall has thickness d w .In the sketch the wall is vertical and parallel to the gravity acceleration vector g .However, in the paper an inclination angle of 45 • is also considered in two cases: g > 0 and g < 0 (Rayleigh-Taylor instability).The last condition corresponds to a thin film falling down below a thick wall.Two situations are considered with respect to the temperatures: T A0 > T A1 and T A0 < T A1 Vol:.( 1234567890) measured with d f the thickness of the liquid layer.The distance in the x-and y-direction by ∕2 = d f ∕ , where the wavelength is a representative longitude of the free surface deformation and the scaling parameter = 2 d f ∕ , is assumed to be  < 1 in the small wave- number (long wavelength) approximation.The time is scaled with d f ( ∕2 )∕ = d 2 f ∕ , velocity with ∕d f and pressure with 2 / d 2 f where and are the density and the kinematic viscosity of the liquid, respectively.The temperature is scaled with ΔT = T A0 -T A1 , where T A0 is the temperature of the atmosphere close to the wall and T A1 is the temperature of the atmosphere close to the liquid free surface.These scalings mean that the variation in the x and y-directions is different from that in the z-direction.This will have consequences in the scaling of the equations of motion and boundary conditions.
The nonlinear evolution equation of the unknown h(x, y, t) is derived by means of an asymptotic method using the following expansions of the dependent variables in terms of .They are The asymptotic method consists in substituting this expansions in the following equations and boundary conditions and solving recursively each resulting linearized set of equations.The scaled non dimensional equations of motion, continuity and energy of the liquid are: (1) (3) (5) The energy equation of the temperature T w of the solid wall is: where G = gd 3 f ∕ 2 is the Galilei number, q is the angle of inclination of the thick wall and Pr = ∕ is the Prandtl number.Here, is the thermal diffusivity of the liquid, c p is the heat capacity of the liquid and w and c pw are the density and heat capacity of the wall, respectively.is the ratio of the heat conductivity of the wall over that of the liquid.It is important to point out that the series expansion Eq. 1 is convergent in the present approximation when the Galilei number is order one or smaller.
The boundary conditions for the velocity and temperature of the liquid are the following.The slip boundary conditions (see [31]) and the condition of impenetrability at the wall are where = b∕d f is the non dimensional slip length of the wall and b is the longitudinal slip length.The no-slip boundary condition is attained setting = 0 at the wall.The ideal slip (friction free surface) limit is reached when → ∞ and u z = v z = 0 at the wall.Notice that the dimensional slip length b may have a range from 10 −9 m up to 10 −4 m for a superhydropho- bic substrate [47].However, there is also a variety of thin liquid film thicknesses from 10 −7 m (see [48]) up to 10 −4 m (see [49]) or more.Thus, the non dimen- sional parameter also has a wide range of magnitudes.Here, the selection of the magnitude of was made to show clearly the behavior of slip found in the analytical results.
The normal stress boundary condition at the free surface is (6) Pr T t + uT x + vT y + wT z = 2 T xx + 2 T yy + T zz .
where S = d f 2 ∕3 2 is the surface tension number, is the surface tension and and that H(x, y, t) is the free surface deformation.Besides, P A is the pressure of the atmosphere outside the liquid.The balance of shear stresses and surface tension temperature gradients on the free surface are the following.The first shear stress boundary condition is: and the second shear stress boundary condition is: where Ma = Ma * ∕Pr and Ma * = −(d ∕dT)ΔTd f ∕ is the usual Marangoni number and is the surface tension.The kinematic boundary condition at the free surface is The boundary conditions of the temperature in the liquid and the wall are as follows.At the free surface is the Biot num- ber and Q h is the heat transfer coefficient at the free surface.
As explained above, the expansions Eq. 1 are substituted into the equations and corresponding boundary conditions and the solutions are calculated recursively at the different orders of .Thus, from the (10) equations of the temperatures, the basic temperature profiles which satisfy the boundary conditions are: The denominator is defined as which is in fact a function of x, y and t through the free surface deformation h(x, y, t).It is very important because it includes the ratio d∕ , which is able to control the stability of the flow, as will be shown presently.Notice that in the lubrication approximation, instead of evaluating the free surface at a specific magnitude of the z-coordinate, the evaluation is made at h(x, y, t), the unknown location of the free surface deformation.As will be demonstrated below, h(x, y, t) will satisfy an evolution equation derived using the kinematic boundary condition.The solution of the pressure is where ∇ ⟂ = ( ∕ x, ∕ y) is the horizontal gradi- ent.It is important to point out that under the present approximation, it is enough to have the main temperatures and pressure profiles given in Eqs. 18, 19 and 20.
The solution of the x and y-components of the velocity at the lowest order, that is, the main velocity profiles, are: The continuity equation is used to obtain the vertical component of velocity at the lowest order.That is: Substitution of these results of Eqs.21, 22 and 23 into the kinematic boundary condition leads, at this order, to the equation: This approximation is used to substitute h t in the places where it appears in the following steps of the asymptotic procedure.The expressions of the solutions of the velocity components in the next order are very large and will be presented in the Appendix.However, it is important to mention the form of the following factors which appear in the thermocapillary terms: After substitution of these results into the kinematic boundary condition Eq. 13, it is possible to obtain the nonlinear evolution equation satisfied by the free surface perturbations in the presence of slip ( ) at the wall.That is: Notice that when q = 0 , this evolution equation reduces to that of Dávalos-Orozco and Sánchez-Barrera [46].In the same way, but with no-slip at the wall, this evolution equation reduces to that of Oron et al. [50] (in the absence of evaporation), to that of Oron [51] and to that of Podolny et al. [52] (when the (24) liquid is not binary).When q > 0 and Ma = 0 , Eq. 25 reduces to Eq. 3.7 of Samanta et al. [53].
In the next section, the linear stability is analyzed using this evolution equation.Then, in the following section, the nonlinear stability is investigated.

Linear stability
The linear instability is investigated linearizing the evolution Eq. 25 and assuming normal modes for the free surface deformation H(x, y, t) = A h exp(ikx + Ωt) found in h(x, y, t) = 1 + H(x, y, t) .Here k is the x-component of the wavenumber vector.Besides, Ω = Γ + i where Γ is the growth rate and is the frequency of oscillation.A h is a constant amplitude of the free surface perturbation.Substitution leads to the dispersion relation for Ω .The result is separated into real and imaginary parts.The frequency of oscillation is: = −G sin(q)(1 + 2 )k .Clearly, the frequency, and therefore the phase velocity, increases with .
The growth rate is obtained from the real part.That is: where L = 1 + Bi(1 + d∕ ) , is the linearized denom- inator (x, y, t) .As can be seen, L includes the ratio of thicknesses ratio d over the heat conductivities ratio in one parameter d∕ , which will play an important role on the stability through thermocapillary effects.The critical wavenumber is obtained from the condition that the growth rate Γ = 0. Two solutions are found.One corresponds to k c = 0 and the other one to (26) Vol.: (0123456789) The wavenumber corresponding to the maximum growth rate is obtained making zero the k derivative of the growth rate.
Notice the relation with respect to the critical wavenumber k c .The maximum growth rate is obtained after substitution of k max into Eq.26.That is As can be seen, Γ max in Eq. 29 represents a parabola in a plot Γ max vs Ma.That is: The focal distance p f is and the vertex is located at: (28) This means that at Ma Γ max =0 the instability is neu- tral against any perturbations.Notice here that both the width of the parabola (see p f in Eq. 31) and the location of the vertex through Ma Γ max =0 are controlled by the ratio d∕ and .This has important conse- quences on the results presented below.It is of interest to mention that when q = 0 , Ma Γ max =0 reduces to that found in [46].In that case Ma Γ max =0 could be negative but here it could be positive.This may have major consequences, mainly in the behavior of the stability with respect to the slip factor and the intersection of the curves of Γ max .
Next, the instability will be discussed based on the analytical results found above.First the case G > 0 will be discussed for the two conditions q = 90 • and 45 • .Then, the case G < 0 (Rayleigh- Taylor instability) will be discussed for the condition q = 45 • .The flow will also be subjected to a destabilizing Ma > 0 and a stabilizing Ma < 0 .Here it is assumed that the change in sign of Ma is due to a change in direction of the temperature gradient across the whole system.The number of parameters is large and therefore it is important to say that two parameters will be fixed through all the numerical calculations.They are Bi = 0.1 and S = 1.

Flow down above a thick wall: G > 0
In Fig. 2 results of two magnitudes of G = 0.1 and 0.3, two magnitudes of d∕ , and two inclination angles q = 90 • and 45 • are presented.The magnitude G = 0.1 was selected with the goal to find intersec- tions of the curves of Γ when increases (see Fig. 2a  and c).As can be seen, for q = 90 • and an increase to G = 0.3 intersections are not possible (see Fig. 2b and  d).The figure also shows the important role played by the ratio d∕ .Figure 2c and d demonstrate how this ratio decreases the instability when increasing from 0.5 to 1.8.In all cases, Ma destabilizes the flow Fig. 2 Growth rate Γ vs k.G = 0.1 : Fig. 2a, c and e. G = 0.3 : Fig. 2b, d and f (notice the different scales).d∕ = 0.5 : Fig. 2a and b. d∕ = 1.8 : Fig. 2c and d.From Fig. 2a-d: q = 90 • and Ma = 1.3 (1), 1.5 (2), 1.7 (3).From Fig. 2e-f: q = 45 • and d∕ = 0.5 .Figure 2e: Ma = 0.8 (1), 1.2 (2), 1.6 (3). Figure 2f: Ma = 1.9 (1), 2.2 (2), 2.5 (3).= 0.0 (solid line), 0.1 Vol.: (0123456789) increasing the magnitude of Γ .Observe here that the solid curves correspond to the no-slip boundary condition and were plotted for the sake of comparison.The intersections found in Fig. 2a and c are very important because they show that an increase of stabilizes the flow after a particular wavenumber.As found for q = 45 • , in Fig. 2e and f, it is still possible to find intersections even for G = 0.3 .The wavenum- ber for intersection between two curves with different , lets say 1 and 2 , can be derived making equal the two corresponding expressions of Γ .The result, for any angle of inclination, is: This means that for a wavenumber k > k + , has a sta- bilizing effect.Notice that k + reduces to that found in [46] when q = 0 which is independent from .Furthermore, in order to have an intersection the radicand has to be positive.Beware that it is possible to have intersection for stable Γ < 0 .That is what is going on in Fig. 2b and d, as can be understood following the tendency and slope of the curves.Therefore, Eq. 33 should be used with care.This intersection was independent of in the previous paper [46].
For example, the intersection of the curves = 0 and 0.1 in Fig. 2a for q = 90 • , Ma = 1.3 and d∕ = 0.5 occurs at k + = 0.1897 .That between the curves = 0 and 0.3 occurs at k + = 0.1933 and that between the curves = 0.1 and 0.3 occurs at k + = 0.1951 .In the case of Fig. 2e for q = 45 • the intersections when Ma = 1.2 and the other param- eters are the same are as follows: the intersection of the curves = 0 and 0.1 occurs at k + = 0.0911 , that between the curves = 0 and 0.3 occurs at k + = 0.0948 and that between the curves = 0.1 and 0.3 occurs at k + = 0.0966 .They are very close but different.
From Eq. 33 notice that always stabilizes when This may be verified in the left column of Fig. 2.This is more easy to verify in Fig. 2e and f for q = 45 • , ( 33) where clearly there is a tendency to move the k + of intersection to the left, very close to the origin, by decreasing Ma.In contrast, always destabilizes when k + = k c .In this case, k c should be evaluated with either of the 's, 1 or 2 .
The critical wavenumber number is plotted against G in the left column of Fig. 3.The Maximum growth rate is plotted against Ma in the right column of Fig. 3. Two angles of inclination q = 90 • and 45 • were used.Figure 3a 3d is for d∕ = 0.5 and G = 0.3 .Clearly, in the first column of Fig. 3, for q = 90 • intersections of the curves of k c occur for cer- tain G.However, this is not the case when q = 45 • , where, in addition, k c decreases with G, for both mag- nitudes of d∕ .The intersections between curves of k c occur at the following G: Notice that the intersections for q = 45 • occur for G + < 0 , using the minus sign in Eq. 35.This case of Rayleigh-Taylor instability is discussed below.It is important to point out that the role of the ratio d∕ is to stabilize the flow when its magnitude increases.It also contributes to modify the point of intersection of the curves of k c through L in G + .The angle of incli- nation q also has important influence on G + .How- ever, the whole role of q on the stability is still more complex and it is discussed in different places below.
The Maximum growth rate in the right column of Fig. 3 shows that for q = 90 • the maximum growth rate Γ max increases monotonically with an increase of to the right of the vertex ( Ma Γ max =0 ) of the parabola.However, the contrary occurs to the left of the vertex.It is of interest that this behavior is reversed with respect to the vertex for q = 45 • .Therefore, from the point of view of the maximum growth rate also has a stabilizing effect at q = 45 • .This may be corroborated at the right column of Fig. 3 where the maximum growth rate is plotted against Ma for two angles of inclination q = 90 • and 45 • .Figure 3b corresponds to G = 0.1 and Fig. 3d to G = 0.3 .Notice how the vertex of the parabola Eq. 30 is displaced with the variation of the parameters (see Eq. 32).Moreover, the width of the ( 35) parabola (latus rectum) also changes as four times Eq. 31.The curves are plotted, to the left and to the right of the figure, up to the place where intersections are found.For q = 90 • , it is easy to observe two intersections.One is located just to the left of the vertex of the parabolas, for Ma < 0 .The second is located farther to the left.Apparently between these two intersections stabilizes.On the contrary, the intersections when q = 45 • are located just to the right of the vertex and farther away.Between these two intersections stabilizes.
It is important to be careful because the intersections of the parabolas made by Γ max for q = 90 • are misleading.The reason is that the intersections occur for Ma < Ma Γ max =0 where Ma < 0 is stabilizing.Then, some of the curves are stable ( Γ < 0 ) and therefore no intersections (that is Γ max ( 1 ) = Γ max ( 2 ) ) are pos- sible.This invalidates the meaning of the intersections found for q = 90 • to the left of the vertex due to the natural stabilizing ( Γ < 0 ) effect of large enough Ma < 0 , not observed in Fig. 3b and d.
The case of q = 45 • is interesting because the stabilizing behavior of in fact occurs between the two intersections.Again, to the left of the vertex it is found that the curves of Γ are stable, in the same way as for q = 90 • , even in the region 0 It is possible to show analytically that Γ < 0 (stable) at Ma Γ max =0 of the vertex.To this goal, it is enough to substitute Ma Γ max =0 of Eq. 32 into Γ in Eq. 26.It is found that at Ma Γ max =0 the growth rate = −Sk 4 (3 + 1) is always negative for k > 0 and any , G and inclination angle q.As a consequence, Γ < 0 for any magnitude of Ma ≤ Ma Γ max =0 .In par- ticular, the result is valid for q = 0 [46].The Ma for an intersection is found assuming that Γ max and Ma are the same for two different , lets say 1 and 2 .Two solutions are found and they are: where These two formulas Eqs.36 and 37 reduce to those derived in [46] when q = 0 • .They give the two limits of the regions of Ma where the behavior of changes from destabilizing to stabilizing from the view point of the maximum growth rate.As mentioned above, this intersections are valid and have physical meaning only for Ma > Ma Γ max =0 .
Therefore, it is of interest to understand in a deeper way the behavior of the vertex of the parabola that is responsible of the intersections among the curves of different .To attain this goal first assume that Ma Γ max =0 is a function of and that it can be written as Ma Γ max =0 ( ) .Then assume that  2 >  1 .Subtract Ma Γ max =0 ( 2 ) -Ma Γ max =0 ( 1 ) .The result has a number of positive factors, but the product of some other factors of the numerator determines if the difference is positive or negative.That part of the numerator is: where is always positive.Notice that by assumption  1 −  2 < 0 and therefore the term inside brackets determines the sign of the difference between the ( 36) two Ma's of the vertexes.To understand the behavior of the vertexes it is possible to calculate two critical values useful to determine the change of sign of Ma Γ max =0 ( 2 ) -Ma Γ max =0 ( 1 ) .One critical value is: which makes zero the numerator.It is also possible to derive a critical cos(q) by substitution of sin(q) 2 = 1 − cos(q) 2 and look at the change in sign from the point of view of the angle of inclination.The result is When G < G * , or else for q < q * , from the difference it is concluded that , that is, the vertex Ma Γ max =0 ( 2) is located to the right of the vertex Ma Γ max =0 ( 1 ) .In case G > G * , or else, for q > q * the contrary occurs and Ma Γ max =0 ( 2) is located to the left of Ma Γ max =0 ( 1 ) .This is also the reason for the displacement of the intersections among the curves of different 's.They change from left to right of the vertex depending on the magnitude of G (or q) with respect to the critical G * (or q * ).In fact, the criti- cal G (or q) of the subtraction of the Marangoni number of the intersections of the curves of maximum growth rate, that is Ma ++ -Ma −− , is exactly the same G * Eq. 38 (or q * Eq. 39).For G < G * , it is found that Ma ++ > Ma −− , that is, Ma ++ is located to the right of Ma −− .The contrary is valid when G > G * .Notice that in the following useful subtractions the critical is also G * .For the subtraction of Ma All these inequalities are reversed when G > G * .From all these inequalities, derived algebraically, it is concluded that: These inequalities are very informative, but to give them a more precise meaning it will be helpful to employ the two extremal cases q = 0 • (which cor- responds to Eq. 40) [46] and q = 90 • (which cor- responds to Eq. 41).If 2 > 1 , it can be shown analytically that for q = 0 • the vertex always has Ma Γ max =0 > 0 , and the intersections of the maximum growth rate are located at Ma ++ > 0 and Ma −− > 0 .
On the other extreme, for q = 90 • the vertex always has Ma Γ max =0 < 0, Ma ++ < 0 and Ma −− < 0. The change of these particular inequalities start to occur when G (or q) approaches and crosses the critical G * Eq. 38 (or q * Eq. 39).That is, in the process some of the terms of the inequalities in Eqs.40 and 41, not all, might be negative.Notice that in the case of G < 0 (Rayleigh-Taylor instability) the sign of the numerator is exclusively determined by 1 -2 < 0. Thus the numerator is negative.Therefore, as will be shown presently, Ma Γ max =0 ( 2) is located to the left of Ma Γ max =0 ( 1 ) .In this particular case, Ma ++ < Ma −− , that is, Ma −− is always located to the right of Ma ++ .Then, for q = 0 • the vertex always has Ma Γ max =0 < 0, and the intersections Ma ++ < 0 and Ma −− < 0 [46].For q = 90 • , that is, for a vertical wall, the results when G < 0 are of no interest because the phenomena are exactly the same as for G > 0.
3.2 Flow down below a thick wall at q = 45 • (Rayleigh-Taylor instability) G < 0: destabilizing Ma > 0 and stabilizing Ma < 0 From now on the angle of inclination is fixed at q = 45 • and the ratio is d∕ = 0.5 .The numeri- cal results for the destabilizing and stabilizing Ma are plotted in Fig. 4. The Fig. 4a and b correspond to a destabilizing Ma > 0 and Fig. 4c and d to stabi- lizing Ma < 0. Besides, Fig. 4a and c correspond to G = −0.1 and Fig. 4b and d correspond to G = −0.3 .
It is important to note the different scales of the growth rate of each sub-figure.Observe that for the magnitudes of the parameters used, the intersections of the curves, where changes its behavior from destabilizing to stabilizing, appear in Fig. 4a and b for Γ > 0 .The wavenumber k + for these intersections is (41) also predicted by the formula Eq. 33.No intersections were found for Γ > 0 in Fig. 4c and d.The critical wavenumber for the Rayleigh-Taylor instability is plotted against G in the first column of Fig. 5. Clearly, intersections of the curves are observed which are described by the formula in Eq. 35 of G + .In this case the minus sign in front of G + should be used.It is found that the wavenumber increases with the magnitude of G but the growth is different for each .
As a consequence, only the right branches of the parabolas are physical and therefore destabilizes from the point of view of the maximum growth rate.Notice that the curves are displaced to the left when d∕ increases from 0.5 to 1.8.
The maximum growth rate against Ma is plotted in the right column of Fig. 5.For the sake of visibility of the separation of the curves, only G = −0.1 and - 0.2 were used in the figure .The ratio d∕ = 0.5 cor- responds to Fig. 5b and d∕ = 1.8 corresponds to Fig. 5d.In each figure, the right hand intersection is found very close to the vertex of the parabolas.Therefore it is of interest to know if there is a possibility to have a real change of the stability behavior of at those particular intersections.To this goal, again use is made of the result which says that at Ma Γ max =0 the growth rate Γ = -Sk 4 (3 + 1) is always negative for any , G and inclination angle q.With this, it is possible to compare if any intersection occurs to the right of Ma Γ max =0 to be of interest.It is found, using Eqs. 36and 37, that the right hand intersection is always to the left of Ma Γ max =0 and therefore the intersections between maxima of the growth rate are not real.In fact, the curves of the growth rate are stable.
A discussion on the physical effects of slip is in order.As can be seen in Eqs. 25, 26 and 29 slip affects in a different way the Galilei number and the Marangoni number through their coefficients.Even more, the presence of gravity alone as a restoring force perpendicular to the wall and as promoting motion along the wall have different slip coefficients and therefore different slip effects which oppose each other from the stability point of view.Besides, when G and Ma oppose to each other a variation of changes the relative magnitude between these two parameters.As a consequence this has the stabilizing effect found in the growth rate Γ after certain wavenumber k + smaller than k c .In the case of the maximum growth rate Vol.: (0123456789) Γ max this effect occurs between a lower and an upper limit Marangoni numbers (which have to be positive to be valid intersections).Between these two limits the magnitude of leads Γ max to decrease when Ma increases.Therefore, the usual behavior of gravity and thermocapillary terms on stability is altered by slip through the modification of the coefficients of the different terms inside the equations.Notice that only the physical behavior of the phase velocity occurs as expected in the presence of slip.That is, the phase velocity just increases with .

Nonlinear free surface profiles
In this section the free surface profiles are analyzed using a normal mode expansion of the height h up to 5 modes.Seven modes were also investigated, but it was found that the free surface deformation profiles were almost the same.
The expansion is:  However, a nine normal modes expansion was needed in the case where a relatively large Marangoni number was used.This h(x, t) is substituted into Eq.25 to obtain a set of ten coupled nonlinear ordinary differential equations for A i (t) and Ac i (t) , i = 1, ⋯ , 5. The initial value is assumed to be A 1 (0) = Ac 1 (0) = 0.01 , which corresponds to a cosine initial perturbation of the form 2cos(kx).Each solution of A i (t) and Ac i (t) is substituted back into Eq.42 for h(x, t), which is plotted in the figures below.The goal is to have the free surface profile of the film flowing down a wall at a particular time.The dependance of the profiles on the parameter d∕ is also investigated.In general, the sample data were selected from the linear problem.In particular, the wavenumber corresponding to the maximum growth rate was selected for each plotted curve in the numerical analysis.The free surface perturbations are initially unstable and then, after a large enough time, they reach nonlinear saturation.

Nonlinear films flowing down above a thick wall
The nonlinear results of the free surface profiles of a thin film flowing down above a thick wall are presented in Fig. 6.These results are based on those of Fig. 2 for linear instability.For each curve, the wavenumber corresponding to the maximum growth rate was selected.Notice that the wall is not drawn for the sake of presentation of the free surface profiles which have a small amplitude in comparison with the non dimensional thickness of the liquid layer.It is important to remember that the phase velocity of the perturbations differs for each given G and .In Fig. 6 the first figures, from Fig. 6a to d, the angle of inclination is q = 90 • and for the last two figures, from Fig. 6e to f, q = 45 • .In the first set of figures the Marangoni number is fixed to Ma = 1.3 .However, Fig. 6a: G = 0.1 and d∕ = 0.5 , Fig. 6b: G = 0.3 and d∕ = 0.5 , Fig. 6c: G = 0.1 and d∕ = 1.8 , Fig. 6d: G = 0.3 and d∕ = 1.8 .For the second set of figures different Marangoni numbers are used.For Fig. 6e Ma = 0.8 with G = 0.1 and for Fig. 6f Ma = 1.9 with G = 0.3 .In general, the free surface profiles are plotted using the same lines of .That is, = 0.0 (solid line), 0.1 (dashed line), 0.3 (dotted line).The running time used for all the curves was 5000 units.
A comparison of Fig. 6a and c shows that and increase of d∕ decreases the waves amplitudes.The same can be said about the curves in Fig. 6b and d.However, this increase stimulates the growth of the subharmonics (see Fig. 6c and d).On the contrary, it is observed that works to decrease the amplitude of the subharmonics, as can be seen starting from the solid curve (no-slip) up to the dotted ( = 0.3 ) one.This may be considered a nonlinear stabilizing effect of when the film flows down a vertical wall.However, in the Fig. 6a-d the amplitudes increase with , which is a destabilizing effect.
When the wall is inclined at q = 45 • and d∕ = 0.5 , Fig. 6e shows no notable change when increases.However, Fig. 6f presents an interesting picture.Here, 9 normal modes were used to describe the free surface profile.It is observed that an increase of stimulates the increase of the amplitude of the subharmonics when going on from the solid curve (no-slip) up to the dotted ( = 0.3 ) one.This may be considered a nonlinear destabilizing effect of when the film flows down above an inclined wall.Observe that in Fig. 6e decreases the amplitude of the free surface perturbation but in Fig. 6d increases the amplitude of the surface deformation.

Nonlinear films flowing down below a thick wall:
Rayleigh-Taylor instability Here the curves are plotted below the x-axis of each sub figure of Fig. 7.That is because in the Rayleigh-Taylor instability h = 1 − H and the film is pulled in the direction parallel to gravity at the same time that it is flowing down and below the inclined thick wall.The angle of inclination is fixed to q = 45 • and d∕ = 0.5 .The Galilei number of Fig. 7a is G = 0.1 and that of Fig. 7b is G = 0.3 .However, both figures include a destabilizing effect of Ma = 1.3 .Besides, Fig. 7c is for G = 0.1 and Fig. 7d is for G = 0.3 , but both include a stabilizing effect of Ma = 0.3.In Fig. 7a, for a destabilizing Ma, it is observed that the harmonics are suppressed by the increase of .However, the amplitude increases with .It is interesting that an increase of the magnitude of G in Fig. 7b changes the behavior of .As can be seen, = 0.1 first stabilizes and then destabilizes with a notably larger amplitude for = 0.3 .Moreover, the appearence of sub harmonics is stimulated with .
In Fig. 7c, the Marangoni number Ma stabilizes and it is seen that stimulates the presence of sub harmonics but works to decrease the amplitude of the free surface deformations.With the increase of the magnitude of G in Fig. 7d, promotes the suppression of subharmonics but increases the amplitude of the free surface perturbation.It is important to note the different amplitudes among the curves in Fig. 7a-d.For the sake of comparison the vertical axes present the same scale length in all the sub figures.

Conclusions
In this paper the linear and nonlinear thermocapillary instability of a thin film flowing down above or below on a thick wall with slip was investigated.It is important to point out that the results about the stabilizing and destabilizing effects of slip are new.The finite thickness and thermal conductivity of the wall showed to have a relevant role on the stability.Depending on the magnitude of the parameters, the growth rate is able to stabilize after a certain wavenumber k + in the presence of slip.The curves of criticality may intersect among them at G + when increases.This intersections are a consequence of the stabilizing effect of slip on the growth rate.From the point of view of the maximum growth rate plotted against the Marangoni number, promising intersections among the curves brought about the possibility of a change in behavior of slip from destabilizing to stabilizing.A very detailed analysis was done from which it was concluded that a change of the stability behavior of slip is only possible to the right of the vertex of the parabola formed in a plot Γ max vs Ma.The reason is that, starting from Ma Γ max =0 , the growth rate Γ is already negative.For Ma < Ma Γ max =0 the growth rate is still more negative and there is no possibility of intersections of the maximum growth rate.The general result is that only intersections found to the right of the vertex  are physically realizable.It was found that only when the film flows down above an inclined thick wall and for angles of inclination of the wall smaller than those of the vertical (q < 90 • ), it is possible to find intersections of the curves of maximum growth rate to the right of the vertex.There were found intersections among the curves of maximum growth rate where the stability behavior of slip may change.
From the nonlinear point of view, for certain magnitudes of the parameters it was found that slip may suppress the appearance of sub harmonics and for other magnitudes slip stimulates their appearance.It was also observed that slip may increase or decrease the amplitude of the nonlinear free surface perturbations, depending of the magnitude of the parameters selected.