Antoine Faugaret

Yohann Duguet

Yann Fraigneau

Laurent Martin Witkowski

Influence of interface pollution on the linear stability of a rotating fluid

Keywords: Rotating flow, free surface, experiment, numerical simulation

The boundary conditions at a liquid-gas interface can be modified by the presence of pollutants. This can in turn affect the stability of the associated flow. We consider this issue in the case of a simple open cylindrical cavity flow where a liquid is set in motion by the rotation of the bottom. The problem is addressed using an experimental set-up, a linear stability code and direct numerical simulation. A robust mismatch between numerical and experimental predictions of the onset of instability is found. We model the possible effect of unidentified pollutants at the interface using an advection-diffusion equation and a closure equation linking the surface tension to their concentration. The chosen closure is inspired by studies of free surface flows with surfactants. Numerical stability analysis reveals that the base flow and its linear stability threshold are strongly affected by the addition of pollutants. Pollutants tend to decrease the critical Reynolds number, however the nonlinear dynamics is less rich than without pollutant. For sufficiently high pollution levels, the most unstable mode belongs to a different family, in agreement with experimental findings.

Introduction

Fluid flow simulations rely on a mathematical formulation associating given governing equations to specific boundary conditions. The choice for the boundary conditions is sometimes not trivial, in particular in the presence of a liquid-gas interface. Beyond the difficulties stemming from a deformable interface, it appears that in practice the correct boundary conditions are not well known even for a perfectly flat interface. The classical boundary condition considered in textbooks is commonly deduced from the balance of tangential stresses at the interface. For a gas-liquid interface, where the dynamic viscosity of the gas is negligible with respect to that of the liquid, this leads to a "free slip" condition, which is simple to implement in simulation codes. Unfortunately this ideal boundary condition is not necessarily representative of realistic experiments, even for liquids as common as plain water. Contamination by pollutants present in the ambiant air can influence the rheology of the interface and drastically impact the effective boundary conditions. Such modifications can deeply alter the flow, as a consequence the numerical predictions with a free surface condition are no longer representative of the true physical flow. Considering how difficult it is to experimentally ensure that a gas-liquid interface is free from any chemical pollution, it is crucial to know how to model the interface, without necessarily knowing all the physical properties in detail. Such issues arise for instance in flow where Marangoni effects (modifications of the surface tension due to e.g. temperature effects) may interfere and impact the flow dynamics. The simplified phenomenology of surface pollutants assumes that, although the precise chemical composition of the pollutants is by essence unknown, their qualitative effect is to reduce the effective surface tension. This suggests an effect akin to that of surfactants added on the free surface. As in studies of free-surface flow with controlled amounts of surfactants, we assume for simplicity that the amount of pollutants is small enough so that they are confined to a Langmuir monolayer located directly at the interface, and do not penetrate the bulk. Advected by the local velocity field tangent to the interface, the pollutants cluster at some given locations, their accumulation being only resisted by weak diffusion. The resulting inhomogeneity of the concentration field at the surface induces a local change in the surface tension. The gradients of effective surface tension lead to additional stresses that modify the global stress balance.

In this investigation we choose a flow case feasible in the laboratory as well as in numerical simulations, where such ideas can be tested. In particular, we focus on simple flow likely to develop instability modes via a classical Hopf bifurcation scenario. The selected most unstable mode, its growth rate and the associated onset Reynolds number serve as quantitative indicators of how reliable a given set of boundary conditions are. The flow consists of a cylindrical cavity partially filled with liquid, in most cases water. The top of the cavity is open while its bottom rotates with a constant angular velocity. The side walls do not rotate and are fixed in the laboratory frame. For simplicity, we restrain ourselves to the parameter regime where the fluid interface remains approximately flat even as the instability develops and saturates. A sketch of the experimental set-up can be found in Fig. 1. The two main parameters for this flow are the geometric aspect ratio G = H/R, where H is the undisturbed liquid height and the inner cylinder radius R, and the Reynolds number Re = ΩR 2 /ν where Ω is the rotation rate and ν is the kinematic viscosity. This flow has been previously studied both numerically and experimentally. The earliest publication we found about this configuration is a numerical investigation of the base flow for an aspect ratio G between 0.1 and 1 [START_REF] Hyun | Flow in an open tank with a free surface driven by the spinning bottom[END_REF]. Under the assumption that the flow remains axisymmetric, the transition to unsteadiness has been studied numerically for G=2 by [START_REF] Daube | Numerical simulation of axisymmetric vortex breakdown in a closed cylinder[END_REF]. The transition point was found near Re=2975. Evidence for an instability breaking the axisymmetry of the base flow was given only later. In [START_REF] Young | Period-doubling route to chaos for a swirling flow in an open cylindrical container with a rotating disk[END_REF], no visualisation of the pattern was shown, however laser Doppler velocity measurements revealed the growth of an instability near Re = 2000 for a G=2 geometry. The first experimental visualisations of non-axisymmetry were performed by Hirsa et al. (2002b); [START_REF] Lopez | Symmetry breaking in free-surface cylinder flows[END_REF] in the same geometry, and compared with numerical results for G=2 and G=1/4 in [START_REF] Lopez | Symmetry breaking in free-surface cylinder flows[END_REF]. For the larger aspect ratio (G=2), numerical predictions and experimental results tend to agree, yet for the shallower configuration, mismatches in critical Reynolds number and azimuthal wavenumber m were reported. In particular, the wavenumber selection was described in these works as sensitive to the presence of surface contaminants. Among recent publications, the experimental work by [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] surveys many aspect ratios ranging from G=0.0179 to 0.107. Larger aspect ratios G from 0.3 to 4 have been studied numerically as well [START_REF] Cogan | Symmetry breaking and instability mechanisms in medium depth torsionally open cylinder flows[END_REF][START_REF] Iwatsu | Analysis of flows in a cylindrical container with rotating bottom and top underformable free surface[END_REF][START_REF] Serre | Coupled numerical and theoretical study of the flow transition between a rotating and a stationary disk[END_REF]. For higher rotation rates, a different regime takes over, with strong deformations of the interface and sometimes mode switching [START_REF] Suzuki | Surface switching of rotating fluid in a cylinder[END_REF][START_REF] Tasaka | Flow transitions in the surface switching of rotating fluid[END_REF]. Polygonal patterns at the deformed interface have been reported by [START_REF] Vatistas | Experiments on waves induced in the hollow core of vortices[END_REF]; [START_REF] Jansson | Polygons on a rotating fluid surface[END_REF]; [START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF] and modelled by [START_REF] Tophøj | Rotating polygon instability of a swirling free surface flow[END_REF].

