Didier Clamond
email: diderc@unice.fr

Denys Dutykh
email: denys.dutykh@ucd.ie

FAST ACCURATE COMPUTATION OF THE FULLY NONLINEAR SOLITARY SURFACE GRAVITY WAVES

Keywords: Solitary gravity waves, Fully nonlinear water wave equations, Petviashvili method

In this short note, we present an easy to implement and fast algorithm for the computation of the steady solitary gravity wave solution of the free surface Euler equations in irrotational motion. First, the problem is reformulated in a fixed domain using the conformal mapping technique. Second, the problem is reduced to a single equation for the free surface. Third, this equation is solved using Petviashvili's iterations together with pseudo-spectral discretisation. This method has a super-linear complexity, since the most demanding operations can be performed using a FFT algorithm. Moreover, when this algorithm is combined with the multi-precision floating point computations, the results can be obtained to any arbitrary accuracy.

Solitary waves play a central role in nonlinear sciences. They appear in various fields ranging from plasmas physics to hydrodynamics and nonlinear optics [10,22]. In some special cases, solitary wave solutions can be found analytically. For example, explicit expressions are known for integrable models such as KdV and NLS equations, but also for some non-integrable equations such as [4,[START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF]. However, no closed form solutions are known for the practically very important case of the free surface Euler equations. In order to construct these solutions, one has to apply some approximation methods. Historically, asymptotic approximations have been proposed first [6] but, these series being divergent [8], the high-order approximations are only valid in the limit a → 0 (a the wave amplitude). In order to avoid all these limitations, various numerical approaches have been also proposed, e.g. [7,15]. We can mention the boundary integral equation method applied to the computation of solitary waves, see [5,[START_REF] Hunter | Solitary and periodic gravity-capillary waves of finite amplitude[END_REF]13,[START_REF] Parau | Nonlinear three-dimensional gravity-capillary solitary waves[END_REF] and the references therein.

One of the most widely used methods nowadays is Tanaka's algorithm [START_REF] Tanaka | The stability of solitary waves[END_REF]. In this note, we propose a simpler and faster method based on Petviashvili's iterations [18] leading to a very efficient numerical scheme for the computation of solitary gravity waves of the full Euler equations in the water of finite depth. Petviashvili's iterations have been successfully used for various simplified water wave equations but, to our knowledge, never for the fully nonlinear problem. The algorithm proposed here is efficient and it allows a very straightforward compact implementation under Matlab, for example. Such a script can be freely downloaded [1].

The note is organised as follows. In section 2, we introduce the hypothesis, the notations and the Babenko equation which is equivalent to the original problem. In section 3, we describe a very simple and fast algorithm to compute the solution, and we give some numerical results. Finally, the main conclusions are outlined in Section 4.

Mathematical model

We consider a steady two-dimensional potential flow due to a solitary gravity wave propagating at the surface of a homogeneous fluid layer of constant depth d. Let be (x, y) a Cartesian coordinate system moving with the wave, x being the horizontal coordinate and y the upward vertical one. Since solitary waves are localised in space, the surface elevation tends to zero as x → ±∞, and x = 0 is the abscissa of the crest. The equations of the bottom, of the free surface and of the mean water level are given correspondingly by y = -d, y = η(x) and y = 0, respectively.

2.1. Governing equation. Let be φ, ψ, u and v the velocity potential, the stream function, the horizontal and vertical velocities, respectively, such that u = φ x = ψ y and v = φ y = -ψ x . It is convenient to introduce the complex potential f ≡ φ + iψ (with i 2 = -1) and the complex velocity w ≡ uiv that are holomorphic functions of z ≡ x + iy (i.e., w = df dz). The complex conjugate is denoted with a star (e.g., z * = x-iy), while subscripts 'b' denote the quantities written at the seabed -e.g., z b (x) = xid, φ b (x) = φ(x, y = -d) -and subscripts 's' denote the quantities written at the surface -e.g., z s (x) = x + iη(x), φ s (x) = φ(x, y =η(x)). We emphasise that ψ s and ψ b are constants because the surface and the bottom are streamlines.

The far field velocity is such that (u, v) → (-c, 0) as x → ±∞, so c is the wave phase velocity observed in the frame of reference where the fluid is at rest at infinity (c > 0 if the wave travels to the increasing x-direction). Note that c = (ψ bψ s ) d due to the mass conservation [12].

The dynamic condition can be expressed in term of the Bernoulli equation

2 p + 2 g y + u 2 + v 2 = c 2 , (2.1)
where p is the pressure divided by the density and g > 0 is the acceleration due to gravity. At the free surface the pressure equals that of the atmosphere which is constant and set to zero without loss of generality, i.e., p s = 0. 1). The conformal mapping yields the Cauchy-Riemann relations x α = y β and x β = -y α , while the complex velocity and the velocity components are

Conformal mapping. Consider the change of independent variable

z ↦ ζ ≡ (iψ s - f ) c, that conformally maps the fluid domain {-∞ ⩽ x ⩽ ∞; -d ⩽ y ⩽ η} into the strip {-∞ ⩽ α ⩽ ∞; -d ⩽ β ⩽ 0} where α ≡ Re(ζ) and β ≡ Im(ζ) (c.f. Fig.
w = -c d z dζ -1 , u = -c x α x 2 α + y 2 α , v = -c y α x 2 α + y 2 α , u 2 + v 2 = c 2 x 2 α + y 2 α . Introducing new dependent variables X(α, β) ≡ x -α and Y (α, β) ≡ y -β, the Cauchy- Riemann relations X α = Y β and X β = -Y α hold, while the bottom (β = -d) and the free surface (β = 0) impermeabilities yield Y b (α) ≡ Y (α, -d) = 0 and Y s (α) ≡ Y (α, 0) = η.
The functions X and Y can be expressed in term of X b solely using some pseudo-differential operators [3]

X(α, β) = cos[ (β + d)∂ α ] X b (α), Y (α, β) = sin[ (β + d)∂ α ] X b (α), (2.2)
so that the Cauchy-Riemann relations and the bottom impermeability are fulfilled identically. At the free surface β = 0, (2.2) yields X s (α) = cos[d∂ α ] X b (α), that can be inverted as X b (α) = sec[d∂ α ] X s (α), and hence the relation (2.2) yields

Y s (α) = η(α) = tan[ d ∂ α ] X s (α), (2.3)
which relates quantities written at the free surface only. The relation (2.3) can be trivially inverted giving, in particular,

(∂ α X) s = C {Y s } ≡ ∂ α cot[d∂ α ] Y s ,
where C is a self-adjoint pseudo-differential operator acting on a pure frequency as

C e ikα = k coth(kd) e ikα (k ≠ 0) 1 d (k = 0) (2.4)
This operator can be efficiently evaluated in the Fourier space using a FFT algorithm.

Babenko's equation.

A Lagrangian for water waves can be obtained considering the kinetic energy minus the potential one, leading, in the conformal plane, to the functional

L = ∫ ∞ -∞ L dα with L being the Lagrangian density [15] L = 1 2 c 2 η C {η} -1 2 g η 2 (1 + C {η}).
The governing equation is the Euler-Lagrange equation of this functional, i.e.

0 = ∂ L ∂η + C ∂ L ∂ C {η} = c 2 C {η} -g η -1 2 g C {η 2 } -g η C {η}, (2.5)
which is the Babenko equation for gravity solitary surface waves [START_REF] Babenko | Some remarks on the theory of surface waves of finite amplitude[END_REF]. Applying the operator C -1 to (2.5) and splitting the linear and nonlinear parts, the equation can be presented as

L {η} = N {η}, L {η} ≡ c 2 η -g C -1 {η}, N {η} ≡ g C -1 {η C {η}} + 1 2 g η 2 ,
(2.6) that is more convenient for the simple and efficient numerical resolution explained below.

Numerical method and results

First, we define the wave by choosing its Froude number F ≡ c √ gd. Second, we choose a KdV-like initial guess of the free surface elevation η 0 (α)

η 0 (α) ≡ d (F 2 -1) sech 2 (κα 2), (a) F 2 = tan(κd) κd (b). ( 3
η n+1 = S 2 n L -1 {N {η n }} , S n = ∫ ∞ -∞ η n L {η n } dα ∫ ∞ -∞ η n N {η n } dα , (3.2)
where S n is the co-called stabilisation factor which can be computed either in physical or Fourier spaces, thanks to the Parseval identity [START_REF] Titchmarsh | The Theory of Functions[END_REF]. In the implementation, we systematically privilege the Fourier space since the operators C and C -1 can be very efficiently computed according to the definition (2.4) using any FFT algorithm. The convergence of the iterative process (3.2) is checked by following the norm of the difference between two successive iterations along with the residual in the ℓ ∞ norm

η n+1 -η n ∞ < ε 1 , L {η n } -N {η n } ∞ < ε 2 ,
where ε n are some prescribed tolerance parameters that are usually of the order of the floating point arithmetics precision.

3.1. Numerical results. The present algorithm is illustrated on the following example of a fairly high amplitude solitary wave. It is validated and compared against Tanaka's algorithm [START_REF] Tanaka | The stability of solitary waves[END_REF]. Tanaka's solution is parametrized by the dimensionless parameter q c = u c √ gd ∈ (0, 1), where u c is the fluid velocity at the wave crest. We choose the special value q c = 0.4 for comparisons. Tanaka's algorithm was implemented in Matlab 1 in its original version without any peculiar optimisations. The interval -30 ⩽ φ ⩽ 30 is discretised using 20, 001 nodes and the iterations are continued until the tolerance 10 -10 is obtained between two successive values of the Froude number F , i.e. Then, we take this Froude number and use it as the solitary wave parameter in the Babenko equation (2.6). We use a (periodised) domain of the same length as in Tanaka's method, and discretise it with N = 16, 384 nodes. The tolerance is chosen to be 10 -15 , i.e. very close to machine accuracy. These parameters are kept in all the computations presented below. The iterative procedure is stopped when the ℓ ∞ norm of the difference between two successive iterations is less than the tolerance (the norm of the residual is checked a posteriori). Our algorithm required 372 iterations to converge and the computations lasted for about 0.5 s, i.e. almost three orders of magnitude faster than Tanaka's algorithm. This flagrant difference in CPU times can be explained by the fact that our method has the super-linear complexity O(N log N) while Tanaka's method has a quadratic complexity O(N 2 ).

Moreover, the present algorithm can be used to compute solitary waves to an arbitrary precision. To illustrate this, we implemented the scheme using a Multi-Precision Computing Toolbox (MPCT) for Matlab [14]. The code modifications needed to implement the arbitrary accuracy are minor; this constitutes one of the major advantages of this MPCT. In these higher-precision experiments, all floating point operations are done with 50 significant digits. The tolerance parameter is set to 10 -50 . To achieve such extreme accuracy, the Petviashvili scheme required 1, 386 iterations. The convergence of the algorithm is shown on Figure 2-b. One can notice that the residual error is slightly higher than the difference between two successive iterations, likely due to the presence of quadratic nonlinearities into the equation. However, both measures of the numerical error decrease monotonically with iterations. The resulting solitary wave profile is depicted on Figure 2-a along with the corresponding Tanaka solution. The two curves cannot be distinguished up to the graphical resolution, thus validating our method. The respective accuracies can be asserted comparing the dimensionless amplitudes a d (Table 1). One can notice that despite the tolerance parameter set to 10 -10 in Tanaka's algorithm, its real accuracy is seriously degraded for this large amplitude wave (only five digits of the wave amplitude are predicted correctly). Conversely, the present method gives fourteen correct digits with a tolerance of 10 -15 in double precision calculation (about sixteen digits), and 49 correct digits using the 50-digit computation (the latter was determined using a 60-digit calculation with a tolerance of 10 -55 with eight times more nodes over an interval twice the size).

The algorithm being validated, it can be used to produce some physically sound numerical results. For example, we compute to high-precision (here 10 -15 , but it is not a limitation) the speed-amplitude relation for solitary waves. The result is shown on Figure 3. We also computed several values in a wide range of amplitudes using Tanaka's algorithm (represented with circles). Again, up to the graphical resolution we have a very good agreement between the two methods. However, the speed of Tanaka's algorithm does not allow to compute many values in a reasonable time.

We note finally that, for simplicity, we defined the wave by its Froude number which varies in the interval 1 ⩽ F ≲ 1.29421. Doing so, we can compute all waves such that 0 ⩽ a ≲ 0.796 that is sufficient for all cases of practical interest. One of the largest solitary waves which can be predicted by this parametrization is represented on Figure 4. However, the solutions of larger amplitudes exist, the limit a ≲ 0.796 being obtained here because F is not a monotonic function of a (i.e., F (a) is multi-valued when a is large). Thus, in order to compute larger waves, another parameter (e.g., a d or q c ) should be used to define the wave. This requires simple modifications in the code [1]. We successfully did it and were able to compute all waves such that a ≲ 0.82. For even higher waves, the present method is not efficient because these extreme waves being (almost) singular at the crest, the Fourier series does not converge fast enough and another spectral decomposition should be used instead.

Conclusions

We have described a fast, accurate and easy to implement method for computing steady gravity solitary surface waves solution of the full Euler equations. The method is based on first rewriting the problem into the equivalent but simpler Babenko's equation and, second, using Petviashvili's iterations to solve it numerically. The resulting scheme allows the computation to any arbitrary accuracy using a multi-precision floating-point arithmetics library [14]. All ingredients of the method are well known, but their combination seems to be new in this context and turns out to be very efficient. Our Matlab implementation is rather compact: The computational numerical core is about a dozen lines of code, the rest being pre-and post-processing. The script [1] can be freely downloaded and used (it can also be easily translated into any programming language).

Figure 1 .

 1 Figure 1. Sketch of physical and transformed domains.

F 2 n+1 -F 2 n<Figure 2 .

 222 Figure 2. Solitary wave profile (a) and convergence of the numerical error during the iterations (b).

Figure 3 .

 3 Figure 3. Speed versus amplitude relation predicted by the present (solid line) and Tanaka's (circles) algorithms.

Figure 4 .

 4 Figure 4. A very high amplitude solitary wave (a d ≈ 0.790).

Table 1 .

 1 Amplitude corresponding to F = 1.274236297977903 (wrong digits are underlined).

	Tanaka (16 digits)	Present (16 digits)	Present (50 digits)
	0.6819448200954614 0.6819656673869021 0.68196566738691572070985777997577037303844037027343

Acknowledgements 7 References 7 Acknowledgements Denys Dutykh acknowledges the support from ERC under the research project ERC-2011-AdG 290562-MULTIWAVE. The authors would like to thank Pavel Holodoborodko for invaluable help with the multi-precision floating point arithmetics.


F Doumenc
email: doumenc@fast.u-psud.fr

E Chénier

B Trouette

T Boeck

C Delcarte

B Guerrier

M Rossi

Free convection in drying binary mixtures: solutal versus thermal instabilities

Free convection occuring during the drying of plane layers of polymer solutions may be due to various mechanisms, based on buoyancy or Marangoni effect, of thermal or solutal origin. This theoretical and numerical work provides all the required tools to analyze thoroughly the problem. In this transient flow, different methods (frozen time, non-normal, nonlinear methods) are proposed to predict critical times for convection onset and threshold values for convection. Nonlinear and non-normal methods give similar results, within the uncertainty inherent to any transient problem. It is shown that, when linear stability analysis indicates the presence of several instability mechanisms, it is necessary to invoke nonlinear arguments to establish the leading mechanism. The proposed methodology is then applied to experimental results from the literature for two polymer solutions (Polyisobutylene/toluene and Polystyrene/toluene).

Introduction

Thermal free convection in a horizontal layer of liquid is still an active field of research since the first experimental studies by H. Bénard a century ago. Driv-ing forces can originate from density or surface tension dependency on temperature, respectively responsible of buoyancy forces or surface tension gradients.

Theoretically, the mechanisms leading to hydrodynamic instabilities have been elucidated first by Rayleigh [START_REF] Rayleigh | On convection currents in a horizontal layer of fluid when the higher temperature is on the under side[END_REF] for the buoyancy-driven flow (Rayleigh-Bénard convection), and by Pearson [START_REF] Pearson | On convection cells induced by surface tension[END_REF] for the surface-tension-driven flow (Bénard-Marangoni convection). Pearson showed that, for the most common pure liquids, the thermal instability threshold was driven by Marangoni effect for layer thicknesses approximately lower than 1 cm, and by buoyancy for higher thicknesses. Another instability mechanism due to surface deformation can also be present: it is known as the long-wave instability (Scriven and Sternling [START_REF] Scriven | On cellular convection driven by surfacetension gradients: effects of mean surface tension and surface viscosity[END_REF], Reichenbach and Linde [START_REF] Reichenbach | Linear perturbation analysis of surfacetension-driven convection at a plane interface (Marangoni instability)[END_REF], Goussis and Kelly [START_REF] Goussis | On the thermocapillary instabilities in a liquid layer heated from below[END_REF]). In liquid mixtures, density and surface tension can also depend on concentration. In that case, the above thermal instability mechanisms have also solutal counterparts. For a review of recent developments of Rayleigh-Bénard-Marangoni free convection see Bodenschatz et al. [START_REF] Bodenschatz | Recent developments in Rayleigh-Bénard convection[END_REF], Schatz and Neitzel [START_REF] Schatz | Experiments on thermocapillary instabilities[END_REF], Colinet et al. [START_REF] Colinet | Nonlinear Dynamics of Surface-Tension-Driven Instabilities[END_REF], Nepomnyashchy et al. [START_REF] Nepomnyashchy | Interfacial Convection in Multilayer Systems[END_REF], Manneville [START_REF] Manneville | Rayleigh-Bénard Convection: Thirty Years of Experimental, Theoretical, and Modeling Work[END_REF].

Free convection induced by evaporation has been studied for a long time, both experimentally and theoretically. In pure liquids or solutions, evaporation at the upper free surface generates a temperature gradient induced by cooling through vaporization latent heat. In the case of solutions, this evaporation also produces a concentration gradient due to difference of components volatility. Both gradients might cause four instability mechanisms at the origin of convective patterns: buoyancy via thermal effects (thermal Rayleigh-Bénard), buoyancy via solutal effects (solutal Rayleigh-Bénard), Marangoni via thermal effects (thermal Bénard-Marangoni) and Marangoni via solutal effects (solutal Bénard-Marangoni). Spangenberg and Rowland [START_REF] Spangenberg | Convective circulation in water induced by evaporative cooling[END_REF] used a schlieren technique to characterize the flow pattern driven by thermal buoyancy at the surface of an evaporating thick layer of pure water (10 cm). Berg et al. [START_REF] Berg | Natural convection in pools of evaporating liquids[END_REF] found a great variety of patterns applying the same technique to several pure substances and solutions of different thicknesses. Certain patterns were induced by buoyancydriven-instability, and others by surface-tension-driven instability. In the latter case, these authors demonstrated the importance of surface contamination, very effective for aqueous systems. Indeed, thermal Marangoni convection can hardly be observed in water, because of the high sensitivity of water surface to contamination (Cammenga et al. [START_REF] Cammenga | On Marangoni convection during the evaporation of water[END_REF]). More recently, Toussaint et al. [START_REF] Toussaint | Experimental characterization of buoyancy-and surface tension-driven convection during the drying of a polymer solution[END_REF] and Bassou and Rharbi [START_REF] Bassou | Role of Bénard-Marangoni instabilities during solvent evaporation in polymer surface corrugations[END_REF] performed experimental studies of free convection occuring during the drying of polymer solutions. Zhang and co-authors used shadowgraphs techniques to characterize the patterns observed during the drying of NaCl/water solutions (Zhang et al. [START_REF] Zhang | Marangoni convection in binary mixtures[END_REF]) and ethanol/water solutions (Zhang et al. [START_REF] Zhang | Novel pattern forming states for Marangoni convection in volatile binary liquids[END_REF]).

In these evaporative convection phenomena, which are generally of transient nature, the determination of conditions corresponding to the onset of free convection motivated a lot of works. Some authors performed a linear stability analysis of a one-dimensional steady basic state in a liquid/gas bilayer, when constant temperatures and/or concentrations were prescribed at the bottom and top boundaries (see for instance Merkt and Bestehorn [START_REF] Merkt | Bénard-Marangoni convection in a strongly evaporating fluid[END_REF], Moussy et al. [START_REF] Moussy | Influence of evaporation on Bénard-Marangoni instability in a liquid-gas bilayer with a deformable interface[END_REF], Machrafi et al. [START_REF] Machrafi | Bénard instabilities in a binary-liquid layer evaporating into an inert gas[END_REF]). When the basic state is unsteady as for most of the evaporation processes, various approaches were implemented. The frozen-time approach has been used by Vidal and Acrivos [START_REF] Vidal | Effect of nonlinear temperature profiles on onset of convection driven by surface tension gradients[END_REF] to determine the time of thermal Bénard-Marangoni convection onset in a shallow layer of propyl alcohol. In this approach, one applies a classical normal mode stability analysis to the unsteady basic temperature profile frozen at each given time t. The amplification method, which fully takes into account the transient nature of the basic state, was used by Foster [START_REF] Foster | Stability of a homogeneous fluid cooled uniformly from above[END_REF][START_REF] Foster | Onset of convection in a layer of fluid cooled from above[END_REF] to determine the onset time of thermal Rayleigh-Bénard convection in a deep layer of water (10 cm). Good agreement was obtained with experiments for an amplification factor around 10-100. More recently, Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF] used a non-normal method, which is also based on amplification ideas, to determine the stability conditions of the thermal Rayleigh-Bénard-Marangoni problem in drying polymer solutions. In the thermal problem, Touazi et al. [START_REF] Touazi | Simulation of transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF] determined the stability threshold using direct nonlinear simulations, and found a good agreement with non-normal results. For the transient solutal problem, nonlinear simulations were also used by Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF] to pinpoint these thresholds in a configuration with a flat interface and variable viscosity. Serpetsi and Yiantsios [START_REF] Serpetsi | Stability characteristics of solutocapillary Marangoni motion in evaporating thin films[END_REF] also studied the stability of the solutal problem, for a constant viscosity but a deformable interface, using the frozen state approach and nonlinear simulations.