In the present investigation, we revisit the primary instability mechanism using a joint experimental and numerical approach. We focus on the primary instability in the case of an approximately flat interface. For small enough angular velocities the centrifugal acceleration remains much smaller than gravity and the curvature of the fluid interface can be indeed neglected in the small Froude number hypothesis. The main aspect ratio under scrutiny corresponds to G=1/14. As shown in Table 1, the experimentally determined thresholds are lower by least 75% than those of [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF]. The various possible reasons for this discrepancy have been reviewed in our experimental set-up with great care, among them residual vibrations, lack of axisymmetry of the cavity, finite curvature of the free surface, presence of a meniscus, ionisation of the water. In all cases these hypotheses were ruled out as quantitatively insignificant. Note that quantitative discrepancies with experimental measurements have been also already reported earlier for this flow for low G. In [START_REF] Kahouadji | Seuils de stabilité pour un écoulement à surface libre engendré dans une cavité cylindrique tournante à petit rapport de forme[END_REF], the stability thresholds in Re determined by LSA were compared with Poncet and Chauve's experimental estimates for varying values of G. In both studies the threshold value Re c increases with decreasing G. While the agreement between numerics and experiments is very satisfying, it suddenly deteriorates for G 0.07-0.08 (see fig. 3 in [START_REF] Kahouadji | Seuils de stabilité pour un écoulement à surface libre engendré dans une cavité cylindrique tournante à petit rapport de forme[END_REF]), with a mismatch on Re c of 100% for G ≈0.04. Following Lopez and co-authors, we assign such a mismatch between experiments and linear stability analysis to the unavoidable presence of pollutants at the interface, and hence to the simplistic free slip model for the boundary conditions at the liquid interface. The present investigation is devoted to a quantitative analysis of the influence of these pollutants, via a simple phenomenological model, on the linear stability threshold of this flow.

The outline of the paper is as follows. In Section 2, we give a brief description of the flow and its primary instability. We detail the experimental methods as well as the numerical methods for the linear and nonlinear stability. Section 3 is devoted to a critical comparison between experimental and numerical results. In Section 4, we introduce a new model for the free surface where surface pollution is taken into account. Section 5 discusses the possible simplification of the model in the limit of high surface pollution. The final Section 6 contains a summary of the present investigation together with open questions and perspectives for future work.

Flow set-up and related investigation techniques

Base flow description

We briefly recall the main features of the base flow as described by [START_REF] Iwatsu | Analysis of flows in a cylindrical container with rotating bottom and top underformable free surface[END_REF] and [START_REF] Yang | Axisymmetric rotating flow with free surface in a cylindrical tank[END_REF]. It is axisymmetric with three non-zero velocity components. Its structure for small aspect ratio G is sketched in fig. 1. We use a classical direct cylindrical coordinate system (r, θ, z), where r is the radial distance, θ the azimuthal angle and z the distance from the rotating bottom. In the vicinity of the instability threshold, the azimuthal velocity profile possesses a simple radial structure almost independent of z except in the boundary layers. In the regimes we focus on, the azimuthal velocity increases with the radial distance from r = 0 to r ≈ 0.67 -0.68R, where R is the radius of the set-up, and decreases to zero as the wall is approached. The latter zone is Above a given rotation rate, this base flow is known to support an instability breaking its axisymmetry. Ignoring in a first stage the geometrical and rheological parameters, a simplistic explanation for this symmetry breaking is as follows : a shear instability, akin to a Kelvin-Helmhotz instability along a curved streamline, develops where the azimuthal velocity profile displays the strongest curvature. Given the cylindrical geometry, a direct analogy with the instability of Stewartson layers in the split disk configuration [START_REF] Stewartson | On almost rigid rotations[END_REF]) case has been suggested in order to justify the relative size of the instability region [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF]. Beyond this naive picture, the bifurcation scenarios leading to the presence of different non-axisymmetric patterns with a dominating wavenumber m = 0 as in fig. 2, are not entirely clear from the literature. [START_REF] Lopez | Symmetry breaking in free-surface cylinder flows[END_REF] describe the bifurcation as a standard Hopf bifurcation. Experiments in [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] reveal the existence of hysteresis, suggesting a possibility for subcritical bifurcation. In the present work, we focus on the emergence of a m=5 mode, the most unstable one as predicted by linear instability theory for the aspect ratio G considered. A competing unstable mode with m=4, though theoretically expected to appear for parameters where the mode m=5 is already unstable, has also been investigated.

Experimental technique

The main element of the experimental set-up is a cylindrical shaped plexiglas cavity. Its internal radius is R = 140.3 ± 0.05 mm, and the thickness of the plexiglas is 6.8 mm. The value of R is used to define the Reynolds number Re = ΩR 2 /ν. The cavity was engineered from a single block, so that the cylinder and the bottom are monolithic, avoiding any risk of leak. Its bottom is drilled along its axis in order to mount a brass foot, that will host the shaft of the rotating disc, itself also made of brass. The radius of the disc is R d = 139.6 mm, its thickness 8.5 mm and its mass 5 kg. The shaft is held in place with two ball bearings, and the sealing is insured by a spring-loaded double-lip seal. An aluminium rigid sleeve coupling, relying on a thrust ball bearing, is used to connect the disc shaft and the motor reducer unit. The motor used is a direct current motor (Parvex RX320E-R1100) with a 1:12 reducer. The rotation speed is controlled using a tachometer closed-loop. Special attention was paid to minimising the gap between the disc edge and the vertical wall of the cavity. The liquids used in this experimental investigation are tap water, de-ionised water, and a water-glycerol mixture. As the cavity is not thermo-regulated, fluid temperature is monitored continuously, with a digital thermometer that allows it to be known with an accuracy of 0.1K. The corresponding kinematic viscosity is then evaluated using an empirical formula [START_REF] Cheng | Formula for the viscosity of a Glycerol-Water mixture[END_REF]. The experimental Reynolds number, based on the angular velocity Ω, the radius and the kinematic viscosity, is hence known within a given accuracy of the order of a percent for the range of parameters investigated. The relative error is expected to increase as the rotation rate decreases.

Flow measurements are made using an LDV device, composed of a Dantec Laser linked to a BSAFlow processor. The liquid is seeded with Dantec 10 micrometers diameter silver coated hollow glass spheres. Because of the cylindrical geometry, as the laser beams are placed for the acquisition of u θ they undergo a deviation that both shifts the location of the focus and impacts the quantitative measurements. This is fixed at the post-processing stage using the technical corrections suggested in [START_REF] Huisman | Applying laser Doppler anemometry inside a Taylor-Couette geometry using a ray-tracer to correct for curvature effects[END_REF]. For visualisations such as in fig. 2, the flow patterns were highlighted by injecting either Kalliroscope or ink into the fluid. The experimental protocol is based on a waiting time of five to ten Figure 3: Instability growth for a m = 5 mode for a water + Kalliroscope mixture initially at rest. G=1/14, Reynolds number Re=16 550. 3a) t = 52s, 3b) t = 60s, 3c) t = 71s, 3d) t = 99s.

minutes after each increase of the rotation speed, followed by LDV acquisition performed over another five minutes duration.

Experimental evidence for m=5 instability

We describe the experimental instability leading to a steadily rotating m=5 mode, using Kalliroscope visualisations or pointwise LDV measurements. The initial configuration of the flow is rest. The angular velocity is directly set to a finite value defining the target Reynolds number Re, the value used in Fig. 3 being Re=16 650. For low enough Reynolds numbers, the flow remains axisymmetric as in figure 3a. With the increase in Re, an annulus characterised by stronger shear appears around r = 0.7R (fig. 3b). An m=5 mode emerges (fig. 3c) and evolves towards a steadily rotating configuration with 5 co-rotating vortices (figure 3d). A similar scenario occurs for other values of m, in particular m = 4 which has been observed for other nearby values of Re. The vortex pattern rotates with a constant angular velocity smaller than Ω. The angular frequency f of the pattern can be deduced using f =2πf d /(mΩ), where f d is the dimensional frequency obtained experimentally using pointwise LDV measurements at a location fixed in the laboratory frame. The main frequency f varies moderately over the range Re=[4230,16300]. A Fourier transform of the time series is shown in fig. 4 in the case m=5 for Re=4230 and Re=16300. The main frequency and related harmonics dominate the spectrum.

Numerical methodology for free slip interfaces

As a complementary part of this investigation, we have used numerical tools based on the incompressible Navier-Stokes in order to investigate both linear and nonlinear aspects of the symmetry-breaking instability. The present section first introduces the numerical methods used. It also features a comparison with the experimental results of Subsection 2.3.

Mathematical model

We adopt the point of view of a single-phase flow. The velocity field u(r, θ, z, t) inside the liquid is governed by the incompressible Navier-Stokes equations in the (non-rotating) laboratory frame :

∇ • u = 0, ( 2
∂u ∂t + (u • ∇)u = -∇p + 1 Re ∇ 2 u. (2.2) Equations (2.1)-(2.
2) have been non-dimensionalised using the length scale R and the time scale Ω -1 , and the dimensionless fluid density is taken as unity.

The flow obeys no-slip at all solid boundaries. This implies

u r = u θ = u z = 0 (2.3)
at the fixed vertical boundary at r = 1, whereas

u r = u z = 0, u θ = r (2.4)
on the rotating disk at z = 0.

The boundary condition at the liquid-gas interface at z = G is classically derived from the stress balance at the interface. As the viscosity of the air is much smaller than the water one, we can neglect the gas phase altogether. We first consider the generic free-slip boundary conditions :

∂u r ∂z = 0, ∂u θ ∂z = 0, u z = 0 at z = G. (2.5)
As in all mesh-based numerical methods, the singularities of the velocity field occurring at both corners (r = 1, z = 0) and (r = 1, z = G) are smoothed out in practice by the finite mesh without the need, as for spectral methods, for regularising functions [START_REF] Serre | Vortex breakdown in a cylinder with a rotating bottom and a flat stress-free surface[END_REF] or singular splitting [START_REF] Duguet | Instability inside a rotating gas cylinder subject to axial periodic strain[END_REF]. This is consistent with the "natural" regularisation occurring in the experiment in the presence of a very thin gap.

Linear stability analysis

In order to determine the critical Reynolds number Re c for the onset of instability, we use an in-house linear stability solver named ROSE, based on a finite difference method in r and z. The technique as well as the equations written in cylindrical coordinates are found in [START_REF] Kahouadji | Thermocapillary instabilities in a laterally heated liquid bridge with end wall rotation[END_REF]. The steady axisymmetric base flow is first determined by solving Eqs. 2.1, 2.2 with the associated boundary conditions using a Newton-Raphson solver. The steady solution is solved for in an (ω, ψ, u θ , c) formulation, where

ω = ∂u r ∂z - ∂u z ∂r , u r = 1 r ∂ψ ∂z , u z = - 1 r ∂ψ ∂r .
(2.6)

The Newton-Raphson solver allows for additionnal scalar fields c(r, z) such as temperature or concentration further discussed in section 4.

Let (U, P ) represent the velocity-pressured field for such a steady axisymmetric solution of Eqs. 2.2, and let (u * , p * ) be a small-amplitude perturbation to (U, P ). The dynamics of the perturbation is governed by the linearised stability equations

∇ • u * = 0,
(2.7)

∂u * ∂t + (U • ∇)u * + (u * • ∇)U = -∇p * + 1 Re ∇ 2 u * . (2.8)
It is associated with Dirichlet boundary conditions u * = 0 on all solid boundaries together with a boundary condition on u * at the interface similar to that for u in Eq. 2.5. All meshes used are Cartesian in the meridional plane (O, r, z). For G=1/14 the mesh consists of 701 × 101 grid points. Numerical comparison with [START_REF] Kahouadji | Seuils de stabilité pour un écoulement à surface libre engendré dans une cavité cylindrique tournante à petit rapport de forme[END_REF] confirms that this resolution is sufficient for an estimation of Re c with an accuracy below one percent.

Direct numerical simulation

For the nonlinear validation of stability thresholds we have used the direct numerical simulation (DNS) code Sunfluidh developed at LIMSI for incompressible flows. It is based on a projection method to ensure a divergence-free velocity field. The equations are discretized on a staggered structured non-uniform grid using a finite volume approach with a second-order centered scheme in space. A second-order Backward Euler Differentiation is used for time discretization. Details can be found in [START_REF] Yang | Large axisymmetric surface deformation and dewetting in the flow above a rotating disk in a cylindrical tank: spin-up and permanent regimes[END_REF]. The interface condition is as in Eqs. 2.3 and 2.4. The code offers the possibility to enforce a given rotational symmetry R m characterised by a fundamental azimuthal wavenumber m 0, such that every velocity field verifies

(R m u)(r, θ, z) = u(r, θ + 2π m , z) = u(r, θ, z) (m = 0), (2.9)
or axisymmetry for m=0. In that case the simulation only needs to be carried out over an angular sector 0 θ 2π/m with azimuthal periodicity. For the simulations without symmetry imposed, we have used a mesh consisting of 180 × 180 × 64 cells in r, θ and z.

Critical comparison of the different approaches

Comparison between the numerical methods

For identical parameters, we report excellent agreement between the base flows computed by the two methods for all Re. Whereas the base flow can be converged for all Re using the Newton method, it is only accessible for Re < Re c using time-stepping. However since the base flow is apparently the only axisymmetric solution of the system, it is also found using DNS for all Re by simply imposing m=0 (2D axisymmetric case) and stepping forward in time. The value of Re c for m=5 is first identified by LSA using a secant method. In the DNS code, the procedure used to identify Re c is different : above and below Re c , an arbitrary perturbation of finite but small amplitude is applied to the

Experiment (LDV) Poncet ROSE Sunfluidh Rec [3160 -4230] [14367 -16420] 17006 [17000 -17100] fm=5 [Noise -0.764] - 0.709 [0.707 -0.699]
Table 1: Critical Reynolds number and angular frequency of the pattern for the m=5 instability for G=1/14. For Sunfluidh, linear interpolation leads to Re c = 17010. When relevant, a lower and an upper bound are given for Re c , with the corresponding values for the frequency. Experimental results by [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] are included for comparison.

system after an initial transient, with the R 5 symmetry imposed or not. This impulse response leads to either exponential decay towards the base flow, or exponential growth towards a nonlinear regime at large times. A linear interpolation of these rates leads to an evaluation of the critical threshold Re c . Both approaches agree quantitatively very well regarding the prediction of Re c for m = 5 since the relative error is close to 0.3% (see Table 1). Interestingly, this comparison, as well as the lack of unstable impulse response for Re < Re c (even for larger amplitude impulses), both suggest that the instability is not subcritical, at least for a clean interface obeying the boundary condition (2.5).

Mean flow structure

Since both numerical approaches yield a truly similar base flow solution, comparison with the experimental base flow measured using LDV would be relevant at this point. As we shall see, measurements below Re c turn out to be experimentally difficult. Another comparison, easier to perform, concerns the mean velocity profiles obtained for Re > Re c by either temporal or spatial average. Such a comparison is displayed in Figure 5. For the eigenmodes computed using ROSE, their average is by construction zero. Hence only the base flow obtained by LSA is included in Figure 5, whereas the spatial average is taken for the DNS data and the time average for experimental LDV data. A common value of Re = 18620 is chosen for the comparison. Although the agreement is by and large satisfying, a noticeable overshoot appears around r ≈ 0.67 in all numerical azimuthal velocity profiles, with no equivalent in LDV measurements despite sufficient measurement accuracy. The same absence of overshoot area was also reported for comparable parameters in [START_REF] Yang | Axisymmetric rotating flow with free surface in a cylindrical tank[END_REF].

Threshold detection

The most dramatic mismatch between numerics and experiments concerns the critical Reynolds number. While both numerical simulations agree on a critical Reynolds number around 17000 (see Table 1), LDV measurements display persistant oscillations in the azimuthal velocity field for Re as low as 4200, with a normalised frequency f 5 = 0.76, indicative of the presence of the m=5 mode. This upper bound on the value of Re c is smaller by a factor of 4 than the previous experimental estimates by [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF]. These values can be found in Table 1. These discrepancies are robust : although the exact same spin-up protocol as [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] was observed, the respective ranges of Re c differ. We note that the threshold detection by [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] is based on Kalliroscope visualisations. Kalliroscope appears in our set-up as a poor diagnostic for Re c for this flow case : the threshold detection is erratic and protocol-dependent. Indeed the estimation of Re c fluctuates between 6200 to 9300. . Its spatial structure is directly comparable to that of the saturated flow above Re

(LSA) c displayed in fig. 3d.
On one hand there is perfect numerical agreement between LSA and DNS about the estimation of Re c , on the other hand there is a troubling match with [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] at odds with the experimental/numerical discrepancy we report. We have hence carried out an exhaustive investigation of the possible reasons for such a discrepancy by focusing on experimental imperfections. A classical reason for discrepancies in rotating machines is the presence of mechanical noise that could force an instability by direct or parametric resonance. The experimental displacement of the disc surface was measured using a pair of LK-G10 sensors, and their associated LK-GD500 controller. Displacements were evaluated to approximately 10 -4 m, with a mean frequency corresponding to the disc rotation, yet no link with the pattern frequency was found. This does not suggest any obvious experimental flaw in our experimental methodology.

Eventually, in order to confirm our experimental approach, we switch temporarily to a different geometry with G=1/4 where a direct and favourable comparison with the experimental results of [START_REF] Lopez | Symmetry breaking in free-surface cylinder flows[END_REF] can be made. For these parameters there is also a robust mismatch between experiments and numerics : LSA predicts the most unstable mode m = 2 with Re c = 3500, whereas Lopez's experiments at Re = 2000 show a mode m = 3, which is predicted using LSA to be unstable only for Re 4600. We have then conducted our own experiments with two different mixtures of water with 20% and 55% glycerol. The motivation for these two different mixtures is to allow for a wider span of rotation speeds : using water the instability would have occurred for rotation speed below 1 rpm where the signal-to-noise ratio in the LDV degrades. In both cases the mode m=3 is detected, either using ink or LDV, for Re = 2160 in the 55% glycerol fluid (see Fig. 2a) and Re = 2520 in the 20% glycerol fluid. All the results are gathered in table 2. This side study confirms, in good agreement with [START_REF] Lopez | Symmetry breaking in free-surface cylinder flows[END_REF] that numerics overestimate the experimental thresholds in Re. The discrepancy reported here for G=1/14 has hence a robust physical origin, which the rest of this paper is devoted to.

Nonlinear dynamics

The mismatch between numerics and experiments for G=1/14 is even more dramatic further above Re c . Although DNS initially displays an azimuthal wavenumber m=5 close to Re c (see fig. 6a at Re=17100), the instability pattern evolves towards m=7 as Re is pushed to 18620, less than 9% above Re c . [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] have also reported an evolution of the modal content of the flow with Re, yet with m decreasing as Re is increased. A similar decrease of m with Re was also observed qualitatively in our experiment for Re sufficiently higher than Re c . However, the wavenumber m=5 remains experimentally stable from Re=4200 to at least Re=18620. The frequency spectrum is shown in fig. 7a for Re=18620. Given such a mismatch, larger values of Re were not investigated neither experimentally nor numerically. Differences in the nonlinear dynamics for Re=18620 also emerge in velocity measurements : while experimental time series display a single frequency, the signals from DNS display a broader spectrum and richer dynamics, see fig. 7. In addition to the mismatch in the modal behaviour between experiments and DNS, the vorticity patterns (figs 6a and 6b) do not match the experimental figures (2c and 3d) very convincingly. This raises doubts about whether the mode predicted in the numerics does indeed correspond to the structure observed experimentally. examined by [START_REF] Peaudecerf | Traces of surfactants can severely limit the drag reduction of superhydrophobic surfaces[END_REF] in a channel flow with superhydrophobic surfaces, in presence of carefully added surfactants. As it is nearly impossible, in standard laboratory conditions, to achieve an experiment with a perfectly clean interface at all times, it is necessary to take additional effects into account in order to properly model the behaviour of the fluid at a realistic liquid-gas interface. Previous publications with a similar experimental set-up, in which the adsorption of pollutants at the interface was carefully controlled, already demonstrated the crucial influence of pollution of the base flow [START_REF] Hirsa | Measurement and computation of hydrodynamic coupling at an air/water interface with an insoluble monolayer[END_REF](Hirsa et al. , 2002a)). There, pollutants were assimilated to a vitamin K 1 monolayer, considered as a surfactant.