For a given problem, such analysis helps to exhibit the mechanism producing the convection onset. This is particularly interesting for binary mixtures, where many such mechanisms can be potentially active. Pearson [START_REF] Pearson | On convection cells induced by surface tension[END_REF] extended his famous stability analysis of a pure liquid, to show that a 0.5 mm thick layer of 5 % of ether in liquid paraffin was unstable for the thermal and solutal surfacetension-driven problems, but stable for buoyancy. For the water/ethanol system, Machrafi et al. [START_REF] Machrafi | Bénard instabilities in a binary-liquid layer evaporating into an inert gas[END_REF] took into account buoyancy and surface tension effects, of thermal and solutal origin, as well as Soret effect in their linear stability analysis.

They concluded that the stability threshold was driven by the solutal Bénard-Marangoni mechanism. In drying polymer solutions, de Gennes [START_REF] De Gennes | Instabilities during the evaporation of a film: Nonglassy polymer + volatile solvent[END_REF] used scaling arguments to estimate the critical thickness of the Bénard-Marangoni instability.

He concluded that the solutal critical thickness was much smaller than the thermal one, so the concentration effects should dominate the thermal effects.

However, drying experiments on polymer solutions performed by Toussaint et al. [START_REF] Toussaint | Experimental characterization of buoyancy-and surface tension-driven convection during the drying of a polymer solution[END_REF] exhibit experimental results in agreement with thermal Rayleigh-Bénard-Marangoni convection. The reasons of this apparent contradiction are derived in this paper.