Limitations of the clean interface hypothesis

In the next section, we model explicitely the presence of pollutants at the interface in the Navier-Stokes equations and investigate its qualitative as well as quantitative influence on the linear stability of the flow.

Modelling of interface pollution

Modification of the effective surface tension

The present modelling of the pollution at the interface is directly inspired by the modelling in [START_REF] Hirsa | Measurement and computation of hydrodynamic coupling at an air/water interface with an insoluble monolayer[END_REF] and [START_REF] Kwan | A mathematical and numerical study of incompressible flows with a surfactant monolayer[END_REF]. Let c d (r, θ, t) be the instantaneous concentration of the pollutants at the interface. We assume that pollutants are advected by the velocity field of the fluid while diffusing with a simple nondimensional diffusion coefficient D s . Moreover, we assume that no transport occurs from the surface to the bulk of the flow, so that the bulk concentration can be neglected [START_REF] Bandi | Hydrodynamic signatures of stationary Marangoni-driven surfactant transport[END_REF]. Starting from the Boussinesq-Scriven surface fluid model for a Newtonian fluid-gas interface [START_REF] Scriven | Dynamics of a fluid interface: Equation of motion for Newtonian surface fluids[END_REF], and under the hypothesis of negligible surface dilatational viscosity and surface shear viscosity [START_REF] Hirsa | Measurement and computation of hydrodynamic coupling at an air/water interface with an insoluble monolayer[END_REF], the boundary conditions can be written as

∂u r ∂z = 1 Ca ∂ σ ∂r , ∂u θ ∂z = 1 Ca 1 r ∂ σ ∂θ , u z = 0, (4.1)
where Ca = µΩR/σ 0 , σ = σ/σ 0 , with σ 0 the reference surface tension of the solvant (for water, σ 0 = 72.8 mN.m -1 ) and µ the dynamic viscosity.

Eq. (4.1) is not closed until the effective surface tension σ is linked to the surface concentration c d . Several closures have been suggested in the literature, all characterised by a decrease of σ with increasing concentration. Some models display a plateau for low concentrations, while some display a monotonic decay. Since our focus is on the strictly monotonic relationship between σ and c, models with plateaus are not considered. The model chosen here reads

σ = σ 0 (1 - α 2σ 0 c 2 d ), (4.2)
which is non-dimensionalized as

σ = 1 - αC 2 0 2σ 0 c 2 . (4.3)
Here C 0 represents the average mass concentration of pollutant at the surface. Note that β can be linked to a Marangoni number, based on C 0 and the diffusion D s such that M a = (αC 2 0 R)/(D S µ), and to the Péclet number P e s = ΩR 2 /D s so that β = M a/P e S . In practice, the input parameter for the LSA is β.

Modelling of pollutant concentration

When all pollutants stay at the interface z = G, their concentration c(r, θ, t) obeys a superficial advection-diffusion equation of the form

∂c ∂t + ∇ s • (cu s ) + c(∇ s • n)(u • n) = 1 P e s ∆ s c, (4.7)
where ∇ s represents the gradient operator in the directions tangent to the interface, ∇ s represents the gradient operator in the directions tangent to the interface and ∆ s is the associated Laplacian [START_REF] Stone | A simple derivation of the time-dependent convective-diffusion equation for surfactant transport along a deforming surface[END_REF]). In Eq. (4.7), the original velocity field u is split into a normal component (u • n)n and the resulting tangential component

u s = u -(u • n)n.
In the simple case where n = e z Eq. (4.7) reduces to

∂c ∂t + ∇ s • (cu s ) = 1 P e s ∆ s c. (4.8)
We consider a decomposition into base flow and perturbation, where the perturbation is written using a complex ansatz of the form e λt+imθ , such that u r = U r + u * r e λt+imθ , u θ = U θ + iu * θ e λt+imθ , and c = C + c * e λt+imθ . For the steady axisymmetric base flow characterised by the velocity field U and the concentration field C, Eq. (4.8) becomes

1 r ∂rCU s r ∂r - 1 P e s ( ∂ 2 C ∂r 2 + 1 r ∂C ∂r ) = 0. (4.9)
By substracting Eq. (4.9) from Eq. (4.8), the equation for the fluctuation concentration (u, c) reads

-λc * = 1 r ( ∂rCu * r ∂r + ∂rU r c * ∂r ) + imU θ c * r - mCu * θ r - 1 P e s ( ∂ 2 c * ∂r 2 + 1 r ∂c * ∂r - m 2 c * r 2
). (4.10) The diffusion coefficient for pollutant D s is usually one or two order of magnitude smaller than that of kinematic viscosity and thus in the present case, superficial diffusion effects remain small with respect to advection effects.

The constraint (4.4) in non-dimensional form reads :

1 π 2π 0 1 0 c(r, θ, t)rdrdθ = 1, (4.11)
and reduces for the steady axisymmetric base flow to For the base flow, axisymmetry implies ∂C ∂r = 0 at the axis. The constraint 4.12 also imposes a zero mass flux at r = 1. For the perturbation field c * , the boundary conditions depend on the value of the azimuthal wavenumber m [START_REF] Kahouadji | Thermocapillary instabilities in a laterally heated liquid bridge with end wall rotation[END_REF]. For m 1 (the case of interest), c * = 0 is imposed at the axis and ∂c * ∂r = 0 at the outer wall. All superscripts * are from here on dropped for simplicity.

Structure of the modified base flow

As demonstrated in [START_REF] Lopez | Coupling between a viscoelastic gas/liquid interface and swirling vortex flow[END_REF], the presence of a surfactant layer at the interface modifies the structure of the base flow. However the potential influence on its linear stability has not been investigated yet. In this subsection, we study the influence of the pollution concentration β, modelled using the surfactant law (4.2), on the base flow for G=1/14 and Re=18620. Increasing β causes a small but monotonous decrease of the length of the meridional recirculation, see figure 8. This is accompanied by the progressive disappearance of the overshoot in U θ , evident in fig. 9. This observation is directly consistent with the experimental measurements, in which no overshoot has been found for z=0.8G. 1.

thresholds Re c for m=4 and 5, with G fixed to 1/14 except when noted. The numerical resolution is unchanged compared to the pollution-free case. The Péclet number, though in principle larger, is hence limited to P e = 10 3 in order to prevent steeper gradients and numerical issues. The focus on m=4 and 5 mirrors the modal selection predicted for the reference case β=0, and is also consistent with experimental findings at onset.

The neutral curves Re c (β) obtained as β is varied are shown in figure 10, where neutral modes m=4 and 5 appear as squares and five-pointed stars, respectively. The most striking result in fig. 10 is the dramatic drop in Re c occurring at β ≈ 0.48 for m=4 and β ≈ 0.14 for m=5. The asymptotic value of Re c predicted for large β and approached for β as small as 0.2, is also below 3000 in much closer agreement with experimental estimations (shown as the blue strip in fig. 10) than the numerical prediction with β=0. These results suggest that a minute amount of surfactants can dramatically impact the flow stability, while additional pollution does not worsen the phenomenon further. In other words the effect of pollution is almost binary : either the interface is perfectly clean and the stability of the flow obeys the classical prediction from Section 3, or it is not and the stability characteristics of the flow are of a fully different kind. This scenario is so far consistent with the experimental reproducibility of Re c .

For the mode m=5, a sharp change of slope is evident for Re c =17259, β = 0.14. This marks the presence of a codimension two point where two different marginal curves for two different modes m=5 intersect in a double Hopf bifurcation : on each side of the corresponding value of β, these are not the same family of eigenmodes that go unstable first, despite a common azimuthal wavenumber m. The crossing of eigenvalues is confirmed in fig. 11 where the pair of eigenvalues is displayed on each side of the crossing. In each case, the branch taking over for larger β does apparently not extend down to β=0 : it corresponds to a new instability not found in the clean interface case. The evolution of the two leading eigenvalues as functions of both β and Re is detailed in fig. 12(a). The least stable eigenvalue for the 'clean' case β=0 is labelled 'C m ', where m=5, whereas the least stable eigenvalues for the 'dirty' case β 1 is simply labelled 'D m '. From fig. 12(a) it appears that the trajectories of the eigenvalues C 5 and D 5 in the complex plane, for variations of Re and β, obey different trajectories : C 5 gets destabilised by increasing Re but stabilised by increasing β, whereas D 5 is destabilised by both increasing β and increasing Re.