The present article focuses on the leading mechanisms and critical conditions for a specific type of evaporative convection: the one occuring in a plane layer geometry and induced by evaporation of a volatile solvent in which a non-volatile solute is diluted. Such a case is of particular interest for coating applications. Indeed, drying of such mixtures (polymer solutions, colloidal suspensions,. . . ) leads to solid-like deposit on the substrate at the end of the process. In some cases, this deposit presents undesirable undulations (self-patterning). Although ruptures of a superficial glassy crust (de Gennes [START_REF] De Gennes | Solvent evaporation of spin cast films: crust effects[END_REF]) or buckling of a glassy elastic skin (Huraux et al. [START_REF] Huraux | Wrinkling of a nanometric glassy skin/crust induced by drying in poly(vinyl alcohol) gels[END_REF]) are proposed explanations for these undulations, free convection cells observed during drying are often invoked to understand the non uniform thickness of the final deposit (Sakurai et al. [START_REF] Sakurai | Two-dimensional undulation pattern on free surface of polymer film cast from solution[END_REF], Weh [START_REF] Weh | Self-organized structures at the surface of thin polymer films[END_REF], Sakurai et al. [START_REF] Sakurai | Control of mesh pattern of surface corrugation via rate of solvent evaporation in solution casting of polymer film in the presence of convection[END_REF], Gorand et al. [START_REF] Gorand | Instabilités de plissement lors du séchage de films polymères plans[END_REF], Bassou and Rharbi [START_REF] Bassou | Role of Bénard-Marangoni instabilities during solvent evaporation in polymer surface corrugations[END_REF], Minaȓík et al. [START_REF] Minaȓík | Preliminary investigation of factors determining self-organised structures preparation in polymer layers[END_REF], Uchiyama et al. [START_REF] Uchiyama | Spontaneous formation of micrometer-scaled cell-like patterns on alkoxide-derived silica gels induced by Bénard-Marangoni convections[END_REF][START_REF] Uchiyama | Spontaneous formation of linearly arranged microcraters on sol-gel-derived silicapoly(vinylpyrrolidone) hybrid films induced by Bénard-Marangoni convection[END_REF], Jun and Lee [START_REF] Jun | The effect of Bénard-Marangoni convection on percolation threshold in amorphous polymer-multiwall carbon nanotube composites[END_REF]). For instance, Bassou and Rharbi [START_REF] Bassou | Role of Bénard-Marangoni instabilities during solvent evaporation in polymer surface corrugations[END_REF] provides strong experimental evidences connecting solutal Bénard-Marangoni cells and deposit patterning. Apart from this latter work, the convection mechanism (thermal or solutal, surface-tension-or buoyancy-driven) is generally not specified, or it is not supported by clear-cut arguments.

In our previous studies, the thermal problem (Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF], Touazi et al. [START_REF] Touazi | Simulation of transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF]) and the solutal problem (Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF]) have been studied independantly. The present work aims at discriminating which is the dominant mechanism for a given set of experimental parameters corresponding to configurations above the thermal and solutal thresholds. When the solutal and the thermal problems are both unstable, it is shown that a stability analysis cannot provide a complete answer. It is then necessary to study the nonlinear regime to decide which mechanism produces the highest velocity in the liquid layer. This article is organized as follows. In section 2, we briefly present the experimental configuration which have been considered to study this transient problem.

In section 3, models and dimensionless numbers characterizing the thermal or solutal problems are detailed. Section 4 describes the different methodologies used to derive the conditions leading to convection onset or stability thresholds.

In particular, the time and wavelength corresponding to the onset of convection are determined for given non-dimensional numbers. In section 5, comparisons of these results with experiments are presented. First the stability threshold is obtained for both the thermal and the solutal problems. Thereafter, we focus on cases where both thermal and solutal problems are unstable. Critical times, and wavelengths predicted by the thermal and the solutal models are compared with experimental data. In section 6, we focus on the quasi-steady regime in Bénard-Marangoni convection for which we propose a scaling analysis. This provides simple relations for the orders of magnitude of the thermal and solutal velocities in the nonlinear regime. Finally experimental results for the solutal case are compared with the above scaling.

Experimental Configuration

The typical behavior observed during the drying of a polymer solution is illustrated in figure C.1 for a 8 mm thick plane layer. First, a decrease of the evaporative solvent mass flux per unit area, Q m , is observed (0 ≤ t t 1 = 1500 s). The second stage (t 1 t t 2 = 4000 s) is characterized by a nearly constant evaporative flux (see Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF] for details). In the final regime (t 2 t), the evaporative flux decreases. First it slowly decreases then it rapidly falls down (t 3 = 26 000 s t).

This behavior can be accounted for by using a relationship between Q m and the saturated vapor pressure P V S of the volatile solvent above the interface. Indeed let us restrict ourselves to the common situation in which evaporation is limited by vapor diffusion in an inert gas, usually air at atmospheric pressure (liquids in contact with their pure vapor can lead to very different behaviors, see for example Uguz and Narayanan [START_REF] Uguz | Instability in evaporative binary mixtures. i. the effect of solutal Marangoni convection[END_REF][START_REF] Uguz | Instability in evaporative binary mixtures. ii. the effect of Rayleigh convection[END_REF]). In that case, the phenomenological law (see Guerrier et al. [START_REF] Guerrier | Drying kinetics of polymer films[END_REF])

Q m = h m [(c g S ) interf ace -(c g S ) ∞ ] with (c g S ) interf ace = P V S (T, ϕ s ) M S RT (1)
applies where (c g S ) interf ace (resp. (c g S ) ∞ ) denotes the solvent vapor concentration in air at the interface (resp. far from the interface). (c g S ) ∞ is zero in our experiments and (c g S ) interf ace is given assuming local thermodynamic equilibrium. M S stands for the solvent molar mass, R the ideal gas constant. The mass transfer coefficient h m depends on the air velocity above the solution which is kept constant during all experiments. It is evaluated from measurements obtained from the evaporation kinetics of a pure toluene layer (Toussaint et al. [START_REF] Toussaint | Experimental characterization of buoyancy-and surface tension-driven convection during the drying of a polymer solution[END_REF], Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF]).

The dependency of P V S on both temperature T and solvent concentration at the interface ϕ s explains the succession of several regimes for evaporation onset, we use two separate models: a thermal model and a solutal one. We therefore ignore thermal and solutal coupling effects. This assumption is consistent with our objective, which is the determination of the leading mechanism at short time, just after convection onset. A complete thermo-solutal model would be desirable for a more detailed analysis (like for instance prediction of 3D patterns), or to investigate longer times, when the thermal problem depends on concentration through the viscosity. In addition, Soret effect is supposed to be negligible (this assumption will be justified a posteriori by scaling arguments) and we disregard the long-wave instability. This latter mechanism may be pertinent in some experiments (see the numerical simulations by Yiantsios and Higgins [START_REF] Yiantsios | Marangoni flows during drying of colloidal films[END_REF][START_REF] Yiantsios | A mechanism of Marangoni instability in evaporating thin liquid films due to soluble surfactant[END_REF] for the thermal and solutal problems) but it has not been observed in the experiments presented in section 2.

flux Q m . The decrease of Q m for 0 ≤ t
Both thermal and solutal models share common assumptions: (a) The liquid layer is located between the bottom at z = 0 and the upper free surface which remains flat at z = d(t). The layer width W is assumed large with respect to its height d(t). (b) Bottom and side walls are adiabatic/impermeable, with a no slip condition for velocity. (c) The solution is Newtonian of viscosity µ.