In the case m = 4, an equivalent codimension two point can be identified in figures 10 and 11, at Re c = 18869, β = 0.17. For this mode the drop in Re c is more dramatic than for m=5, and does not occur immediately after the codimension two point. Instead, the new branch (in green in fig. 10) continues to increase until at Re c = 20951, β = 0.48 where it turns back. Again the trajectory of the corresponding eigenvalues C 4 and D 4 is documented fig. 12(b). The trajectories of C 4 in the complex plane are similar to those of C 5 . However the scenario for D 4 differs from that for D 5 : an increase in Re stabilises the corresponding eigenmode whereas an increase in β destabilises it. Interestingly, the asymptotic value of Re c for m = 4 as well as the corresponding value of β at which the lowest values of Re c (β) are reached, seem to match that for m=5. This suggests that the value of Re c does not, for large β, depend on the value of m, at least for this value of G. Preliminary computations for G=1/4 have not confirmed this observation. A parametric study of Re c as a function of both G and m would shed light on this question, but this lies outside the present scope.

The evolution of the angular frequency of the pattern at Re = Re c is displayed in figure 13 for both m=4 and 5 as β is varied. Direct comparison of this figure with fig. 10 shows that each jump to a new branch corresponds to a discontinuity in angular frequency. Again the quantitative match with the experimental angular frequency is much more satisfying at finite β than around β=0. For large β, the data approaches the blue strip in fig. 13 within 1% or less.

Visualisations in physical space of the different modal families for a common wavenumber m are again consistent with experimental findings. The superficial axial for G=1/14, obtained using LSA for the modes m=4 (squares) and m=5 (stars). The blue strip corresponds to the experimental value obtained using LDV (cf Table 1).

vorticity ω z (r, θ, z = G) for the first marginal eigenmode m=5 is displayed for two values of β=0 and 5 in figures 14a and 14b. From figs. 10 and 13 it is now clear that these two eigenmodes correspond to two different modal families. For the β=0 case the marginal eigenmode corresponds to the one found in [START_REF] Kahouadji | Thermocapillary instabilities in a laterally heated liquid bridge with end wall rotation[END_REF], it features sharp gradients close to r ≈ 0.67R and thin counter-vorticity tongues on each side of this layer of strong shear. The structure of the marginal eigenmode for β=5, in turn, does not feature such sharp gradients, and displays much more regular cellular structures in a zone closer to the outer wall (see fig. 14b). This is consistent with the different All the results above support that the numerical prediction for finite β (nonclean interface) is consistent, both regarding the base flow and its marginally unstable eigenmodes, with the experimental findings, whereas the clean interface (β=0) hypothesis is not.

Frozen interface condition

Search for a simpler parameter-free interface condition

The results from the previous section have shown that a simple surface pollution model can capture qualitatively and quantitatively well the main features of the instability under investigation without any of the physical processes related to the surface contamination. While it is possible to make the model quantitatively closer to the real case by adding more parameters, we search in this section for an even simpler model for the interfacial conditions. In particular, we would ideally like to have an analytically simple boundary condition for the velocity at the liquid/gas interface that is parameterfree and does not request simulating additional concentration fields. This would make the implementation of such a model easy to achieve in practice in existing numerical codes, without depending on the precise (and usually unknown) details on the adsorption at the interface. The results of Section 4 suggest the presence of a well-established asymptotic regime for large β, and the synthetic interfacial condition sought for is requested to match the large β limit. Several authors have already reported that the presence of pollutants in two-dimensional geometries removes slip at the interface, in condtradiction with the traditional hypothesis of free slip. In particular, whereas in our flow case the larger azimuthal velocity remains weakly affected by pollutants, the radial component of the velocity is severely diminished, making the hypothesis of vanishing u r at the interface plausible [START_REF] Spohn | Recirculating flows in a cylindrical tank[END_REF][START_REF] Lopez | Surfactant-influenced gas-liquid interfaces: Nonlinear equation of state and finite surface viscosities[END_REF]. This is achieved in the numerical codes by changing the free-slip boundary conditions at the interface from (2.5) into :

u r = 0, ∂u θ ∂z = 0, u z = 0 at z = G. (5.1)
The prime advantage of such an interfacial condition is its simplicity : as requested it is parameter-free, chemistry-free, it does not request coupling with an equation for the concentration, and it does not rely on any closure for the effective surface tension. In the following, we assess numerically whether imposing this interfacial condition U r = 0 for the base flow is a satisfying hypothesis.

Base Flow

We estimate first how much the 'frozen' condition (5.1) is consistent with large values of β by assessing the spatial structure of the base flow. Figure 15 velocity components, respectively U r (r) and U θ (r) evaluated at z = G. The length of the radial interval where the radial velocity U r is non-zero decreases with increasing β, and the minimum value of U r also approaches zero, suggesting absolute convergence to a homogeneous U r =0 profile. This justifies, for β large enough, U r = 0 as an interfacial boundary condition, in agreement with previous experimental observations. The major difference in the concentration curves (15c) is the non-zero concentration of pollutants for every radial position when β 5. For smaller values of β, the concentration is advected towards the axis by the meridional recirculation, allowing for a small clean area to remain close to the outer wall. This leads to completely different vorticity profiles at the surface : for smaller values of β the vorticity is gathered around the radial position where the drop of concentration occurs. However for β ≈ 5 and above, the vorticity is stretched over a larger radial range, and converges to the frozen interface case. Interestingly, despite the fact that U r is not exactly zero at the interface, all base flow profiles between β = 5 and the frozen surface condition appear identical, as is visually clear from Figs 16 and 17. In practice β is limited by chemistry considerations : high values of β correspond to a highly polluted surface, and we may, under these conditions, exceed the critical micellar concentration (CMC) above which the validity of the superficial model breaks down entirely.

Nonlinear dynamics

The condition (5.1) is here explicitly imposed in the nonlinear DNS calculations too. The simulation has been conducted with the same spatial resolution as for in Section 2. The newly computed instability pattern is shown in figure 18 past the initial transient. Without imposing any rotational symmetry, we see that above Re c the most unstable mode emerges with an azimuthal wavenumber m=5. Compared with simulations based 7b), the present nonlinear regime is much more predictable : this m=5 mode persists for the whole observation time of up to t=1400 time units) and the frequency spectrum remains limited to multiples of the fundamental frequency. The comparison between the numerical and experimental pointwise spectra is displayed in fig. 19 for the same value of Re, and deserves to be compared with fig. 7. From such a cross-comparison, it is a non-ambiguous fact here that the frozen condition leads to a much better spectral reproduction of the experimental flow. Besides, although such comments are subjective, we report that the aspect of the pattern in fig. 18 From figure 3d we can estimate the width of this annular stripe in the experiment to [0.62-0.92]. This rapid quantitative comparison confirms that the instability pattern with a frozen surface is much closer to the experimental pattern than the pattern from the free surface simulation.

Critical Reynolds number and least stable mode

We observe in figure 10 that Re c (β) hardly evolves once β is large enough (larger than e.g. 5). This behaviour is confirmed for larger β, where the increase in Re c remains limited. The instability pattern appears similar to those shown in figure 14, which suggests that the most unstable mode of the frozen condition again belongs to an unstable branch different from the 'clean' unstable mode identified with the free slip Table 4: Re c for m=5 and four sets of boundary conditions for the couple (U ,u * ).

Conclusions on the frozen surface condition