Actually non-Newtonian effects are expected in polymer solutions, but only for high enough shear rates, larger than the inverse of the viscoelastic relaxation time of the solution. For test case 2, the shear rate was estimated to be of order of 1 s -1 . For PIB/toluene at the initial concentration of test case 2 and the same PIB molar mass, no viscoelastic effect were detected in rheometer measurements up to 50 s -1 shear rate (Gorand et al. [START_REF] Gorand | Instabilités de plissement lors du séchage de films polymères plans[END_REF]). So viscoelastic effects are unlikely to occur in the beginning of the drying, even if they cannot be excluded at later stage. (d) Liquid density ρ is assumed constant, except in the buoyancy term. (e) Local thermodynamic equilibrium is assumed at the upper evaporative surface. (f) Both models are based on a one-layer approximation in which the heat (mass) transfer in the vapor phase is given by heat (mass) transfer coefficients. (g) For the sake of simplicity, only two-dimensional flows are studied: for the linear problem, the two-dimensional assumption is not a constraint. In addition, 2D and 3D computations gave identical answers on the thermal nonlinear problem (Trouette et al. [START_REF] Trouette | Numerical study of convection induced by evaporation in cylindrical geometry[END_REF]), when evaluating the orders of magnitude attached to the different processes as performed here.

Since we focus on the beginning of the drying (a few minutes), the liquid height and viscosity variations can be neglected. This assumption has been used in the thermal model. In the solutal model, which has been developped in a previous work to simulate solutal convection over longer times (Trouette et al. [48]), liquid height and viscosity variations are taken into account. For a detailed description of the validity domain of the thermal model, see Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF], Touazi et al. [START_REF] Touazi | Simulation of transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF]. The assumptions and the validity domain of the solutal model can be found in Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF].

Thermal model

The velocity field, v = v x e x + v z e z , is governed by the Navier-Stokes equations in the context of Boussinesq approximation: density ρ is taken to be the density at T = T ∞ except for the buoyancy term where one sets

ρ(T ) = ρ(T ∞ ) [1 -α th (T -T ∞ )] (2)
with α th the thermal expansion coefficient. Surface tension σ(T ) is a linearly decreasing function of temperature T .

σ(T ) = σ(T ∞ ) -γ th (T -T ∞ ) (3)
where γ th > 0. Finally, the fluid is characterized by a thermal diffusivity κ. At the upper surface z = d, the balance of tangential forces reads

µ ∂v x ∂z = -γ th ∂T ∂x at z = d, (4)
and the conservation of energy flux is ensured through the phenomenological equation (for a derivation see Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF])

-λ ∂T ∂z + H th (T ∞ -T (z = d)) = LQ m (T ∞ ), H th = h th + L ∂Q m ∂T (T ∞ ) (5)
where λ denotes the thermal conductivity of the liquid, h th the heat transfer coefficient between the gas and the liquid.

The governing equations are put in a dimensionless form using the constant layer thickness d, the thermal velocity V = κ/d and thermal diffusion time d/V = d 2 /κ for length, velocity and time scales. Pressure p is scaled by κµ/d 2 . Finally one introduces the dimensionless quantity θ(x, z, t)

= (T (x, z, t) -T ∞ )/∆T st where ∆T st ≡ LQ m H th . (6)
denotes the steady temperature difference in the final steady state regime (t 1 t). The superscript * is used for dimensionless variables of the thermal model.

In Cartesian coordinate system, the equations for dimensionless velocity v * , dimensionless deviation from the hydrostatic pressure p * and dimensionless temperature θ read

∇ • v * = 0, (7) 1
P r ∂v * ∂t * + (v * • ∇)v * = -∇p * + ∇ 2 v * + Ra th θ e z , (8)
∂θ ∂t * + (v * • ∇)θ = ∇ 2 θ, (9) ∂v
* x ∂z * + M a th ∂θ ∂x * = 0 at z * = 1. ( 10
)
∂θ ∂z * + Bi θ + Bi = 0 at z * = 1, (11)
∂θ ∂z * = 0 at z * = 0, ∂θ ∂x * = 0 at x * = 0, A, (12)
v * z = 0 at z * = 1, v * x = v * z = 0 at z * = 0 or x * = 0, A. (13)
This system depends on five dimensionless parameters: Rayleigh, Marangoni, Biot, Prandtl numbers and an aspect ratio

Ra th ≡ α th gρd 3 ∆T st µκ , M a th ≡ γ th d∆T st µκ , Bi ≡ H th d λ , P r ≡ µ ρκ , A ≡ W d . ( 14
)
The analysis is restricted to Prandtl number P r ≥ 1. This assumption is clearly valid for most of the liquids, if one excepts liquid metals. It means that the thermal diffusion time scale is always larger than the viscous diffusion time scale. In all the simulations, the aspect ratio is large (A = 20 for the thermal problem). It was checked that, when large enough, A does not affect our results.

Solutal model

In the solutal model, the temperature is assumed constant (T ≃ T ∞ ) and the fluid is a binary solution with a constant diffusion coefficient D. The system is characterized by the solvent volume fraction ϕ s (x, z, t) or the polymer volume fraction ϕ p (x, z, t) = 1 -ϕ s (x, z, t). During drying, these concentrations are inhomogeneous and unsteady up to the very end when the film is totally dry

(ϕ s = 0).
The velocity field is governed by the Navier-Stokes equations. The difference between pure polymer and pure solvent densities being around 6% for PIB/toluene and 20% for PS/toluene, density ρ is assumed constant and equal to ρ i the density at the initial solvent volume fraction ϕ si , except for the buoyancy term where one sets

ρ = ρ 0 (1 + α sol (1 -ϕ s )) ( 15
)
where ρ 0 denotes the density of the pure solvent and α sol > 0 the solutal expansion coefficient. The interface possesses a surface tension σ which is a linear

function of ϕ s σ = σ 0 + γ sol (1 -ϕ s ), (16)
with σ 0 the surface tension of the pure solvent and γ sol > 0.

(µ) = a 0 + a 1 Y + a 2 Y 2 + a 3 Y 3 + a 4 Y 4 (17)
with Y = log 10 (ϕ p ) and a i are fixed coefficients (for details see Appendix A).

We assume that the evaporative flux remains constant. As mentioned in section 2, this is only approximately valid during the period t t 1 (typically 30 percent variation due to evaporation induced cooling) and experimentally verified during the period t 1 t t 2 . It is thus possible to simulate the drying process up to time t 2 when the solvent volume fraction reaches about 0.4 (Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF]) after which the flux is no more steady. The evaporative flux being constant and the surface remaining flat, the layer thickness d(t) = d(0)-v ev t decreases at a constant velocity v ev . At the upper surface, the balance of the tangential stresses imposes

µ ∂v x ∂z = -γ sol ∂ϕ s ∂x at z = d(t) (18)
The solvent and polymer mass conservation leads to two additional boundary conditions (see Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF] for a derivation)

v z = 0 at z = d(t), (19)
-D ∂ϕ s ∂z = v ev (1 -ϕ s ) at z = d(t), (20)
Dimensionless equations are obtained by scaling coordinates x and z, velocity v, dynamic pressure p and time t respectively by d i , D/d i , µ i D/d 2 i and d 2 i /D where d i is the initial thickness and µ i the initial viscosity. Superscripts + are used for dimensionless variables of the solutal model. In addition to v + , p + , one introduces, based on the solvent volume fraction ϕ, the quantity

φ s = ϕ s -ϕ si ∆ϕ
where ϕ si is the initial solvent volume fraction and ∆ϕ is based on the concentration gradient near the interface and the layer thickness

∆ϕ ≡ -d i ∂ϕ s (z = d i , t = 0) ∂z int = P e int (1 -ϕ si ) ( 21
)
where ϕ si is the initial solvent volume fraction and P e int ≡ v ev d i /D is a Péclet number based on the interface velocity (Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF]).

This scaling is customary when the evaporation flux is imposed at the boundary. The dimensionless equations for v + , p + and φ s then read

∇ • v + = 0, (22) 1 Sc ∂v
+ ∂t + + (v + • ∇)v + = -∇p + + ∇ • µ µ i ∇v + + ∇ T v + + Ra sol φ s e z , (23)
∂φ s ∂t + + (v + • ∇)φ s = ∇ 2 φ s , (24)
∂v

+ x ∂z + + µ i µ M a sol ∂φ s ∂x + = 0 at z + = 1 -P e int t + . ( 25
)
∂φ s ∂z + + 1 -P e int φ s = 0 at z + = 1 -P e int t + , (26)
∂φ s ∂z + = 0 at z + = 0,

∂φ s ∂x + = 0 at x + = 0, A, (27)
v + z = 0, at z + = 1 -P e int t + v + x = v + z = 0 at z + = 0 or x + = 0, A. (28)
This problem depends on the initial solvent volume fraction ϕ si through the concentration-dependent viscosity µ/µ i (see Eq. ( 17)), on the Péclet number P e int , as well as on the Rayleigh, Marangoni, Schmidt dimensionless parameters and on the initial aspect ratio A

Ra sol ≡ α ′ sol gρd 3 i ∆ϕ µ i D , M a sol ≡ γ sol d i ∆ϕ µ i D , P e int ≡ v ev d i D , Sc ≡ µ i ρ i D , A ≡ W d i (29) with α ′ sol = α sol /[1 + α sol (1 -ϕ si )]
. In all simulations, the aspect ratio is large (5 ≤ A ≤ 30 for the solutal problem). It was checked that, when large enough, A does not affect our results.

Convection Onset and Marginal Stability Determination Methods

The configurations studied are known to display a transition from a purely diffusive state to convective patterns. Since the undisturbed pure diffusive state is unsteady, determination of convection onset amounts at defining a time for convection onset and a corresponding wavenumber. In the subsection below, we present linear and nonlinear methods to determine when does the convection start for a given set of control parameters ( 14) or ( 29), and at which wavelength. This leads also to the determination of the marginal stability curve for this unsteady problem.

Time and wavenumber corresponding to Convection Onset: Nonlinear Approach

This first method necessitates nonlinear two-dimensional simulations to be performed. For both thermal and solutal problems, they are based on finite volume schemes (see Touazi et al. [START_REF] Touazi | Simulation of transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF] and Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF] for details). In order to take into consideration the moving upper surface in the solutal problem, this numerical method is based on moving grids (Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF]). In both models, the initial velocity field is set to zero. For the thermal problem, the initial dimensionless temperature field θ(x * , z * , t * = 0) is a random perturbation added at each discretized spatial location. This random perturbation is of zero mean and uniformly distributed between -r * /2 and r * /2 [START_REF] Trouette | Numerical study of convection induced by evaporation in cylindrical geometry[END_REF][START_REF] Chénier | Sensitivity of diffusive-convective transition to the initial conditions in a transient Bénard-Marangoni problem[END_REF]. For the solutal problem, the initial dimensionless solvent volume fraction φ s (x + , z + , t + = 0) is also a similar random field of zero mean and uniformly distributed between -r + /2 and r + /2.