The nonlinear dynamics captured in DNS using the frozen surface condition is in excellent qualitative match with experimental measurements both from the point of view of the dynamics and the modal content of the saturated flow. This is again confirmed by the good agreement between the amplitude spectra shown in 19. The base flow with the frozen surface condition is hardly distinguishable from the base flow obtained using the pollutant model for β 5. Nevertheless, once again the comparison of the values of Re c is not favourable, as for the frozen surface Re c is 260% (10555 vs 2934) higher than for β=5, which discredits the frozen condition as a direct substitute to the free-slip condition. This negative conclusion is further confirmed for G=1/4 (the value considered by [START_REF] Lopez | Symmetry breaking in free-surface cylinder flows[END_REF]), where the most unstable mode remains poorly representative of experimental visualisations and Re c is pushed even further up. These quick tests reveal how sensitive the instability threshold is to the choice of boundary conditions.

The quantitative discrepancy in Re c can eventually be resolved by introducing a new boundary condition of a mixed type. Whereas the hypothesis of vanishing U r is very satisfying for the base flow and large enough values of β, the above argument suggests that in terms of time-dependent perturbations, vanishing u r is inconsistent with the temporal oscillations of the concentration field. We suggest to separate the velocity field u into its base flow component U and its perturbations u * and to apply a different set of boundary conditions to U and u * . The situation is summed up in Table 4, where the four possible combinaisons of boundary conditions are considered and Re c has been re-computed using ROSE for each case. Consistently with the previous arguments, we focus on the mixed type boundary condition where the base flow obeys a frozen condition whereas the perturbation obeys the free slip condition. The value of Re c is now 2776, much lower than the fully frozen threshold value of 10555. The quantitative mismatch is now reduced down to 6% (2779 vs 2934) and 31% (2779 vs 4230) compared to β=5 and experimental Re c , respectively, while the spatial structure (see figs. 14 and 20 for the eigenmodes) is compatible with the frozen case.

Discussions and perspectives

The hydrodynamic instability occuring inside a fixed, cylindrical cavity with a rotating bottom has been investigated for a small form factor depth/radius G=1/14. The selection of a least stable mode with azimuthal wavenumber m=5, predicted by linear instability analysis, is verified experimentally, as well as numerically using DNS assuming no stress at the liquid interface. Using water -the most widespread liquid -as the experimental fluid, a robust quantitative and qualitative mismatch is however evidenced between our experiments and our numerics. The mismatch concerns the presence of an overshoot in the azimuthal velocity profile and, crucially, Re c for the development of the instability is overestimated in the numerics by a factor of more than 4. Regarding the mismatch between our experimental estimation of Re c and the literature, the use of Kalliroscope as a marker, used precisely for the identification of instability thresholds, is our best suspect to explain the discrepancy. The results obtained here using LDV, quantitatively safer, demonstrate that thresholds formerly deduced from visualisation using markers, are over-evaluated. The linear and nonlinear numerical approaches, on their side, report a robust threshold Re c , that still differs strongly from the experimental one. After a cautious search for possible experimental flaws, the standard free slip interface condition (σ=cst, h=cst, free slip interface) used in the simulations emerged as the most credible source of mismatch, in line with former investigations by [START_REF] Spohn | Recirculating flows in a cylindrical tank[END_REF] and Hirsa et al. (2002a) : in experiments, such an ideal interfacial condition cannot be matched due to residual ambient air pollution. The inevitable presence of pollutants at the interface modifies the surface tension of the flow and, as a consequence, impacts the velocity field of the base flow and shifts the instability threshold. A pollution model has been implemented into the linear stability solver, based on a modification of the effective surface tension by the presence of a superficial concentration of unknown pollutants. Using a quadratic closure between the tension surface and the surperficial pollutant concentration inspired by surfactant studies, yields results quantitatively consistent with experiments : for sufficiently large value of β (the parameter that pilots the surface contamination), Re c drops by more than 80% and the mismatch on Re c goes down roughly from 400% down to 30%. Interestingly, the instability mode selected for finite pollutant concentrations does not belong to the same modal families as predicted by linear stability theory in the clean interface case : new branches of "dirty" modes destabilize for small yet finite concentration levels and take over as least stable modes. The corresponding eigenmodes are very stable for clean interface conditions and have not been identified before. In terms of bifurcations, the robust mismatch in Re c , angular frequency and flow structure between numerics and experiments can hence be explained, at least for the case of the modes m=4 and 5 investigated here, as a jump from one modal family to another one as the contamination of the surface increases. The spatial structure of the base flow is also more consistent with LDV measurements : the meridional recirculation length is reduced and the overshoot in azimuthal velocity vanishes. For β 5, the flow at the interface verifies an approximate no-slip condition for the radial velocity component. As a consequence, in an effort to deliver a simpler parameter-free model boundary condition for unclean liquid/gas interfaces, the "frozen condition" u r =0 was also simulated. Whereas it displays better qualitative agreement as well as simpler nonlinear dynamics consistent now with experiments, at least for low G, the threshold value Re c remains too high compared to experiments. This new quantitative mismatch is eventually resolve once and for all by considering an interfacial boundary condition of mixed type : frozen for the base flow and free-slip for the perturbation.

Interfacial experiments involving water have long had the reputation of being 'difficult' in the sense that Marangoni effects linked with variations of the surface tension are hard to tame. The present hypothesis of a modification of the surface tension by pollution effects is one such illustration. The surprising effect of this pollution is, despite relatively small modifications of the structure of the base flow, an important quantitative impact on the stability thresholds. Besides, not only does the instability mode change its growth rate, it also belongs to another family of destabilised modes compared to the clean interface case. From such a simple conclusion it is tempting to critically revisit the discrepancies between experimental studies and to deduce that higher values of Re c (as reported in [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF]) are linked to a cleaner interface due to different experimental conditions. While this is a priori possible (and very difficult to assess rigorously), it does not remove the caveat that Kalliroscope visualisations are poorly reliable in terms of measurements. Besides an interaction of the marker itself with the solvent cannot be excluded for high Kalliroscope concentrations. We are hence not in a position to conclude about the values reported by various teams using Kalliroscope or other markers, and encourage instead the use of non-intrusive techniques such as LDV for more reliable estimations. Other experimental improvements could here be useful, such as measuring simultaneously several velocity components, including smaller components such as the axial one and the radial component near the interface. This could allow for a critical evaluation of the model interfacial condition suggested in Section 5.

While the present study is essentially a proof of concept that the stability characteristics of a given flow case depend heavily on the surface pollution, the simplicity of the analytical model in Section 4 must be kept in mind. The advantage of such a simple model is a straightforward identification of the mechanisms altering the spatial structure of the base flow. The main difficulty lies in the mathematical parametrisation of a chemically complex phenomenon. Adsorption of pollutants by the interface is an unsteady process that depends on the precise chemical composition of the ambient particles in the air and of the exact properties of the liquid. None of these hypotheses have been included in the present model, resulting in a simple law parametrised by one unique real parameter β and unable to cope with different chemical compositions. Besides, as is made clear in Hirsa et al. (2002a), the rheology of the interface does not even need to be considered Newtonian, which implies the arduous determination of further rheology constants. Further chemical complications can occur for increased concentration levels, notably bulk diffusion (necessitating a volumetric model rather than a superficial one) and later the formation of micelles inside the bulk of the fluid. Ideally, the modelling of pollution effects should be compared with an experimental set-up where the surperficial concentration of each pollutant can be quantitatively controlled and properly modelled. It is not excluded that each different pollutants contribute differently to the final surface tension rather than all obeying the quadratic law of Eq. (4.2).