Let us denote by < q(z 0 , t) > the mean of quantity q(x, z 0 , t) spatially averaged over the horizontal plane z = z 0 . For the thermal problem, these simulations provide the difference < T (z = 0, t) > -< T (where superscript "nl" stands for "nonlinear") can be exhibited that corresponds to the moment where the norm of the dimensionless velocity becomes greater than a threshold. This threshold which characterizes the deviation from pure diffusive state is here defined such that the dimensional velocity becomes greater than the heat diffusion velocity κ/d. In dimensionless units, it corresponds to v * 2 > 1 (Touazi et al. [START_REF] Touazi | Simulation of transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF]). Since it depends on the initial perturbation amplitude and on the threshold, there is some blurredness in the evaluation of time t * (nl) on . On the contrary, the saturated nonlinear regime reached after t * (nl) on is independent on the initial state (see figure C.4). This subsequent regime is called quasi-steady in the following since the temperature difference varies much less during this period than during the period close to time of convection onset. Note that, at time t * (nl)

on , one may obtain the average cell-to-cell distance providing a wavenumber k * (nl) on . The same qualitative behavior is observed in the solutal problem for < ϕ s (z = 0, t) > -< ϕ s (z = d(t), t) > (see figure C.5 for the equivalent experimental configuration). The method is then similar, the threshold corresponding to

v + 2 > 1.

Time and wavenumber corresponding to Convection Onset: Linear Approaches

Other methods which are used to determine convection onset are quite different being based on linear arguments. They are here presented only for the thermal problem with infinite aspect ratio. The linear theory accounts for the dynamics of infinitesimal perturbations v * p (x * , z * , t * ), θ p (x * , z * , t * ), p * p (x * , z * , t * ) superposed near the purely diffusive basic state θ BS (z * , t * ) (which is simulated setting Ra th = M a th = 0). They are assumed of the form

(v * p , θ p , p * p ) = (v x (z * , t * ), vz (z * , t * ), θ(z * , t * ), p(z * , t * )) exp(ik * x * ) ( 30
)
where k * denotes the nondimensional wavenumber along the x * direction. The complex amplitude of the perturbations are governed by the linear system

ik * vx + ∂v z ∂z * = 0, (31) 1 P r
∂ ∂t * vx + ik * p - ∂ 2 ∂z * 2 -k * 2 vx = 0, ( 32
) 1 P r ∂ ∂t * vz + ∂ p ∂z * - ∂ 2 ∂z * 2 -k * 2 vz -Ra th θ = 0, (33)
∂ ∂t * θ + vz ∂θ BS ∂z * - ∂ 2 ∂z * 2 -k * 2 θ = 0, ( 34
) vz = 0, ∂ z * vx + M a th ik * θ = 0, ∂ z * θ + Bi θ = 0, at z * = 1, ( 35
) vx = vz = 0, ∂ θ ∂z * = 0 at z * = 0. ( 36
)
In this linear framework, various alternatives are possible to evaluate the time of convection onset. First, one may use a frozen-time approximation which considers at each time t * the unsteady temperature basic profile θ BS (z * , t * ) and applies the classical normal mode stability method to this frozen state θ BS (z * , t * ). The first time for which the stability spectrum of all wavenumbers contains an eigenvalue with a zero real part, is defined as the onset time t * (f r) on .

This naturally leads to a critical wavenumber k * (f r) on as well.

Another linear approach explicitly takes into account the unsteady character of basic profile θ BS (z * , t * ). It is based on the non-normal method (for details see Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF]). For the non-autonomous system (31)- [START_REF] Roots | Cooperative diffusion in semidilute polystyrene solutions at good and solvent conditions[END_REF], the standard notion of growth rate based on eigenvalues is no more valid. One should resort to amplification gains as a function of time to characterize the flow stability. Given an initial disturbance profile, a time t * , and a wavenumber k * , one defines two different norms E V and E T for the perturbation amplitudes. The first norm is based on the kinetic energy of perturbations

E V (t * , k * ) ≡ (v x (z * , t * , k * )v # x (z * , t * , k * ) + vz (z * , t * , k * )v # z (z * , t * , k * ))dz * (37)
where superscript # denotes complex conjugation. The integration is performed over the entire layer height and perturbations are obtained after integrating the above linear non-autonomous system (31)-( 36) over the time period [0, t * ]. The second norm is based not on the velocity field but on the temperature field:

E T (t * , k * ) ≡ θ(z * , t * , k * ) θ# (z * , t * , k * )dz * (38)
The amplification gain at time t * and for each wavenumber k * is then evaluated by computing factor

E V (t * , k * )/E V (0, k * ) or E T (t * , k * )/E T (0, k * ). Given
a norm, the non-modal analysis solves an optimization problem based on an iteration procedure and the introduction of an adjoint problem (Farrell and Ioannou [START_REF] Farrell | Generalized stability theory part I: autonomous operators[END_REF], Andersson et al. [START_REF] Andersson | Optimal disturbances and bypass transition in boundary layers[END_REF], Luchini [START_REF] Luchini | Reynolds-number-independent instability of the boundary layer over a flat surface: optimal perturbations[END_REF], Schmid and Henningson [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF]) which allows: (a) to determine, for given wavenumber k * and time t * , the maximum en-

ergy amplification Ĝ(t * , k * ) ≡ M ax[E(t * , k * )/E(0, k * )
] over all possible initial perturbations profiles; and (b) to exhibit the optimal perturbation mode defined by its initial z * -profile, which actually reaches this upper bound (Schmid and Henningson [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF]).

Using the value Ĝ(t * ; k * ), it is feasible to extend to unsteady flows the usual concepts of classical stability analysis. For instance, Ĝ(t * ; k * ) can be maximized over wavenumber k * providing a maximum amplification G max (t * ). This value is effectively reached for a specific wavenumber k * max (t * ) and for a specific initial perturbation structure in z * . These latter two quantities play the role of the most amplified wavenumber and of the most amplified mode for the standard analysis but at a given time t * . Let us now consider a threshold G thres . When the system is capable to reach such a value G thres , the basic flow is considered as unstable. The first time when G max (t Based on this remark, in most of the cases considered in this work, we will mainly use the nonlinear approach to determine the conditions of convection onset. Moreover direct numerical simulations also give a more complete description by providing velocity fields in the nonlinear regime and time evolution of the wavelengths.

Critical Control Parameters

For some sets of control parameters (M a th , Ra th , Bi, P r), the critical times t * (nn) on , t * (f r) on or t * (nl) on cannot be defined. For instance, this happens in the nonnormal mode approach when G max (t * ) always remains below G thres . In that case, convection is predicted not to occur and the system to be stable through the instability mechanism considered. Critical Rayleigh or Marangoni numbers can be then defined as the smallest Rayleigh or Marangoni numbers for which times of convection onset appear. In the control parameters space, a marginal stability curve such as Ra c th (Bi, M a th , P r) can also be computed, each point of the curve being associated to critical optimal time t * c = t * on (Ra c th , M a th , Bi, P r) and a critical wavenumber k * c = k * on (Ra c th , M a th , Bi, P r). The results of marginal stability curve has been extensively analyzed in previous studies. See Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF], Touazi et al. [START_REF] Touazi | Simulation of transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF] for the critical conditions of the thermal problem, and Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF] for the solutal problem. Here it is introduced because of the comparison with experiments performed in section 5.1.

Comparisons of numerical predictions with results of PIB/toluene

experiments. The frontier delimiting the observation of convection in PIB/toluene experiments is in agreement with the region of transition from stable to unstable configurations for the thermal Rayleigh-Bénard-Marangoni problem (two continuous lines in figure C.8). Indeed no convection was observed for configurations in which solutal convection is expected (empty squares in figure C.8). This seems paradoxical since, based only on marginal curves, the solutal Bénard-Marangoni mechanism appears to be the most unstable (an exception would be the lowest viscosities corresponding to solutions very close to pure solvent, where solutal convection is not defined). In addition, the predominance of solutal with respect to thermal instability given by theoretical model agrees with previous works (de Gennes [START_REF] De Gennes | Instabilities during the evaporation of a film: Nonglassy polymer + volatile solvent[END_REF], Machrafi et al. [START_REF] Machrafi | Bénard instabilities in a binary-liquid layer evaporating into an inert gas[END_REF]). This paradox is explained in section 5.3.

Comparisons of critical control parameters

In the subsections below, we focus on configurations above the thermal and solutal stability thresholds. to determine in a clear cut manner the nature of the observed flow.

Comparisons of

Comparisons of wavelengths and velocities when thermal and solutal problems are both unstable

The wavelengths obtained by numerical simulation or experimental visualization are compared in figure C.9. For both test cases, the observed wavelengths are in very good agreement with the thermal model results while the solutal wavelengths are significantly lower. This is true near the onset but remains so in the later saturated regime.

The preponderance of the thermal mechanism for the two configurations is confirmed by the quantification of velocities (see figure C.10). Velocities induced by thermal effects predominate over the one of solutal origin: thermal velocities are an order of magnitude larger than solutal ones for test case 1, half a decade for test case 2. When the thermal convection is not active, velocities induced by solutal gradients exist but are too small to be detected by visualization used in the mentioned experiments (see section 6.2 for a more general argument). This explains the paradox mentioned in subsection 5.1. As a conclusion, the analysis of these two test cases clearly highlights the need for a complete characterization of convective regimes to determine the dominant mechanism. Beyond threshold estimation, the saturating amplitude of the instability is also an important parameter.

6. Quasi-steady regime in the Bénard Marangoni convection.

Previously, we have seen that estimating the velocity field in the quasisteady regime may be a way to determine which is the instability mechanism that predominates when both mechanisms are simultaneously unstable. Here we focus on the quasi-steady regime where it is possible to establish scaling laws for velocity field and thus to compare thermal to solutal contributions.

This generalization is performed for the Bénard-Marangoni problem taking into account both thermal and solutal effects but neglecting buoyancy (Ra th = Ra sol = 0). This is a relevant assumption for thin films encountered in coating applications.

Scaling laws and comparison with numerical simulations

Based on some simplifying assumptions, scaling laws are obtained by "solving" equations of thermal and solutal models in terms of order of magnitude.

All the equations below and in Appendix B must be hence understood as orders of magnitude only since these scalings do not give prefactors. Note that this approach is similar to the analysis of a transient free convection problem described in Bejan [START_REF] Bejan | Convection heat transfer[END_REF], and has already been applied to the thermal Bénard-Marangoni problem in an evaporative liquid (Touazi et al. [START_REF] Touazi | Simulation of transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF]). In the following, the notation "x y" (resp. "x y") stands for "x ∼ y or x ≪ y" (resp. "x ∼ y or

x ≫ y"). The main assumptions are as follows:

H1 Below the free surface, we assume the existence of a hydrodynamic, a thermal as well as a solutal boundary layer of respective thicknesses δ H , δ T and δ S .

H2 Time derivative terms are neglected (quasi-steady regime).

H3 Temperature (resp. solvent concentration) variations across a convective cell in the vertical and horizontal directions are of same order of magnitude. They are denoted by ∆θ (resp. ∆φ s ).

H4 The wavelength of convective structures scales with the layer thickness.

Hence the order of magnitude of the dimensionless characteristic length in the horizontal direction is one.

H5 The analysis is restricted to fluids characterized by a Prandtl number P r 1 (resp. Schmid number Sc 1).

H6 Even for solutal model, viscosity and layer thickness are assumed constant and equal to their initial values. This is valid since we focus on the beginning of drying.

Scaling analysis of the thermal problem and its validation by comparison with numerical simulations have been presented in a previous paper (see Touazi et al. [45]). Scaling laws of the solutal configuration are derived in Appendix B. As a final comment, this analysis assumes a priori the existence of convection patterns and of quasi-steady regime. Frontiers obtained by this analysis must be compatible with this assumption. Region A (no convection) cannot be thus determined using simply this approach and necessitates the use of a stability analysis. This has been done for the thermal problem for 10 -3 < Bi < 10 3

(Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF]). For the solutal model, numerical simulations for P e int 1

gives M a sol ≃ 100 (see Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF]) which is consistent with the scaling

M a sol ∼ 1 for the A-B frontier (see figure C.11). Thermal problem Solutal problem domain δ H δ T ∆θ scaling laws δ H δ S scaling laws A no convection B ∼ 1 ≪ 1 ≪ 1 δ T ∼ (M a th .Bi) -1/3 ∼ 1 ≪ 1/P e int δ S ∼ M a -1/3 sol ∆θ ∼ M a -1/3 th .Bi 2/3 ∆φs ∼ M a -1/3 sol v * x ∼ (M a th .Bi) 2/3 v + x ∼ M a 2/3 sol C ∼ 1 ≪ 1 ∼ 1 δ T ∼ M a -1/2 th ∼ 1 ∼ 1/P e int ∆φs ∼ M a -1 sol .P e 2 int v * x ∼ M a th v + x ∼ P e 2 int D ≪ 1 ≪ 1 ∼ 1 δ T ∼ M a -1/3 th .P r -1/6 ≪ 1 ∼ 1/P e int ∆φs ∼ M a -1 sol .Sc -1/2 .P e 3 int δ H ∼ M a -1/3 th .P r 1/3 δ H ∼ Sc 1/2 .P e -1 int v * x ∼ M a 2/3 th .P r 1/3 v + x ∼ P e 2 int E ≪ 1 ≪ 1 ≪ 1 δ T ∼ (M a th .Bi) -1/4 .P r -1/8 ≪ 1 ≪ 1/P e int δ S ∼ M a -1/4 sol .Sc -1/8 δ H ∼ (M a th .Bi) -1/4 .P r 3/8 δ H ∼ M a -1/4 sol .Sc 3/8 ∆θ ∼ M a -1/4 th .Bi 3/4 .P r -1/8 ∆φs ∼ M a -1/4 sol .Sc -1/8 v * x ∼ (M a th .Bi) 1/2 .P r 1/4 v + x ∼ M a 1/2 sol .Sc 1/4
Table 3: Quasi-steady regime in the BM configuration -Scaling laws.

Comparison of thermal and solutal velocities

From comparison of scaling analysis with numerical simulations, the domain explored in typical drying experiments like test case 1 is shown to correspond to B domains. This is true in thermal and solutal problems. Let us use the scalings of table 3 to compare thermal and solutal velocities.

domain domain boundary equation boundary equation thermal problem solutal problem

A B Bi × M a th ∼ 1 M a sol ∼ 1 A C Bi ∼ M a th unknown B C Bi 2 × M a -1 th ∼ 1 M a sol ∼ P e 3 int B E Bi 2/3 × M a 2/3 th ∼ P r M a sol ∼ Sc 3/2 C D M a th ∼ P r P e int ∼ Sc 1/2 D E Bi 6 × M a -2 th ∼ P r M a sol ∼ P e 4 int .Sc -1/2
v * x v + x ∼ M a th .Bi M a sol 2/3 , (39)
or in dimensional quantities:

v th v sol ∼ Le -1/3 γ th Lc -1 γ sol 2/3 ϕ -2/3 pi Q m ρv ev 2/3 (40)
with Le ≡ κ/D the Lewis number, c the heat capacity (c = λ/(ρκ)) and Q m the evaporation flux at the beginning of the drying. The two first terms depend on the system properties only, the two last ones on experimental conditions.

The Lewis number compares the solutal diffusion characteristic time to the thermal one. It is of order of 10 3 for both PIB/toluene and PS/toluene solutions.

The second dimensionless number (between square brackets) compares the two driving mechanisms (≃ 5.2 for PIB/toluene, 3.5 for PS/toluene). Notice that the difference between the evaporation flux at the beginning of the drying Q m and during the regime at constant evaporation flux (ρv ev ) is due to evaporation cooling effect. Nevertheless, they are very close (the ratio is about 1.3 for the test case 2, cf. figure C.1), so that the last term of Eq. 40 will be equal to one as an order of magnitude.

Remarkably enough, the ratio of the thermal and solutal velocities does not depend on the sample thickness (see Eq. 40). For a quantitative comparison, note that the prefactors obtained through numerical simulations are close to 0.2 for thermal (see Touazi et al. [START_REF] Touazi | Simulation of transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF]) and solutal problem (figure C.12). Hence, for both solutions (PIB or PS in toluene), the ratio of the thermal and solutal velocities thus reads

v th v sol ≃ 0.3 ϕ -2/3 pi (41)
The above result has been derived for pure Marangoni configurations, i.e. neglecting buoyancy, and assuming that both thermal and solutal problems are unstable. In such a configuration, it leads to the following conclusion: when the thermal and solutal mechanisms are both involved, comparison of stability thresholds or critical times is not sufficient. It is also necessary to consider the velocities induced by the two phenomena. For dilute solutions (e.g. ϕ pi = 0.006 as in experimental test case 1), the thermal velocity is one order of magnitude higher than the solutal one, which is consistent with the result presented in figure C.10. The two velocities would reach the same order of magnitude for ϕ pi ≃ 0.2. Higher polymer volume fraction would induce very high viscosity (see figure C.3), so we can conclude that, for the two polymer solutions under investigation, if thermal convection is present, it could hardly be dominated by solutal convection.

Comparison with PS/toluene experiments

The set of experiments performed by Bassou and Rharbi [START_REF] Bassou | Role of Bénard-Marangoni instabilities during solvent evaporation in polymer surface corrugations[END_REF] concerns thin PS/toluene films with initial thickness d i = 1.4 mm. Experimental configurations and corresponding dimensionless numbers are detailed in table 1 (cases 3 to 5). Velocity measurements have been performed using the particle tracking method, with an optical microscope connected to a CCD camera. This method allows the measurement of small velocities (a few µm.s -1 , cf. Bassou and Rharbi [START_REF] Bassou | Role of Bénard-Marangoni instabilities during solvent evaporation in polymer surface corrugations[END_REF] for experimental details). Using the critical Rayleigh and Marangoni numbers from Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF] for the thermal problem, and from Trouette et al. or by uncertainties over system properties (viscosity, diffusion coefficient,...) or else by the choice of velocity L2 norm to express surface velocity from numerical simulations. The reason of such a discrepancy should hence be attributed to some model assumptions. Soret effect is neglected, but this assumption seems correct: using the scaling laws from table 3 (see Appendix C), it is shown that Soret effect generates a polymer flux which is about two orders of magnitude lower than the flux driven by the concentration gradient due to evaporation.

Turning to the modelling of the free surface, solution surface tension is expressed as a function of bulk concentration near the interface (see Eq. 16).

Surface tension actually depends on the surface concentration. In order to get a univoque relationship between the bulk concentration near interface and the surface concentration itself, local thermodynamic equilibrium between the bulk and the interface must be satisfied. Otherwise (out of equilibrium surface), one should consider the dynamic of the solvent transport at the free surface, as it has already been done with soluble surfactants (see for example Yiantsios and Higgins [START_REF] Yiantsios | A mechanism of Marangoni instability in evaporating thin liquid films due to soluble surfactant[END_REF], Shklyaev, S. and Nepomnyashchy, A.A. [START_REF] Shklyaev | Marangoni instability of a heated liquid layer in the presence of a soluble surfactant[END_REF]). This would necessitate a more elaborate model, but also experiments performed with well characterized polymer solutions for its validation.

Conclusion

In this study, we present a detailed analysis of free convection occuring during the drying of plane layers of polymer solutions. Four mechanisms, based on buoyancy or Marangoni effect, of thermal or solutal origin, are taken into account. We show that when a configuration is unstable for several mechanisms, the comparison of stability thresholds may be insufficient to predict which mechanism dominates: a nonlinear approach is then needed. This theoretical and numerical work combined with previous studies (Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF], Touazi et al. [START_REF] Touazi | Simulation of transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF], Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF]) provides all the required tools to analyze thoroughly the problem. [START_REF] Toussaint | Experimental characterization of buoyancy-and surface tension-driven convection during the drying of a polymer solution[END_REF]). The coupling between the two problems should also be needed to perform a more detailed analysis, like for instance 3D pattern evolution, especially when the thermal and solutal mechanisms lead to comparable velocities. Another problem requiring coupled models concerns the case of thermal and solutal mechanisms acting in opposite way, one stabilizing the system.

Appendix B. Quasi-steady regime in the Bénard-Marangoni solutal convection: scaling laws

In this appendix, scaling laws for the Bénard-Marangoni solutal convection are derived in the quasi-steady regime. Since we consider the beginning of the quasisteady regime, we assume d ∼ d i , µ ∼ µ i and ϕ s ∼ ϕ si in the bulk (i.e. below the solutal boundary layer). From these assumptions and after the removal of time partial derivatives (quasi-steady regime), governing equations ( 22), ( 23) and ( 24) read (for clarity sake, supercripts + are omitted on dimensionless variables):