Finally, the present set-up is still academically simple in the sense that no deformation of the interface needs to be considered at such low rotation rates. While this is technically much more involved, especially on the numerical side (see e.g. [START_REF] Yang | Large axisymmetric surface deformation and dewetting in the flow above a rotating disk in a cylindrical tank: spin-up and permanent regimes[END_REF]), it is suspected that pollution effects can also affect the thresholds in the large deformation regime investigated by many others [START_REF] Vatistas | Experimental confirmation of Kelvin's equilibria[END_REF][START_REF] Tophøj | Rotating polygon instability of a swirling free surface flow[END_REF]). These regime involve not only finite deformations of the fluid interface but also partial dewetting which makes the dynamics of the concentration field more complex by involving moving triple contact lines. Eventually, it will be interesting to see how the trend evidenced in the present study (the decrease of Re c by ambient pollution despite a calmer nonlinear regime) can be extended to other unstable flow configurations. This work was supported by the French Agence Nationale de la Recherche under the ANR ETAE Project No. ANR-16-CE08-0011. HPC resources from GENCI-IDRIS (Grant No. 2017-2a10308) are also acknowledged. Declaration of Interests : The authors report no conflict of interest. The authors thank Ivan Delbende, François Gallaire, Wietze Herreman, François Lusseyran, Wen Yang for fruitful discussions, the technical team CTEMO at LIMSI, as well as Fréderic Moisy and Marc Rabaud for their help on the experimental facility.

Figure 1 :

 1 Figure 1: Sketch of the axisymmetric base flow for small aspect ratio G = H/R

Figure 2 :

 2 Figure 2: Instability patterns breaking the axisymmetry of the flow. Photograph taken from above (ink visualisation) in our experimental set-up. From left to right : modes m = 3, 4 and 5 obtained for different aspect ratios and different values of Re above the effective Re c .

Figure 5 :

 5 Figure5: Velocity profiles of u θ (r) below the free surface (z = 0.8G) for G=1/14 and Re = 18620. Comparison between the base flow, spatially averaged DNS and the temporal average for LDV (experiment). LDV acquisition timespan is much larger than the instability period.

Figure 6 :

 6 Figure 6: DNS axial vorticity fluctuation for the free surface condition, Re=17100 (sligthly above Re c (6a), and Re=18620 (6b). With the increase of Re, the mode m=5 selected at Re c evolves into a modulated m=7 pattern.

Figure 7 :

 7 Figure 7: Comparison of frequency amplitude spectra for u θ (t) measured at r = 0.76, z = 0.8G, for the saturated m=5 regime at Re = 18620. Left: experimental, maximum peak at f = 4.15. Right: DNS with free surface condition, same parameters. The frequencies are not normalized by the azimuthal wavenumber.

c

  is non-dimensionalized as c = c d /C 0 . α is a constant coming from the Taylor expansion around c=0 of the model in[START_REF] Kwan | A mathematical and numerical study of incompressible flows with a surfactant monolayer[END_REF]. Eqs 4.1 can hence be rewritten as :

4. 4 .Figure 8 :Figure 9 :

 489 Figure 8: Evolution of the stream function ψ for the base flow with increasing concentration β, Re = 18620 and G=1/14. The same contour values are chosen for all cases. Negatives and positives contours use different scales to highlight the weak recirculation bubble. Negatives contour values (dashed): 4 equispaced levels in [ψ minψ min /5]. Positives contour values : (solid lines): 9 values equispaced levels in [ψ max /10ψ max ]. Zero contour level (solid black lines). ψ min = -8.1.10 -5 and ψ max = 2.3.10 -3 .

Figure 10 :

 10 Figure 10: Neutral curve Re c (β) estimated from LSA for G=1/14 and for the modes m=4 (squares) and m=5 (stars). The blue strip corresponds to the experimental value of Re c , independent of the model based on β. The thickness of the strip is based on lower and upper bounds for Re c from Table1.

Figure 11 :

 11 Figure 11: Zoom on codimension two points of Fig. 10 in the (β,Re c ) plane, for G=1/14, m=4 and m=5.

Figure 12 :Figure 13 :

 1213 Figure 12: Evolution of the eigenvalues associated with the most dangerous modes ('clean' mode C and 'dirty' mode D) as a function of Re and β for G = 1/14. (a) m = 5 and (b) m = 4. The two evolutions start at Re = 15 000 and β = 0. The C and D eigenvalues follow the lines indicated by the red arrows as Re is increased from 15 000 to 17 200 while β = 0. The eigenvalues evolve with varying β (at constant Re=15000, blue arrows) with equally spaced steps from 0 to 0.18 for m =5, and from 0 to 0.5 for m = 4.

Figure 14 :

 14 Figure 14: Vorticity ω z (r, θ) at the fluid interface (normalised by its maximum) for the least stable eigenmode m=5 at Re=18620. a) β=0, b) β=5.

Figure 15 :

 15 Figure 15: Comparison of several fields of the base flow at the surface, for few values of β, at Re=18620. 15a) radial velocity U r , 15b) azimuthal velocity U θ , 15c) pollutant concentration C, 15d) azimuthal vorticity ω s . The results for the frozen surface condition are also included in the figures.

  contains various radial profiles at the interface for the base flow, as β is increased beyond the values shown previously. Values of β up to 10 2 or 10 4 have been considered in order to monitor the dependence of the base flow on β. Fig.15aand 15b display the radial and the azimuthal

Figure 16 :

 16 Figure 16: Comparison of ψ, between β=5 and the frozen surface condition, at Re = 18620 (computed with ROSE). Iso-contour levels as in fig. 8.

Figure 17 :

 17 Figure 17: Meridian (r, z) cut for β=5 and the frozen surface condition for Re = 18620(computed with ROSE). Iso-contour levels as in fig. 9.

Figure 18 :

 18 Figure 18: Axial vorticity of the fluctuations at the interface, DNS with frozen condition u r = 0 at the interface for Re = 18620.

  is visually closer to the experimental one than those of figure 6b for the free surface condition. More quantitatively, the approximate radius range where the axial vorticity fluctuations are concentrated in figure 18 is [0.59-0.91], while in figure 6a they are limited to [0.64-0.72].

Figure 19 :

 19 Figure 19: Comparison of amplitude spectra of u θ measured at (r=0.76, z=0.8G) :experimental data (blue) versus DNS for the frozen surface condition (orange), Re=18620.

Figure 20 :

 20 Figure 20: Vorticity ω z (r, θ) at the fluid interface (normalised by its maximum) for the least stable eigenmode m=5 at Re=18620. a) frozen surface, b) mixed condition. The mixed condition corresponds to a frozen surface condition for the base flow, and free surface condition for the perturbation field.

Table 2 :

 2 Critical Reynolds number for G=1/4. Comparison between experiments by Lopez et al. (2004), present experiments and LSA. The percentage of glycerol indicated is a weight percentage.

Table 3 :

 3 Critical Reynolds number estimated by LSA for G = 1/14.

	on free slip (see fig.

† Email address for correspondence: laurent.martin_witkowski@sorbonne-universite.fr