∂ x v x + ∂ z v z = 0, (B.1) 1 Sc (v x ∂ x v x + v z ∂ z v x ) = -∂ x p + ∂ 2 x v x + ∂ 2 z v x , (B.2) 1 Sc (v x ∂ x v z + v z ∂ z v z ) = -∂ z p + ∂ 2 x v z + ∂ 2 z v z , (B.3) v x ∂ x φ s + v z ∂ z φ s = ∂ 2 x φ s + ∂ 2 z φ s . (B.4)
If one eliminates the pressure gradient term from equations (B.2) (B.3), an equation for vorticity ω is derived

v x ∂ x ω + v z ∂ z ω = Sc(∂ 2 x + ∂ 2 z )ω with ω ≡ ∂ z v x -∂ x v z (B.5)
which expresses a balance between inertia (left-hand side) and friction (righthand side) in the hydrodynamic boundary layer. Finally boundary conditions become at the bottom

v x = v z = 0, ∂ z φ s = 0 at z = 0 (B.6)
and at the upper surface

∂ z φ s + 1 -P e int φ s = 0 (B.7) ∂ z v x + M a sol ∂ x φ s = 0. (B.8) v z = 0. (B.
≡ φ s (x, z = 0) -φ s (x, z = 1)
the variation of solvent volume fraction and velocity can be found:

v x ∼ δ H M a sol ∆φ s . (B.12)
The transport equation (B.4) expresses the balance between the solvent supplied by convection along the x-axis, the solvent supplied by convection along the zaxis, and the diffusion in the liquid. These three terms are respectively of order

Appendix C. Soret effect

Our objective is to show that Soret effect is negligible compared to the diffusion flux driven by the concentration gradient. We consider the configuration of Soret effect is negligible if the condition j s j f ≪ 1 is satisfied. The order of magnitude of j f is readily obtained using the scaling laws of B region in table 3:

j f ∼ Dρ ∆ϕ d i ∆φ s δ s ∼ Dρ ∆ϕ d i (C.2)
A higher bound of j s can be derived from the temperature gradient of the diffusive thermal state (indeed, solutal convection can only decrease the temperature gradient by mixing). For the configurations under consideration, Bi 1 and the dimensionless temperature gradient reads ∂θ ∂z * ∼ Bi (see the study of the thermal basic state by Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF]). We get:

j s DρS T ϕ p ϕ s ∆T st d i Bi (C.3)
The desired relation follows:

j s j f Le -1 (S T Lc -1 ) ( Q m ρv ev ) ϕ si (C.4)
with Le ≡ κ/D the Lewis number, L the latent heat of vaporization and c the heat capacity. From Appendix A we get Le ∼ 10 3 and (S T Lc -1 ) ∼ 10. The two last terms are of order of 1, so the final result is:

j s j f 10 -2 ≪ 1. (C.5)
This result clearly supports the assumption of negligible Soret effect. [17], Mark [START_REF] Mark | Polymer Data Handbook[END_REF] for PIB/toluene, Bassou and Rharbi [START_REF] Bassou | Role of Bénard-Marangoni instabilities during solvent evaporation in polymer surface corrugations[END_REF] for PS/toluene).

  (z = d, t) > between mean temperatures at free surface and bottom as a function of time. The typical evolution is displayed in figure C.4. A deviation from the pure diffusive state (thick line which is obtained by setting Ra th = M a th = 0) is taken as a signature of free convection onset. Quantitatively a time t * (nl) on

  Experimental data and previous numerical results for thermal and solutal models which provided stability thresholds are gathered in figure C.8. The diagram is drawn in the initial thickness/viscosity plane for the PIB/toluene solution (the initial viscosity is related to the initial polymer volume fraction through equation 17). The theoretical thresholds separate the diagram in three regions (a) a stable region for both solutal and thermal mechanisms; (b) a region where solutal Bénard-Marangoni mechanism generates an instability; (c) a region where both solutal and thermal mechanisms generate an instability.

[ 48 ]

 48 for the solutal problem, it is possible to show that these configurations are all stable for the thermal and solutal Rayleigh-Bénard problem as well as the thermal Bénard-Marangoni problem, and unstable for the solutal Bénard-Marangoni problem. This confirms that the flow observed in experiments can only be driven by solutal surface tension gradients, as stated byBassou and Rharbi in their article. The mechanism of cells formation in these experiments being clearly known, the measured velocities can be used for comparison with the numerical simulations of the solutal Bénard-Marangoni problem. Velocities measured close to the free surface at the beginning of the drying in the nonlinear quasi-steady regime are plotted in figure C.12 (the error bars give the minimum and maximum values measured during the first 1000 s). There are not enough experimental data to conclude on the validity of the scaling presented in section 3, but comparison with numerical simulations reveals that the model overestimates the experimental results by a factor 3. Even if this discrepancy does not put into questions the conclusions of the paper based on comparison of order of magnitudes, it is interesting to understand its origin. This difference seems too high to be explained by experimental uncertainties (velocity measurement)

First, different methodologies

  are proposed in a transient problem to predict the critical time for convection onset. Nonlinear and non-normal methods give similar results, within the uncertainty inherent to any transient problem. We then focus on two configurations unstable for both the thermal and the solutal problem. The comparison of critical times, which is the information that a linear approach could give, does not allow to discriminate. Nonlinear 2D simulations show that temperature gradients produce higher velocities than concentration gradients. We have thus exhibited configurations where thermal mechanism is dominant, while a simple comparison of stability thresholds would have lead to the reverse conclusion. When comparing to experimental results, a very good agreement between simulated thermal wavelengths and experimental ones are obtained. Finally, using scaling arguments, we derive a general description of the nonlinear quasi-steady regime for the Bénard-Marangoni configuration. Correlations have been obtained for the thermal and solutal problem and validated by numerical simulations. One important result is the estimation of the ratio of the thermal velocity to the solutal velocity (Eq.40), which gives a simple relation to easily estimate the relative importance of the two mechanisms. It would be useful to derive the same type of relation for other systems than polymer solutions, for instance water/ethanol, where both components are volatile. In the present work thermal and solutal mechanisms have been decoupled, in order to highlight their respective roles for experiments where both mechanisms lead to unstable configuration. Next step should concern the development of a complete thermal-solutal model to put forward coupled effects. This would be necessary to analyse longer times, when the increase of viscosity due to the decrease of solvent concentration induces a coupling between the thermal and solutal problems (this can result in the formation of a crust, observed in experiments by Toussaint et al.

section 6 . 3 (

 63 PS/toluene solution of initial thickness d i = 1.4 mm, test cases 3 to 5), where the solutal Bénard-Marangoni problem is the only one to be unstable. The polymer diffusive flux j is the sum of the flux j f = -Dρ∇ω p driven by the concentration gradient (Fick's law) and the flux j s = -DρS T ω p ω s ∇T driven by the temperature gradient (Soret effect). S T is the Soret coefficient, ω p and ω s are the mass fractions of polymer and solvent, respectively. With the assumption of constant density we get: j = j f + j s = -Dρ∇ϕ p -DρS T ϕ p ϕ s ∇T (C.1)

Figure C. 1 :

 1 Figure C.1: Experimental evaporation flux as a function of time (experimental conditions corresponding to test case 2 in table1). The right figure is a focus on the initial stages. Times t 1 , t 2 , t 3 (vertical dashed-dotted lines) separate the different regimes. The horizontal dashed line represents the plateau for t 1 t t 2 .

Figure C. 2 :Figure C. 3 :

 23 Figure C.2: Formation of thermal convection cells observed by infrared camera (experimental conditions of test case 2 in table 1). The dimensional time corresponding to each image is indicated in the lower right corner.

Figure C. 4 :

 4 Figure C.4: Thermal model for test case 2 in table 1: difference < T (z = 0, t) > -< T (z = d, t) > of the horizontal mean temperatures as a function of time. The thick line represents the pure diffusive case (obtained by setting α th = 0 and γ th = 0 i.e. Ra th = M a th = 0), other lines represent convective cases for different initial amplitudes r * for temperature perturbations.

Figure C. 5 :

 5 Figure C.5: Solutal model for test case 2 in table 1: difference < ϕs(z = 0, t) > -< ϕs(z = d(t), t) > of horizontal mean solvent volume fractions as a function of time. The thick line represents the pure diffusive case (obtained by setting α sol = 0 and γ sol = 0 i.e. Ra sol = M a sol = 0), other lines represent convective cases for different initial amplitudes r + for solvent volume fraction perturbations.

Figure C. 11 :

 11 Figure C.11: Solutal model: quasi-steady regime in the BM configuration. Domains derived from the scaling analysis (the figure is drawn for Sc = 2.3 × 10 5 ).

Figure C. 12 :

 12 Figure C.12: Solutal Bénard-Marangoni problem: average dimensionless velocity at the free surface. Black open symbols: L2 norm from numerical simulations (•: PIB/toluene, ϕ pi = 0.047, A = 5, constant viscosity ; : PIB/toluene, A = 10, variable viscosity ; ⋄: PS/toluene, A = 30, variable viscosity ; Straight line: scaling law (B domain, prefactor fitted over numerical simulations) ; Error bars: experimental data from Bassou and Rharbi [2].

Table 1 :

 1 Experimental configurations with d i the initial layer thickness, ϕ pi the initial polymer volume fraction, µ i the initial solution viscosity, H th the heat transfer coefficient, ∆Tst the characteristic temperature variation, vev the evaporation velocity. The corresponding dimensionless parameters for both theoretical models are also given (see definitions in sections 3.1 and 3.2). Physical properties of PIB/toluene and PS/toluene solutions are listed in Appendix A. Experiments are performed at room temperature.

	t 1 is associated to a thermal transient

  By contrast to the thermal model, this model was used previously to study solutal convection over time horizon t = t 2 , so it takes into account the variation of viscosity µ with respect to the the polymer volume fraction ϕ p . For the two polymer solutions (PIB/toluene and PS/toluene) used to compare models to experiments, viscosity

varies over several orders of magnitude (see figure C.3). This dependency is represented by an empirical law log 10

Table 2 :

 2 t * on when thermal and solutal problems are both unstable Time (in seconds) corresponding to the onset of convection: comparison of 2D numerical simulations and experiments for cases 1 and 2 experimental visualizations are compared in table 2. The time intervals given for models correspond to different initial perturbation amplitudes. In experiments, the time visualization at which the pattern becomes clearly visible is estimated similar by iriodin particles or infrared camera. It corresponds to an upper bound of the time of onset.

	From figure C.8, one can deduce that test cases 1 and 2 in table 1 (indicated
	by green circles in figure C.8) are clearly unstable for both thermal and solutal
	problems. The times of convection onset obtained by numerical simulations or

Despite the large Lewis number (Le ≡ κ/D ∼ 10 3 ), models show that solutal convection appears at least concomitantly and generally before the thermal one. This is related to the fact that the two cases are far from the solutal Marangoni stability frontier but close to the thermal one (figure C.8). Time of onset for the solutal model is consistent with the experimental estimation for both test cases. The thermal model is also consistent with test case 1 experiment, and the discrepancy with test case 2 is very weak. This analysis is hence not sufficient

Table 3

 3 Several configurations have been considered, all in the B domain. They are based on PIB/toluene or PS/toluene systems, assuming variable or constant viscosity and different values of aspect ratio A. Figure C.12 shows that in all cases, results of numerical simulations compare very well with the scaling law, over three decades of M a sol . The free surface velocity prefactor deduced by fit is about 0.2.

summarizes the overall results obtained for velocity and temperature or concentration variations as a function of the relevant dimensionless numbers. Different domains are delineated: domain A (no convection) and domains B, C, D, E which correspond to different cases for δ H , δ T and δ S . The frontiers delimiting such domains are given in table 4. An illustration of the different domains is given in figure C.11 for the solutal problem and Sc = 2.3×10 5 (value of Sc corresponding to test case 2).

Solutal scaling laws have been tested via numerical simulations by getting, at the beginning of the quasi-steady state, the L2 norm of the free surface velocity.

Table 4 :

 4 Quasi-steady regime in the BM configuration -boundaries between the different domains.

  For consistency, this equation is only valid when S c /v x ≪ 1. When S c /v x 1, the hydrodynamic boundary layer thickness saturates at δ H ∼ 1, and inertia is no more involved in the problem.

	magnitude. The continuity equation (B.1) coupled to hypotheses H1 and H4
	(see section 6) imposes
	v z ∼ δ H v x	(B.10)
	From hypoteses H1 and H4 and equation (B.5) and assuming δ H ≪ 1 one easily
	gets
	v x ∼	Sc H δ 2	(B.11)
	From equation (B.8), a relation between ∆φ s
			9)
	Scaling laws are obtained by solving equations B.1 to B.9 in terms of order
	of magnitude. From now on, all quantities are thus intended to be orders of

  Case 1' implies vz ∼ v z . The balance between the different terms of (B.13) plus equation (B.10) leads to the same equation (B.14). Case 2 implies vz ∼ δ S v z . The balance between the different terms of (B.13) plus equation (B.10) leads again to equation (B.14). Sc 1, this means that δ S δ H (for δ H ∼ 1, this is obviously true). Note as well that this implies that case 1' is not possible. Due to the fact that φ s ≤ 0 and ∂ z φ s ≤ 0, the boundary condition (B.7) For the other cases, the method is similar. The synthesis of all the cases are presented in tables 3 and 4. The different domains are displayed in figure C.11.

		v x ∼ S 1/4 c M a	1/2 sol	(B.21)
		δ H ∼ S 3/8 c M a	-1/4 sol	(B.22)
	When equation (B.14) is compared with equation (B.11) (valid only when δ S ∼ S -1/8 c M a -1/4 sol (B.23)
	δ H ≪ 1), this imposes	∆φ s ∼ S -1/8 c	M a	-1/4 sol .	(B.24)
	δ S /δ H ∼ S -1/2 c These relations should be compatible with the assumptions which imposes when δ H ≪ 1 (B.15)
	δ H ≪ 1 Since we assume that expresses the balance between only two terms: ⇒ S 3/2 c ≪ M a sol δ S P e int ≪ 1 ⇒ P e 4 int S -1/2 c ≪ M a sol	(B.25) (B.26)
		∆φ s
	v x ∆φ s ; vz where vz denotes the magnitude of vertical velocity within the solutal boundary ∆φ s δ S ; ∆φ s δ 2 S (B.13) layer. This quantity can be different from v z the magnitude of vertical velocity within the hydrodynamic boundary layer. By definition δ H 1 and δ S 1. Three cases can then be considered δ Case β : δ S P e int ∼ 1 (B.19)
	Case 1 : δ H ≪ 1 We must now consider the four combinations of cases 1 or 2 and cases α or β and δ S δ H
	Case 1 ′ : δ H ≪ 1 taking into account that δ S δ H is always valid. In the following, we only and δ H δ S
	Case 2 : δ H ∼ 1 derive explicitly the single case (1 and α). Let us assume that cases 1 and α are and δ S δ H
	valid. This means
	Case 1 implies δ S	δ H and thus vz ∼ δS δH v z . The balance between the δ H ≪ 1 and δ S Pe int ≪ 1 (B.20)
	different terms of (B.13) plus equation (B.10) leads to Using equations B.11, B.12, B.14 and B.18, it is possible to get expressions for
	v x , δ S , δ H and ∆φ s	v x ∼ δ -2 S				(B.14)

S ∼ (1 -P e int φ s ) (B.16) Assuming φ s ∼ 0 in the bulk (quasi-steady assumption), we get ∆φ s ∼ -φ s and the boundary condition (B.7) reads ∆φ s ∼ (δ S + P e int δ S ∆φ s ). (B.17) This equation cannot be verified if P e int δ S ≫ 1, so only two cases are possible Case α : δ S P e int ≪ 1 then ∆φ s ∼ δ S (B.18)

  ). The right figure is a focus on the initial stages. Times t 1 , t 2 , t 3 (vertical dashed-dotted lines) separate the different regimes. The horizontal dashed line represents the plateau for t 1 t t 2 .

Acknowledgments: TB and MR acknowledge financial support from the Deutsche Forschungsgemeinschaft in the framework of the priority program "Transport Processes at Fluidic Interfaces" (grant Bo 1668/6). FD and BG acknowledge financial support from the European Union via the FP7 Marie Curie scheme (grant PITN-GA-2008-214919 (MULTIFLOW)). FD, BT, EC, CD and BG acknowledge financial support from the Fédération TMC (Ministère de l'Enseignement supérieur et de la Recherche, France).

Appendix A. Physical properties of PIB/toluene and PS/toluene solutions

Viscosity

To take into account the variation of the solution viscosity with solvent concentration, we use the empirical law interpolated from measurements (data from Gorand et al. [START_REF] Gorand | Instabilités de plissement lors du séchage de films polymères plans[END_REF], Mark [START_REF] Mark | Polymer Data Handbook[END_REF] for PIB/toluene, Bassou and Rharbi [START_REF] Bassou | Role of Bénard-Marangoni instabilities during solvent evaporation in polymer surface corrugations[END_REF] for PS/toluene):

with Y = log 10 (ϕ p ), (a 0 , a 1 , a

Mutual diffusion coefficient

To the best of our knowledge, there is no available data in the literature for the mutual diffusion coefficient of PIB/toluene solutions at room temperature, for the concentration range covered by our experiments. Nevertheless, a rough estimate D ≃ 10 -10 m 2 s -1 can be obtained by extrapolating measurements performed in the concentrated regime by Doumenc et al. [START_REF] Doumenc | Mutual diffusion coefficient and vapor-liquid equilibrium data for the system pib/toluene[END_REF].

On the contrary, the mutual diffusion in PS/toluene system has been widely [START_REF] Rauch | Collective and thermal diffusion in dilute, semidilute, and concentrated solutions of polystyrene in toluene[END_REF] for D, the two last references also provide S T ). The uncertainty range of D represents the scattering over the different sources. The differences between S T values from Zhang et al. [START_REF] Zhang | Thermal and mass diffusion in a semidilute good solvent-polymer solution[END_REF] and Rauch and Kohler [START_REF] Rauch | Collective and thermal diffusion in dilute, semidilute, and concentrated solutions of polystyrene in toluene[END_REF] are lower than the experimental uncertainty, which is around 10 -3 .

Other physical properties

Some properties of PIB/toluene and PS/toluene solutions are assumed to be equal to those of pure toluene (from Riddick et al. [START_REF] Riddick | Organic solvents, physical properties and methods of purification[END_REF], Monteil and Postel [START_REF] Monteil | Techniques de l'ingénieur[END_REF])

Specific properties of the PIB/toluene solution are the following ones (cf. [START_REF] Toussaint | Experimental characterization of buoyancy-and surface tension-driven convection during the drying of a polymer solution[END_REF][START_REF] Mark | Polymer Data Handbook[END_REF] for details):

PIB molar mass: M w = 500 kg.mol -1 , ρ P IB ≃ 920 kg.m -3 , α sol = 5.82 ×

Specific properties of the PS/toluene solution are the following ones (cf. [START_REF] Bassou | Role of Bénard-Marangoni instabilities during solvent evaporation in polymer surface corrugations[END_REF][START_REF] Mark | Polymer Data Handbook[END_REF] for details):

PS molar mass:

Critical time

Frozen-time Non-normal approach -G thres =1

Non-normal approach -G thres =10 2

Non-normal approach -G thres =10 4 Nonlinear approach given by the frozen-time approach ; t * (nn) on given by the non-normal approach with E T norm and different thresholds G thres ; t * (nl) on given by the nonlinear approach for different perturbation amplitudes (10 -7 ≤ r * ≤ 10 -2 for M a th = 2 × 10 4 , 10 -9 ≤ r * ≤ 10 -4 for M a th = 9 × 10 4 and M a th = 2.5 × 10 5 ). The other dimensionless parameters are such that Ra th = 0, Bi = 0.01. For non-normal and frozen-time methods, Prandtl number is set to P r = ∞, for the nonlinear approach, it is set to P r = 100.

). The critical zone of the thermal Rayleigh-Bénard-Marangoni problem corresponds to the extreme boundaries obtained by the non-normal method formed by curves G thres = 1 or G thres = 100 using E V or E T norms (see Doumenc et al. [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni convection due to evaporation : a linear non-normal stability analysis[END_REF]). The critical conditions of the solutal Bénard-Marangoni problem and the solutal Rayleigh-Bénard problem have been obtained by the nonlinear approach with r + = 10 -6 /∆ϕ (see Trouette et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF]). Changing r + by ± 2 decades results in a maximum variation of the critical thickness of the order of 30 % (not shown on the plot). Experimental test cases 1 and 2 are indicated by green circles and arrows.