Giulio Facchini
email: facchini@irphe.univ-mrs.fr

Benjamin Favier

Patrice Le Gal

Meng Wang

Michael Le Bars

The linear instability of the stratified plane Couette flow

We present the stability analysis of a plane Couette flow which is stably stratified in the vertical direction orthogonally to the horizontal shear. Interest in such a flow comes from geophysical and astrophysical applications where background shear and vertical stable stratification commonly coexist. We perform the linear stability analysis of the flow in a domain which is periodic in the stream-wise and vertical directions and confined in the cross-stream direction. The stability diagram is constructed as a function of the Reynolds number Re and the Froude number F r, which compares the importance of shear and stratification. We find that the flow becomes unstable when shear and stratification are of the same order (i.e. F r ∼ 1) and above a moderate value of the Reynolds number Re 700. The instability results from a resonance mechanism already known in the context of channel flows, for instance the unstratified plane Couette flow in the shallow water approximation. The result is confirmed by fully non linear direct numerical simulations and to the best of our knowledge, constitutes the first evidence of linear instability in a vertically stratified plane Couette flow. We also report the study of a laboratory flow generated by a transparent belt entrained by two vertical cylinders and immersed in a tank filled with salty water linearly stratified in density. We observe the emergence of a robust spatio-temporal pattern close to the threshold values of F r and Re indicated by linear analysis, and explore the accessible part of the stability diagram. With the support of numerical simulations we conclude that the observed pattern is a signature of the same instability predicted by the linear theory, although slightly modified due to streamwise confinement.

Introduction

Shear and density stratification are ubiquitous features of flows on Earth and can strongly affect the dynamic of different fluids like air in the atmosphere or water in the ocean. More generally the interest for the stability of parallel flows dates back to the second half of the nineteenth century [START_REF] Helmholtz | XLIII. on discontinuous movements of fluids[END_REF][START_REF] Kelvin | XLVI. hydrokinetic solutions and observations[END_REF]) and the first crucial statement came with [START_REF] Rayleigh | On the stability, or instability, of certain fluid motions[END_REF] who gave his name to the famous inflexion point theorem proving a necessary criterion for an inviscid homogeneous parallel flow to be unstable. Contemporarily the first laboratory experiments performed by [START_REF] Reynolds | An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels[END_REF] showed that also inflexion-free flows can run unstable at sufficiently high Re, thus highlighting the need for a viscous analysis. Still more than a century ago [START_REF] Orr | The stability or instability of the steady motions of a perfect liquid and of a viscous liquid. part I: A perfect liquid[END_REF] provided a viscous equivalent of the Rayleigh principle. Nonetheless, as reviewed by [START_REF] Bayly | Instability mechanisms in shear-flow transition[END_REF], providing a solution of the Orr-Sommerfeld equation at large Re number turns out to be exceedingly difficult and has drawn since then, the attention of many studies [START_REF] Heisenberg | Über Stabilität und Turbulenz von Flüssigkeitsströmen[END_REF][START_REF] Schlichting | Zur enstehung der turbulenz bei der plattenstrmung[END_REF][START_REF] Lin | The theory of hydrodynamic stability[END_REF]. Interestingly even for the simplest profile of parallel flow, a conclusive answer as been lacking for almost a century as reported by [START_REF] Davey | On the stability of plane Couette flow to infinitesimal disturbances[END_REF]: 'It has been conjectured for many years that plane Couette flow is stable to infinitesimal disturbances although this has never been proved [...] We obtain new evidence that the conjecture is, in all probability, correct'. Since then the stability analysis of the plane Couette (PC hereafter) flow continues to be of deep interest in studying the transition to turbulence via non-linear mechanisms [START_REF] Barkley | Computational study of turbulent laminar patterns in Couette flow[END_REF] but its linear stability is nowadays no more questioned [START_REF] Romanov | Stability of plane-parallel Couette flow[END_REF]). In the present work we show that by adding a vertical linear (stable) density stratification, the PC flow becomes unstable, at strikingly moderate Re numbers, typically Re 700. The observed instability relies on the same resonance mechanism showed by [START_REF] Satomura | An investigation of shear instability in a shallow water[END_REF] for shallow water waves, here extended to the case of internal gravity waves. An interesting feature of this finding is that density stratification is generally thought to be stabilising as it inhibits vertical motion. Nonetheless our counterintuitive result does not come as a prime novelty. In the close context of rotating-stratified (and sheared) flows [START_REF] Molemaker | Instability and equilibration of centrifugally stable stratified taylor-couette flow[END_REF] and [START_REF] Yavneh | Non-axisymmetric instability of centrifugally stable stratified Taylor-Couette flow[END_REF] questioned the other Rayleigh celebrated criterion [START_REF] Rayleigh | On the dynamics of revolving fluids[END_REF]) and showed that Rayleigh-stable Taylor Couette flows may become unstable when adding linear density stratification. The Strato-Rotational-Instability as successively named by [START_REF] Dubrulle | An hydrodynamic shear instability in stratified disks[END_REF] was observed in the laboratory a few years later (Le [START_REF] Bars | Experimental analysis of the stratorotational instability in a cylindrical couette flow[END_REF] and is still the subject of experiments [START_REF] Ibanez | Observations of the stratorotational instability in rotating concentric cylinders[END_REF]. The stability analysis of parallel flows where shear coexists with stratification has also a long tradition. The most famous shear instability, i.e. the Kelvin-Helmoltz instability was found indeed in the context of a two layers fluid endowed with different velocity and density [START_REF] Helmholtz | XLIII. on discontinuous movements of fluids[END_REF][START_REF] Kelvin | XLVI. hydrokinetic solutions and observations[END_REF]. This work was extended to the three density layers configuration, with constant shear in the middle one, by [START_REF] Taylor | Effect of variation in density on the stability of superposed streams of fluid[END_REF] and [START_REF] Holmboe | On the behaviour of symmetric waves in stratified shear layers[END_REF], who identified two different instability mechanisms, and later by [START_REF] Caulfield | Multiple linear instability of layered stratified shear flow[END_REF] who isolated a third possibility. [START_REF] Miles | On the stability of heterogeneous shear flows[END_REF] and [START_REF] Howard | Note on a paper of john w. miles[END_REF] gave the stability criterion of the Kelvin-Helmoltz instability, for the case of continuous linear stratification. Since then, most of the studies have focused on the configuration where density gradient and shear are parallel. On the contrary only a few (e.g. [START_REF] Deloncle | Three-dimensional stability of a horizontally sheared flow in a stably stratified fluid[END_REF][START_REF] Candelier | Shear instability in a stratified fluid when shear and stratification are not aligned[END_REF][START_REF] Arratia | Non-modal instability mechanisms in stratified and homogeneous shear flow[END_REF]) recently considered the case of non alignment as reviewed by [START_REF]Stabilité d'un écoulement stratifié sur une paroi et dans un canal[END_REF], who also showed [START_REF] Chen | Instability of a boundary layer flow on a vertical wall in a stably stratified fluid[END_REF]) that a free-inflexion boundary layer profile is linearly unstable when linear stratification is added. Restricting ourselves to the case of a vertically stratified and horizontally sheared plane Couette flow, [START_REF] Bakas | Gravity waves in a horizontal shear flow. part ii: Interaction between gravity waves and potential vorticity perturbations[END_REF] already considered the problem of the linear stability while investigating the interaction between gravity waves and potential vorticity perturbations, but curiously did not explore the linearly unstable region. We repeat their linear stability analysis using a pseudo-spectral method (i.e. with the same approach as [START_REF] Chen | Instability of a boundary layer flow on a vertical wall in a stably stratified fluid[END_REF]) and find that exponentially growing modes appear at moderate Re number Re ∼ 700 and F r 1, for non vanishing vertical and horizontal wave number k x /k z ∼ 0.2. Results are confirmed by fully non-linear direct numerical simulations (DNS). We also analyse the laboratory flow produced by a shearing device immersed in a rectangular tank filled with salty water linearly stratified in density. We verify that a fairly parallel PC flow can be generated and observe that beyond a moderate Re number Re 1000 and for F r number close to 1, a robust velocity pattern appears in the vertical mid-plane parallel to the shear, that is where no motion is expected for a stable PC flow. In particular perturbations grow in an exponential manner and looking at how their saturation amplitude varies in the (Re, F r) space, we find that an abrupt transition is present close to the marginal stability limit predicted by linear stability analysis. The quantitative agreement of the observed spatio-temporal pattern with the linear theory is only partial, which we claim to be a consequence of the finite streamwise size of our device. This hypothesis is largely discussed and supported by the results of additional DNS confirming that the finite size of the domain weakly affects the base flow, but does modify the shape of the perturbation pattern. We conclude that the observed instability indeed corresponds to the linear instability of the vertically stratified PC flow modified by finite size effects and that a redesigned experiment may reproduce more faithfully the spatio-temporal pattern predicted by the linear theory.

The paper is organized as follows. In section 2 we define the observed flow with its governing equations and describe the linear stability approach. In section 3 we report the results of linear analysis and in section 4 those of direct numerical simulations. The experiments are described in section 5 and the experimental results compared with the linear theory and direct numerical simulations in section 6. In section 7 we summarise our study and briefly discuss possible applications and future development of the present work.

Theoretical frame

We consider the plane Couette flow generated by two parallel walls moving at opposite velocity for a fluid which is stably stratified in density as sketched in figure 1. We denote x the stream-wise direction, ŷ the cross-stream direction (i.e. the direction of the shear) and ẑ the vertical direction (i.e. the direction of the stratification). The vector g denotes gravity while red arrows sketch the shape of the constant shear profile U (y) and red shading mimic vertical stratification ρ(z). In the Boussinesq approximation we obtain the following system of equations:

∂u ∂t + (u • ∇)u = - ∇p ρ 0 - ρ ρ 0 g ẑ + ν∇ 2 u, (2.1)
∇ • u = 0, (2.2) ∂ρ ∂t + (u • ∇)ρ - N 2 g ρ 0 w ẑ = k∇ 2 ρ , (2.3)
where we decompose the pressure and density fields p and ρ in a perturbation p and ρ and a stationary part p = p 0 + ρ 0 gz -N 2 z 2 ρ 0 /2 and ρ = ρ 0 (1 -N 2 z/g), with p 0 and ρ 0 two constant reference values. We indicate with N = -∂ z ρ(g/ρ 0 ) the background Brunt-Väisälä frequency, while ν and k denote viscosity and salt diffusivity.

Linear stability analysis

We perform the linear stability analysis of the equations (2.1)-(2.3) in a Cartesian box of dimensions (L x , L y , L z ) centered in x = y = z = 0. To this aim we introduce the buoyancy b = ρ /ρ 0 g and decompose the velocity perturbation u in a perturbation u and a base solution U = -U 0 y x. Boundary conditions are periodic in the stream-wise and vertical directions and no-slip, i.e. u = 0 at the rigid walls y = ±L y /2. Buoyancy where we use again symbols u, p and b to simplify notations. Substituting in equations (2.1)-(2.3) and retaining only the first order terms we obtain:

-iωu = ik x uy + v -ik x p + 1 Re ∆ y u, (2.5) -iωv = ik x vy - dp dy + 1 Re ∆ y v, (2.6) -iωw = ik x wy - b F r 2 -ik z p + 1 Re ∆ y w, (2.7) 0 = ik x u + ik z w + dv dy , (2.8) -iωb = ik x by + w + 1 ReSc ∆ y b, (2.9)
where we denote with ∆ y the Laplacian operator

∆ y = d 2 /dy 2 -k 2 x -k 2 z .
The system of equations above is solved using a pseudo-spectral approach similarly to [START_REF] Chen | Instability of a boundary layer flow on a vertical wall in a stably stratified fluid[END_REF], the only difference is that discretisation is made on the Gauss-Lobatto collocation points of the Chebychev polynomials (i.e. instead of Laguerre) because this choice is well adapted to a two-side bounded domain. The generalized eigenvalue problem Af = ωBf for f = [u, v, w, b, p] is solved with the QZ algorithm. In parallel we also consider the inviscid approach which consists in neglecting both viscous dissipation and salt diffusion, thus reducing the system (2.5)-(2.9) to one equation for the pressure: (2.10) where ω * = ω + k x y. The equation above is analogous to that provided by [START_REF] Kushner | Coupled Kelvin-wave and mirage-wave instabilities in semigeostrophic dynamics[END_REF] who previously studied the stability of a vertically stratified PC flow in the presence of rotation f . In the limit of no rotation (f → 0) we verify that the two equations are the same but contrarily to [START_REF] Kushner | Coupled Kelvin-wave and mirage-wave instabilities in semigeostrophic dynamics[END_REF] we could not find a meaningful limit in which our equation (2.10) becomes autonomous in y. As a consequence we cannot provide a compact form for the dispersion relation. Nonetheless, looking at equation (2.10) is still extremely instructive. First, one observes that the second term in equation (2.10) possibly diverges at y = 0 when considering stationary modes which are marginally stable (e.g. ω = 0). This corresponds to the existence of a barotropic critical layer, which happens to be regularised because, from the symmetry of the base flow, we expect ∂p/∂y to be null in y = 0 for a stationary mode. Similarly the third term of (2.10) becomes critical in y * = ±1/k x F r when ω = 0. These are baroclinic critical layers, i.e. the locations where the Doppler shifted frequency ω * of internal waves matches the Brunt Väisälä frequency N. In different contexts critical layers can be excited and have been observed in experiments [START_REF] Boulanger | Tilt-induced instability of a stratified vortex[END_REF]) and numerical simulations [START_REF] Marcus | Three-dimensional vortices generated by self-replication in stably stratified rotating shear flows[END_REF]. However in our configuration the most unstable mode is always observed to be stationary and at wave numbers k max < 1/F r, which implies that the corresponding critical layers y * max are always situated outside the numerical domain, |y * max | > 1.

∂ 2 p ∂y 2 - 2k x ω * ∂p ∂y + k 2 z ω 2 * 1 -F r 2 ω 2 * -k 2 x p = 0,

Linear Stability results

We have already mentioned that for unstratified fluids (i.e. F r = ∞) the PC (unperturbed) profile is linearly stable for any value of the Reynolds number Re, thus we expect the flow to be potentially unstable only at finite values of the Froude number. The values of the Schmidt number for common salty water (i.e. in our experiments) is Sc ∼ 700 thus we preliminarily consider the limit Sc = ∞ and discuss the quality of this approximation at the end of this section.

As a first result we report that one stationary growing mode (i.e. Im(ω) > 0, Re(ω) = 0) appears at F r 1, wave numbers k x ∼ 0.8, k z ∼ 5 and remarkably moderate Reynolds number Re 700. In figure 2 we report the value of the imaginary part and the real part of the most unstable eigenmode for Re = 1000 and F r = 1 as a function of k x and k z . Looking at the imaginary part (left) one sees that the flow is unstable over a narrow elongated region centered in k x ∼ 0.8, k z ∼ 5 and stable elsewhere. Correspondingly the real part (center) is zero whenever the flow is unstable and non zero elsewhere. In figure 2 (right) we also report the values of the temporal frequency ω for all the eigenvalues and a various number of collocation points of N y = 129, 257 and 513 for the most unstable wave numbers (k x = 0.815, k z = 4.937). Physical eigenvalues correspond to the points where three different symbols are superposed, all other points corresponding to spurious numerical modes. The inset close to the origin of the diagram shows that a unique eigenvalue is present for which Im(ω c ) > 0. The value of ω c is stable to the variation of N y the number of collocation points, which indicates that the mode we observe is a physical one. We remark that this feature makes the systematic analysis of the space (Re, F r) particularly simple, for example differently from [START_REF] Chen | Instability of a boundary layer flow on a vertical wall in a stably stratified fluid[END_REF] the problem of non physical eigenvalues with positive Im(ω) is not present.

When increasing Re the unstable region in the (k x , k z ) space increases and unstable modes exist over a larger range of Froude number. On the top of figure 3 we report the value of Im(ω) and Re(C) for the most unstable mode at Re = 10000 and F r = 1, where we define C as C = ω/k x . In the same figure we report similar graphs for Im(ω) at Re = 10000, F r = 5 (left) and F r = 0.2 (right). Note that at these values of F r number there is no unstable mode at Re = 1000.

One sees that at F r = 1 (left panel) the diagram is now richer: besides the original unstable branch constituted by stationary modes (a), new unstable branches appear at larger k x which correspond to oscillatory (b, c, e) and stationary (d ) modes, as visible from the value of Re(ω)/k x ) (right). Note that within a same branch the value of Re(ω)/k x varies very little, while the value of Im(ω) shows a maximum and smoothly decreases to zero at the branch boundaries. The quantity Re(ω)/k x then characterises each different branch. A new oscillating branch (f ) also appears at smaller k x but it is still very weak and poorly visible at this Re number. On the bottom of figure 3 we report the value of Im(ω) in the space (k x , k z ) at F r number smaller and larger than one. When the Froude number is diminished to F r = 0.2 (left) we recover almost the same scenario, even if different branches look now more spaced one from another and appear at larger k x and k z similarly to other kinds of shear flows [START_REF] Deloncle | Three-dimensional stability of a horizontally sheared flow in a stably stratified fluid[END_REF][START_REF] Park | The stably stratified TaylorCouette flow is always unstable except for solid-body rotation[END_REF]. On the contrary when the Froude number is increased to F r = 5 (right) the unstable region is drastically reduced, as well as the growth rate, which is dropped by an order of magnitude. Also the most unstable mode moves toward lower value of k x while k z only slightly changes. As a general remark we observe that the unstable branches, i.e. the continuous regions defined by Im(ω) > 0, show an elongated shape. Precisely, unstable regions appear extended when moving along the curve k x k z = const while they are quite narrow in the orthogonal direction. We stress that this result is independent of the Froude number which suggests a self-similar behaviour as already observed by [START_REF] Deloncle | Three-dimensional stability of a horizontally sheared flow in a stably stratified fluid[END_REF]. Also unstable modes always appear at k z , k x = 0, i.e. the flow is linearly unstable only to three-dimensional perturbations, which is different from the studies of [START_REF] Deloncle | Three-dimensional stability of a horizontally sheared flow in a stably stratified fluid[END_REF] and Lucas et al. (2017), performed on different vertically stratified and horizontally sheared flow (the hyperbolic tangent shear profile and Kolmogorov flow respectively).

Stability Diagram

We explore the (Re, F r) parameter space over two decades around F r = 1 and for Re from 500 to 50000. For each combination (Re, F r), we solve the system (2.5)-(2.9) in the discretised wavenumber space k x ∈ [0, 2], k z ∈ [0, 30], and look for all the possible linear growing (Im(ω) > 0) modes. The (k x , k z ) domain is suitably moved toward lower (higher) wave numbers when the F r number is significantly higher (lower) than 1. In figure 4 we report the stability diagram. Each point in the diagram corresponds to the most unstable mode, whose relative k x and k z generally vary. One observes that at Re = 1000 the unstable region is relatively constrained around F r = 1 (i.e. 0.5 F r 2) but already covers two decades in F r at Re = 10000. This indicates that instability first (i.e. at low Re number) appears where density stratification and horizontal shear are comparable, i.e. N ∼ σ, but is likely to be observed in a sensibly wider range of the ratio σ/N provided that the Re number is large enough.

The critical Reynolds number (Re c ∼ 700) appears quite moderate compared to other unstratified parallel flows like the plane Poiseuille flow (Re c = 5772 according to [START_REF] Orszag | Accurate solution of the OrrSommerfeld stability equation[END_REF]). The value we find is comparable with that found by Chen (2016) for a plane Poiseuille flow in the presence of vertical stratification, but still sensibly lower than that indicated by the same authors [START_REF] Chen | Instability of a boundary layer flow on a vertical wall in a stably stratified fluid[END_REF] for the boundary layer (vertically stratified) profile Re c ∼ 1995. The growth rate is moderate even at high Re number, indicating that the observed instability is not only constrained in the (k x , k z ) but also relatively slow to establish.

Finally we want to discuss how the most unstable mode changes as a function of Re and F r numbers separately. In figure 5 we analyse how the growth rate Im(ω) of the most unstable mode (left) changes with the Re number at fixed F r = 0.4. One sees that Im(ω) rapidly saturates to a constant value. This result was confirmed by solving the eigenvalues problem (right) at very high Re number (up to 10 8 ) with (k x , k z ) fixed. In figure 6 (left) we report the value of k x and k z for the most unstable mode as a function of Re at fixed F r = 0.4. One sees that both k x and k z tend to a constant value. Let us recall that our approach demands to discretize the space (k x , k z ) which explains why the rate of this convergence may appear disturbingly abrupt. Thus we conclude that the observed instability must rely on an inviscid mechanism and that the inviscid approximation is sufficient to capture the spatial (k x , k z ) and temporal ω feature of the most unstable mode.

In figure 6 (right) we report the value of k x and k x k z for the most unstable mode as a function of the F r number at Re = 10000. The first panel shows that k x is always slightly lower than 1/F r (dashed line) which means that, for the most unstable mode, baroclinic critical layers (i.e. y = ±1/kF r) fall close to the boundaries but slightly outside the domain boundaries y = ±1, and are likely not involved in the instability mechanism. In the second panel we see that all solutions seem to collapse on the curve A/F r where A = k x k z | F r=1 which provides a rule for the spatial pattern of the most unstable mode and an interesting limit for further analysis of the pressure equation (2.10). Finally one should remark that, according to this relationship, in exploring the stability diagram (Re, F r), the discretization of the wave number (i.e. the step size of the grid kx, kz) becomes critical at low F r, while the size of the domain (k x , k z ) becomes critical at high F r.

The instability mechanism

So far we have only focused on the features of the most unstable mode for a given combination of the dimensionless numbers Re, F r and a typical domain in the wave number space (k x ,k z ). This characterizes the instability from an operational point of view but does not say anything about the underlying mechanism. To this end we now analyse the shape of unstable modes. We have seen that the asymptotic behaviour of the instability at large Re number indicates that it relies on an inviscid mechanism. In the inviscid limit the pseudo-spectral approach is far less intelligible because the solution of the eigenvalue problem contains a large number of spurious modes with Im(ω) > 0, which makes the detection of genuine unstable modes extremely difficult. The idea is then to consider a finite Re number to keep the eigenvalue problem manageable but also large enough to capture all the possible features of the instability diagram. It turns out that the choice Re = 10000 fairly responds to these criteria, thus we focus on the case F r = 1 and Re = 10000 as a reference one. In figure 7 we report the eigenfunctions of the most unstable mode at F r = 1 and Re = 10000, which corresponds to the wave numbers k x = 0.767 and k z = 4.937. One observes that the perturbations of the vertical velocity w and buoyancy b are more important close to the boundaries y = ±1 while at the center of the domain y = 0, the velocity perturbation is mainly horizontal. We consider now a sample mode for each different unstable branch, for example corresponding to the red spots we labelled with a,b,c,d, e and f in figure 3. In figure 8 we compare the pressure eigenmode for all different branches.

One observes that the shape of the eigenmodes is significantly different in each panel. Not surprisingly modes from the two stationary branches (a) and (d) are symmetric in the cross-stream direction y. Conversely, travelling modes (b,c,e,f ) are asymmetric but always appear in pairs, at ω ± = ±Re(ω) + iIm(ω), each mode in a pair being the y-mirrored of the other one with respect to y = 0. Also in panel a we superpose the pressure eigenfunction of the most unstable mode at k x = 0.767 and k z = 4.937, i.e. the same as figure 7 (last panel). One remarks that two pressure eigenmodes belonging to the same branch (i.e. the branch (a)) have basically the same shape. The scenario we described above is strikingly similar to that presented by [START_REF] Satomura | An investigation of shear instability in a shallow water[END_REF] (see e.g. his figure 6) who analysed the stability of a non-stratified PC flow in the shallow water approximation. In this case the pressure p is replaced by the elevation of the free surface h in the analogous of equation (2.10). The author suggested that the instability is produced by the resonance of two Doppler-shifted shallow water waves. In this picture the wave (stream-wise) phase speed C = ω/k x of a shallow water wave which travels close to one boundary can be approximated to that of a shallow water wave in a fluid at rest plus a Doppler shift, say U d , which has the sign of the velocity of the considered boundary. Two distinct counter propagating waves situated at opposite boundaries can then have the same phase speed and become resonant. Moreover the resonant wave numbers constitute a discrete spectrum because rigid walls make the dispersion relation of (non sheared) waves discrete. More recently, the same mechanism was also detailed to be responsible for linear instabilities in stratified, rotating plane Couette [START_REF] Vanneste | Unbalanced instabilities of rapidly rotating stratified shear flows[END_REF]) and stratified Taylor-Couette (Park & Billant 2013) flows. We suggest and show below that this interpretation remains valid in our case if we replace shallow water gravity waves with internal gravity waves. The dispersion relation of the latter is also discrete and one has:

C (n) ± = ± 1 k x F r 1 - k 2 z k 2 x + k 2 z + n 2 π 2 , (3.1)
where we use the same notation as [START_REF] Satomura | An investigation of shear instability in a shallow water[END_REF], i.e. C = ω/k x . Subscript + (-) refers to waves propagating in the positive (negative) direction of the x axis, while superscript n labels different channel modes. Note that here the velocity C does not correspond to the phase velocity of the wave nor to its horizontal component, but it is still the relevant quantity to describe the resonance mechanism.

In the first panel of figure 9 we report the value of C

(n)

± as a function of k x at k z = 10.77 and F r = 1. In the second panel we show both Doppler shifted velocities

C d = C (n) ± ∓ U d
where we consider prograde and retrograde waves moving upstream close to opposite boundaries and transported by the local mean flow. At this stage the value of U d is an adjustable parameter -1 U d 1 and was fixed to U d = 0.6. One remarks that in the fixed frame (i.e U d = 0), prograde waves (solid lines) are well separated from retrograde waves (dashed lines). On the contrary for Doppler shifted waves, there exists a discrete set of resonant k x where two curves of different type cross each other. At F r = 1 the first crossing (resonance) happens at C d = 0 and close to k x = 0.4, which is consistent with the appearance of the first stationary mode (a) in the stability map of figure 3. The next two resonances happen at a non zero value of C d which coherently recovers the appearance of the first two oscillating modes (b) and (c) at larger k x in figure 3. The following crossing happens again at C d = 0, which confirms the appearance of a fourth (stationary) unstable branch in figure 3 when moving along k x and at constant k z = 10.77. Looking back at figure 8 the mode (f) appears as a half of the mode (a) thus we speculate that the corresponding resonance originates from the crossing of a Doppler shifted wave and a non transported wave (i.e. one for which U d = 0) situated at the center of the domain (i.e. the dashed dotted lines in figure 9). Modes (e) do not originate from a resonance, consequently they are not indicated in figure 9. A closer inspection of the velocity field suggests that in this case the baroclinic critical layers are excited, and the instability relies on a different mechanism. This hypothesis is consistent with the fact that the mode (e) (see figure 3) belongs to a region which mainly extends at k x > 1/F r where critical layers can fit within the domain. Now that we have possibly explained the origin of all the distinct modes as resulting from a degeneracy of the Doppler shifted frequency, we want to show that this picture allows to fully capture the shape of the unstable branches in figure 3. First one should recall that for a given channel mode (i.e. n = const), the dispersion relation of internal gravity waves (3.1) is a function of two variables k x and k z , hence a surface. It follows that degeneracy occurs indeed on the intersection of two surfaces (i.e. not two curves) which is a curve (i.e. not a single point). The latter explains why the shape of the unstable branches in figure 3 appears elongated in one direction and constrained in the orthogonal one. In the particular case of a stationary mode, one can easily deduce the equation of such a curve from the dispersion relation (eq. 3.1) modified by the Doppler shift U d . We find:

k z = F(k x , n, F r, U d ) = n 2 π 2 U 2 d F r 2 k 2 x -k 2 x -n 2 π 2 + 1 F r 2 U 2 d .
(3.2)

In figure 3 we have superimposed the value of F to the map of the growth rate Im(ω) at different F r numbers and for n = 1 and 2. One observes the agreement is not only qualitative, for example F reproduces the trend k x k z ≈ const observed before, but also quantitative, because fixing a unique value of U d = 0.6, we are able to predict the position of almost all the unstable stationary branches. Finally we show that the mechanism we describe above allows to predict the boundaries of the unstable region. If we look back at figure 9 one observes that instability appears at a finite value of k x , say k inf , where

C (1) + (k x = k inf x , n = 1, k z , F r) = U d (3.3)
and must disappear when the envelopes of prograde and retrograde modes (red lines) cross each other, at k sup x = 1/F rU d . Note that the latter upper boundary is independent of k z . Conversely the lower boundary can be arbitrarily reduced, for example k inf x → ∞ provided that k z → ∞. Nonetheless any finite Re number will likely inhibit an instability happening at large wave number k z . We conclude that according to the proposed resonance mechanism, the instability is triggered by perturbations which are not streamwise invariant (i.e. k x = 0), and at stream-wise wave number k x < 1/F rU d . Looking at the growth rate diagrams of figure 3 one actually sees that k sup x tends to overestimate the upper bound of the unstable region. A closer estimation of the latter may then be given by the crossing of the prograde (retrograde) n = 1 Doppler shifted mode with the envelope of the retrograde (prograde) waves, which happens at k x < k sup x , in the region we highlighted with a red circle in figure 9. The idea is that multiple degeneracy may inhibit the resonance mechanism. Lower panels in figure 9 illustrate the same resonance mechanism for F r = 0.2 and F r = 5. The bottom left panel (F r = 0.2) confirms that the instability range is extended and pushed at larger k x for small F r number (i.e high stratification). Conversely the bottom right panel (F r = 5), shows that the region where resonances take place, both shrinks and is constrained to smaller k x . Ultimately k inf x and k sup x collide in the limit F r → ∞ and the instability likely disappears or at least reduces to an infinitely narrow range in k x . Note that the results above suggest that the upper boundary of the unstable region in figure 4 is intrinsic to the instability mechanism, while lower boundary is controlled by the Re number: at small F r instability appears at larger k x , thus larger k z and is then more sensitive to viscous dissipation. To conclude this section we recall that if the growth rate varies with the Re number and different branches appear at different Re numbers, the value of C = ω/k x on a same branch is approximately constant and almost does not vary with the Re number. This supports the hypothesis of a resonance and confirms that the appearance of the most unstable stationary and oscillating modes relies on an inviscid mechanism.

Effect of the Schmidt number

All the results we presented above correspond to solutions of the eigenvalue problem where mass diffusivity was completely neglected, that is Sc = ∞. We have modified the eigenvalue problem and tested the relevance of a finite Sc number for the reference case Re = 966 and F r = 0.82 which will serve as a comparison between linear analysis, experiments and direct numerical simulations. All the simulations are performed at the wave numbers k x = 0.96, k z = 5.16, where the most unstable mode appears in the Sc = ∞ case. The results are reported in figure 10. First we report that at Sc = ∞ (circles) and Sc = 700 (crosses) the eigenvalues are well superimposed. This suggests that our non diffusive approximation is qualitatively and quantitatively adequate to compare linear theory with experiments performed with salty water, for which Sc = 700. Second we remark that at Sc = 7 (squares) there is still an unstable mode and close to the origin the distribution of eigenvalues has the same form. For example looking at the close up on the right of figure 10, one sees that all the eigenvalues at Sc = 7 (in red) are located close to a non diffusive eigenvalue. This result makes possible the comparison between linear analysis, experiments, and direct numerical simulations which will be performed at Sc = 7. Finally we observe that at Sc = 1 (diamonds and stars) the eigenvalues are distributed on three distinct Y-shaped branches which is consistent with the previous study of [START_REF] Bakas | Gravity waves in a horizontal shear flow. part ii: Interaction between gravity waves and potential vorticity perturbations[END_REF] and Chen (2016) who found analogous branches at Sc = 1 in the case of the PC flow and the plane Poiseuille flow, respectively. We also remark that at Sc = 1, Re = 966 (diamonds) there is no unstable mode, nonetheless instability is promptly recovered at Re = 2000 (stars). We conclude that increasing mass diffusion, the threshold of the instability is not severely affected as long as Sc 7 while it may change when Sc is of the order of unity.

Direct Numerical Simulations

In addition to the linear stability analysis, we have performed Direct Numerical Simulations (DNS) of the full set of equations (2.1)-(2.3). The aim of a complementary DNS approach is to validate the linear theory and characterize the flow when retaining all the non-linearities. Equations are solved in a rectangular box of dimensions (L x , L y , L z ). The boundary conditions are periodic in both the stream-wise and vertical directions and rigid no-slip insulating boundaries in the cross-stream direction, i.e. u = 0 and db/dy = 0 at y = ±1. In order to keep the computational time reasonable, while still focusing on the high Sc number regime of the experiment described in section 5, we fix Sc = 7. We have seen in section 3.3 that this particular choice does not affect qualitatively the results, and in any case, ad-hoc solutions of the linear problem at Sc = 7 can be considered for a quantitative comparison. In order to ensure that the linear instability is well captured by the numerical simulation, we choose a box of size (

L x = 2π/k x , L y = 2, L z = 2π/k z ),
where k x and k z are the most unstable wave numbers as predicted by the linear stability analysis presented above. We performed DNS using the spectral element solver Nek5000 [START_REF] Fischer | An overlapping Schwarz method for spectral element solution of the incompressible Navier-Stokes equations[END_REF][START_REF] Fischer | Simulation of high-Reynolds number vascular flows[END_REF][START_REF] Paul | nek5000 Web page[END_REF]. The use of spectral elements instead of more classical pseudo-spectral methods will be justified later (see section 5) where we add the effect of the stream-wise confinement to mimic the experimental setup. The global geometry is partitioned into hexahedral elements, with refinement close to the moving boundaries. Velocity, buoyancy and pressure variables are represented as tensor product Lagrange polynomials of order N and N -2 based on Gauss or Gauss-Lobatto quadrature points. The total number of grid points is given by EN 3 where E is the number of elements. For all the results discussed in this paper, the number of elements is E = 6336 and we use a polynomial order from N = 7 up to N = 11 for the highest Reynolds number case. Time integration is performed with a third-order explicit scheme for the advection and buoyancy terms while viscous and dissipative terms are integrated using an implicit third-order scheme. The simulations are initialized with a small amplitude buoyancy perturbation and with an established linear PC flow.

In order to validate the eigenvalue problem we choose the reference case Re = 966, F r = 0.82 which will serve later as a comparison to experiments. In figure 11 (top left) we report the time evolution of the vertical kinetic energy density (thick line) which is defined as:

(w 2 ) 1/2 = 1 V V w 2 dV 1/2 (4.1)
where V refers to the volume of the simulation box. The quantity (w 2 ) 1/2 is appropriate since (w 2 )

1/2 t=0 = 0 for the base flow. One observes that (w 2 ) 1/2 increases exponentially and superposing the exponential growth predicted by the linear analysis one obtains an excellent agreement, with a relative discrepancy on the growth rate σ c of less than 1%. In the same figure (bottom left) we also report the spatio-temporal diagram of the One observes that a stationary pattern appears close to t = 500. Right: instantaneous 3D map of the buoyancy perturbation once the flow has become unstable (DNS). As predicted by the linear analysis the selected mode is mainly modulated in the vertical direction but still not stream-wise invariant (i.e. kx = 0). Note also that perturbations are concentrated near the boundaries y = ±1.

horizontal perturbation u at x = y = 0. One observes that a stationary pattern has established around t = 500 which has a well defined vertical wavelength. In figure 11 (right) we report a visualisation of the buoyancy perturbation b once the instability has saturated. One observes an weakly inclined layering of the density field which is a common feature in stratified turbulent shear flows (see Thorpe 2016, for a review). Again we have a very good agreement with the linear theory: a distinct spatial pattern appears and both vertical and horizontal wavelengths correspond to the predicted values. One can notice that the spatial pattern perfectly fits in the simulation domain. This condition is indeed necessary to observe the instability, and no relevant growth of the vertical kinetic energy is observed when none of the unstable wave numbers fits inside the simulation domain.

Experiments

Experimental apparatus

Now, we want to study whether or not this linear instability of the stratified plane Couette flow does appear in a "real" configuration, and to do so, we look for some of its signatures in an experimental set-up, intrinsically limited in size. The flow is produced with a shearing device which is placed inside a transparent tank (50 cm x 50 cm x 70 cm) made of acrylic. The tank is filled with salty water linearly stratified in density. The shearing device is sketched in figure 12 (left). The device consists of a PVC transparent belt (0.8 mm thick) which is closed on a loop around two vertical entraining cylinders made of dense sponge (we use standard spares entraining cylinders for commercial swimming-pool robots). Two additional pairs of cylinders (inox, 2 cm diameter) constrain the two sides of the loop to be parallel and at a controlled distance d. All cylinders are mounted on a system of acrylic plates which allows to vary the distance between the entraining cylinders (i.e. to tighten the belt) through two pairs of coupled screws (i.e. one pair for the bottom and one for the top). The top acrylic plates also prevent the existence of a free surface which would affect any imaging from the top. The motion of the belt is provided by a motor which is mounted on the top of the device and joined to the axis of one of the entraining cylinders. Finally two PVC rigid plates are mounted vertically in front of the two entraining cylinders in order to reduce at most any perturbations coming from the entrainment system. The distance between the edges of the plates and the belt is a few mm. Thus we look at the flow in the area shaded in light grey (figure 12, left). In the present work we consider two values of the gap width d = 5.8 and 9.8 cm, while the distance between the PVC plates D was respectively 34 cm and 24 cm, leading to a value of the aspect ratio D/d of 5.7 and 2.4 respectively.

The tank is filled with salty water of variable density. As a general rule a water column of height H = 10 to 20 cm linearly stratified in density always occupies the volume delimited by the belt and the confining barriers, while above and below the density stratification was generally weaker or negligible. The density profile is obtained by the double-bucket method [START_REF] Oster | Density gradients[END_REF]. To measure the density profile we collect small samples of fluid (∼ 10 ml) at different heights and analyze them with a density-meter Anton Paar DMA 35. The Brunt-Väisälä frequency N is constant for each experiment with a value between 0.5 rad/s and 3.0 rad/s. We measure the stratification before and after each experiment. The shearing motion clearly affects the stratification especially through the small scale features of the rotating part of the device which necessarily produces some mixing. Also, in our highest Re experiments we observe optical distortion which may indicate the presence of high density gradient zones and thus density layering, for example similarly to that observed in turbulent stratified experiments performed in Taylor-Couette devices [START_REF] Oglethorpe | Spontaneous layering in stratified turbulent TaylorCouette flow[END_REF]). Nevertheless we observe that the density profile at the end of an experiment is weakly perturbed and the relative discrepancy in the area of interest is around 5%. Finally we assume the viscosity to be ν = 10 -6 m 2 /s, and neglect any change associated to variable salt concentration.

The fluid is seeded with (10 µm -diameter) hollow glass spheres and two laser sheets illuminate the particles in the vertical plane y = 0 and the horizontal plane z = 0 as shown in figure 12. The flow is then recorded from the side by a 4 Mpx camera at a frame rate of 8 fps and from the top by a 2 Mpx camera at a frame rate of 30 fps. The velocity field is obtained with a Particle Image Velocimetry (hereafter PIV) crosscorrelation algorithm [START_REF] Meunier | Analysis and treatment of errors due to high velocity gradients in particle image velocimetry[END_REF]. Note that the mid vertical plane y = 0 is the appropriate place to detect the possible onset of an instability because in the ideal PC regime, the velocity should be zero there. The current setup permits only one by one enlighting-recording of the flow, thus movies from the top and from the side are always taken at different times.

Base flow

First we report that a PC flow can be observed in the region confined between the belt and the PVC barriers. On the top left of figure 13 we superpose 40 images of the z = 0 plane exactly as captured by the camera. Only the contrast was altered to exalt streamlines. Both the intersections of the belt with the laser sheet and the left barrier edge can be easily recognized as brighter lines. One also sees that streamlines close up near the PVC barriers and recirculations are present. This is confirmed by the velocity field given by the PIV algorithm and shown just below. The velocity plot is obtained by averaging over 40 PIV fields (∼ 1.3 s), also we plot only one arrow over four in the horizontal direction, to make the diagram readable. One remarks that up to 10 cm far from the center the flow is nicely parallel and the velocity gradient is linear. Both streamlines and PIV fields refer to an experiment where the base flow was already stationary. Now, any experiment necessarily implies a transient phase where the flow evolves from a first stationary phase, e.g. the whole fluid is at rest, to a second stationary phase which is the forced parallel flow. We expect the base flow to establish via viscous entrainment starting from the fluid layers which are close to the walls, thus the viscous time T ν = d 2 /ν seems to be an appropriate time scale for the transient. In order to verify this we need some more quantitative prediction and consider the transient flow generated by two infinite walls treated by [START_REF] Acheson | Elementary Fluid Dynamics[END_REF]. As a first step, recirculations are neglected. If the flow is initiated at t = 0 the horizontal velocity has the form U (y, t) = U 0 (y) -U T (y, t), where U 0 (y) is the asymptotic base flow and the transient part U T (y, t) reads:

U T (y, t) = (U 0 -U i ) ∞ j=1 2 π
(-1) j j e -π 2 j 2 t/Re sin jπy (5.1)

where U i is the velocity of the belt at t = 0, for example U i = 0 if the experiment is started with the fluid at rest. In figure 13 we compare the value of U (y, t) as expected from equation (5.1) with the average value of the horizontal velocity as observed in a typical experiment. The value of U is plotted as a function of y at four different times. First, the velocity profile collapses on the expected PC flow (dashed line) around t = T ν , which confirms that the base flow establishes via viscous entrainment. Also at t = T ν /3 (circles), the value of the average horizontal velocity is already very close to the PC flow. Secondly one remarks an excellent agreement of the experimental observations with the infinite walls approximation, which suggests that the recirculation does not affect significantly the shape of the transient flow. With regard to this, one should notice that knowing the time the base flow needs to establish becomes crucial when determining the growth rate of the instability, which will be discussed later.

Once the PC profile is established we want to detect possible deviations from the base flow. For this aim we mainly focus on the mid plane y = 0 where no motion is expected for the base flow. As a standard protocol we initiate the flow at low shear rate σ and then increase σ by a small fraction (typically 15%). Top views of the plane z = 0 are also taken to verify the shape of the parallel base flow. In each experiment the flow was observed for at least one viscous time T ν = d 2 /ν which may be taken as an upper-bound for establishing the base flow. First we report that starting from very moderate Reynolds number Re 300 the observed fluid oscillates coherently at a well defined frequency. In figure 14 we report the spatio-temporal evolution of the horizontal velocity perturbation u along the vertical line y = 0, x = 0 and the two horizontal lines y = 0, z = 0 and z = 0, x = -d. As visible in each diagram, parallel periodic structures appear which are fairly homogeneous over the spatial domain. This suggests that the entire fluid bulk oscillates in a coherent way. In addition all the diagrams show a well defined temporal frequency. We then perform the temporal fourier transform of the vertical average of u(x = y = 0, z, t) and denotes with f box the peak in the frequency spectrum. In the last panel of figure 14 we report the value of f box as a function of the imposed shear σ for a collection of three experiments at different F r numbers, where only the imposed shear σ is changed. One sees that the observed global frequency f box increases linearly with σ, and rescaling the observed period with the belt revolution time T rev , we find (see inset) that the observed oscillation period T box is very close to T rev /4. A possible explanation for this unexpected observation is that the two pairs of confining cylinders (see figure 12) divide the path of the belt in four (almost) equivalent sections. We stress that these cylinders constitute one of the biggest source of noise because they are rigid and each of them tends to perturb the belt once per revolution time, when they touch the roughness of the belt junction. Correspondingly we expect that a T rev /4 resonant periodic forcing may establish, and give the velocity field the observed temporal pattern. As a summary we report that deviation from a zero velocity field are observed in the horizontal velocity perturbation u, from a very moderate value of the Re number. Nonetheless the observed motion shows a trivial spatial pattern (i.e. k x = k y = k z = 0), looks like a bulk oscillation, and seems to be connected to the shearing device. Thus, in the presence of only spurious bulk oscillations (even if robust), we assess that the corresponding Re, F r pair is a stable point in the stability diagram. In the following we discuss how perturbations become more finely structured at higher Re number, and we give a criterion to distinguish these initial deviations from a truly unstable pattern.

Instability

When the Re number is sufficiently high Re 1000 and for Froude number F r ∼ 1, an exponentially growing motion is observed to form in the mid plane y = 0. The horizontal velocity perturbation u shows a well defined spatial pattern where horizontal and vertical wavelengths λ x , λ z can be fairly detected, with λ x /λ z ∼ 8. Results for this reference case are summarised in figure 15.

On the top of figure 15 we plot the horizontal velocity perturbation u at x = y = 0 as a function of the time t and vertical direction z for a reference (unstable) experiment. At t ∼ 60 s the imposed shear has changed from a lower value of σ = 0.34 s -1 to 1.15 s -1 and at t ∼ 600 s one observes the appearance of a vertical wavelength. One also observes that the t -periodic and z -invariant bulk motion described in the previous section is present since the very beginning and is still visible at large time, superposed to the instability pattern. In figure 15 (bottom left), we consider the time evolution of the order parameter (u 2 ) 1/2 for the same unstable case as above (black line) and for another case where (red line) the imposed shear, σ is smaller by a fraction 1/8. The average square of the horizontal perturbation u 2 is computed at the center vertical line x = 0, y = 0 in three steps starting from the spatio-temporal diagram of u. At each time we take the average of u over a short interval ∼ 4 σ -1 , subtract the linear regression in z, compute the square and finally average over the vertical direction. We stress that subtracting the linear regression allows us to get rid of the bulk oscillations and of any possible top-bottom anisotropy due to a non perfect verticality of the laser sheet. First of all one remarks that the unstable case (black line) shows a clear growth event of the order parameter which does not happen for the stable case (red line), thus indicating the appearance of an instability. Focusing on the unstable case one clearly sees that a first increase of (u 2 ) 1/2 occurs during the interval of ∼ 0.3T ν after the change in the imposed shear, where the . Bottom left: Evolution of the mean horizontal perturbation (u 2 ) 1/2 as a function of time for the same experiment (black line) and for a stable experiment where both F r and Re are diminished by a fraction 1/8 (red line). The inset shows log (u 2 ) 1/2 for the unstable case. At each time u 2 is obtained averaging over a short interval ∼ 4 , taking the square and finally averaging over the vertical direction. Bottom right: Snapshot of the horizontal velocity perturbation u in the plane y = 0, here t ∼ 900 .

value of the viscous time is

T ν = d 2 /ν ∼ 1000 σ -1 .
At larger time, u 2 increases again, now in an exponential way (see the inset in semi-log scale), and finally saturates to a constant value. The exponential growth rate is approximately ω ∼ 0.06 σ -1 (although the noise makes difficult a precise measurement of the growth rate). We claim that the first growing phase coincides with the progressive onset of the base PC flow at the imposed shear while the second growing phase corresponds to the onset of a linear instability. Note that for the stable case (red line) the first growing phase is less visible, because the shear σ is imposed starting from a slightly lower value. Finally on the bottom right of figure 15 we present a snapshot of the u field in the plane y = 0. One observes a regular periodic pattern characterized by a vertical wavelength λ z ∼ 0.7 L 0 and an horizontal wavelength λ z ∼ 5.5 L 0 , where we recall that L 0 is half the width of the channel. Below we consider the stability of our experimental flow in the (Re, F r) space which is the same for the linear stability analysis performed in section 3.1. To this aim we need to define a common protocol to assess the presence or not of the instability. One criterion may be the appearance of a vertical length. Unfortunately the latter is a smooth process, for example a vertical wavelength was often visible at a shear σ lower than what we assess to be the unstable case. Nonetheless the associated signal was generally weak and no growth process was observed. The existence of the latter seems to be the most reliable criterion, but demands longer experiments and generally imposes to start from very small σ, which imposes to slow down the flow after each experiment. As a general rule we rather look at the saturated amplitude of the order parameter (u 2 ) 1/2 as a function of σ and detect if an abrupt change occurs, as it is clearly visible in figure 15 at large time. In particular (u 2 ) 1/2 is computed once the instability has saturated or alternatively after a time of the order of the viscous time T ν after the actual value of the shear is imposed, for example to get rid of the base flow transient. We choose the control parameter σ (i.e. the imposed shear) as the most suitable one, because it can be varied continuously, simply controlling the speed of the entraining motor. As a drawback both Re and F r are linear in σ, thus the stability diagram must be explored moving on tilted straight lines for which the ratio F r/Re = ν/N d 2 is constant. Any change in the vertical stratification N and gap width d is considerably more laborious, which constrains the exploration of the (Re, F r) space to a few different F r/Re = const lines. In figure 16 we report the value of (u 2 ) 1/2 as a function of the control parameter σ for a collection of 6 series of experiments at different vertical pattern establishes but appears less structured and less robust than in the other cases. Also in the case N 3 the transition stable-unstable is anticipated by one position because a growth event was clearly visible at the corresponding value of σ.

Discussion

In figure 17 we superimpose the experimental results of figure 16 on the stability diagram obtained with the linear stability analysis (figure 4). Experimental points which belong to the same line correspond to a series of observation at fixed F r/Re ratio.

Closed symbols and open symbols correspond to stable and unstable cases respectively, as in figure 16. One observes that the transition of the quantity (u 2 ) 1/2 (i.e. from close to open symbols) happens close to the marginal contour where linear growing modes appear according to linear stability analysis. This strongly supports the claim that we experimentally observe the signature of the instability predicted by the linear analysis.

Below we compare the temporal behaviour of the observed instability with the linear analysis. Besides the growth rate, that precisely characterizes the instability onset, we want to discuss first what happens during the transient phase that necessarily comes with each experiment. This constitutes a difference with the linear analysis where the base flow is always constant, and may affect the estimation of the observed growth rate. In other words, one may wonder at which time since the beginning of an experiment the instability is expected to grow. The question becomes particularly relevant when considering that the expected growth rate is comparable with and even smaller than the viscous time. Also we want to rule out the possibility that the appearance of the unstable pattern is due rather to the transient profile of our flow, for example a non constant shear profile like that considered by [START_REF] Chen | Instability of a boundary layer flow on a vertical wall in a stably stratified fluid[END_REF]. The temporal diagram of figure 15 shows that the exponential growth starts at t ∼ T ν /3, which seems to be consistent with the Here we choose the same Re, F r, and Ui (i.e. the belt velocity at t = 0) as the reference experiment. The eigenvalues problem is solved for the combination (kx, kz) which is the most unstable according to linear analysis and at seven different transient times. Different symbols correspond to different times. The dashed line marks the limit for instability.

description of the base flow given in section 5.2. In order to give a more quantitative explanation we solved a modified eigenvalue problem where the base flow is now given by the expression (5.1) for a collection of different times and with the same set of parameters as the reference experiment presented in figure 15. The results are presented in figure 18, where we report the eigenvalues of the most unstable mode focusing close to the transition region Im(ω) = 0. One remarks that no unstable eigenvalue is present for t < 0.2T ν while one unstable mode appears for t 0.3T ν , thus confirming that the base flow must be sufficiently established for the instability to develop. This result was confirmed by specific DNS where the initial condition is not the PC flow, but the flow is at rest and the shear profile is progressively established through the no slip boundaries. Also in this case the growth of perturbations is delayed to the moment when the shear profile has become almost constant. We then conclude that what we observe is associated with a constant shear plane Couette profile.

Besides instability threshold in the Re, F r space we also want to compare the shape of the mode selected in our experimental device with the unstable mode predicted by linear theory. If we focus again on the reference case described in figure 15, we remark that both k x and k z are larger than what is predicted by the linear theory for the most unstable mode, and the perturbation is oscillating in time while linear theory predicts that instability appears as a stationary mode. One remarks that oscillations are quite regular and relatively slow, with a typical period T = 43 ± 3 σ -1 , where the uncertainty is taken as the width at mid height of the peak in the average temporal spectrum. We recall that according to linear analysis oscillatory branches also exist (see figure 3) which appear at higher Re number, typically Re 2000. Interestingly the period associated to the first oscillating branch is always long, for example with the F r of the reference experiment N 1 and Re = 2500 one has T ∼ 42 σ -1 . Moreover this branch happens at larger wave numbers more compatible with the observed ones. In this scenario what we look at may be either a single propagative mode, for example of the type (b), as it is visible in the lower part of the spatio-temporal diagram, or a standing wave generated by two counterpropagating modes of type (b) as it is visible in the upper part of the same diagram. These elements suggest that experimentally, the instability is activated close to the absolute threshold (i.e. that where first stationary modes appear) predicted by linear analysis, but a different non stationary mode is selected in the end.

In any case one should recall that our linear and experimental problems are different, thus unstable modes are not expected to share the same features. The major difference between the theoretical system and our experimental setup is the finite size of the domain. In principle, to mimic periodic boundary conditions imposed in linear calculations, one wants to take the horizontal and vertical aspect ratio D/d

1 and H/d 1 while our best realization (i.e. narrow gap) of this hypothesis was D/d = 5.7 and H/d = 2.4. We observe that the impact of physical confinement is twofold. First from the point of view of modal analysis only the wavelengths which fit in the domain may have a chance to develop. This was confirmed by periodic DNS that do not show any instability whenever the box size does not fit the spatial shape of unstable modes. We notice that the eigenvalue problem is solved assigning an arbitrary value of k x and k z , thus when comparing to DNS and experiments one should retain that the (k x , k z ) grid of figure 2 is coarsed-grain, especially at low k x and k z . Thus the ideal constraint for the aspect ratio are D/d λ x and H/d λ z , where λ x and λ z are the non dimensional wavelengths of the unstable mode we want to observe. A second problem appears in the stream-wise direction because the streamlines must turn and close up when getting close to the walls that close the domain in the stream-wise direction, as it is clearly visible in the snapshot reported in figure 13. We mentioned before that this feature does not modify significantly the shape of the base flow in the bulk, but it may locally destabilize the flow (i.e. close to the corners) and successively affect the stability of the whole domain.

Simulations in a finite domain

In order to closer investigate finite size effects, we performed new DNS where the computational domain is now closed in the stream-wise direction by two solid walls with no-slip insulating boundary conditions as sketched in figure 12 (right), while boundary conditions remain periodic in the vertical direction. Note that compared with the periodic case discussed in section 4 the mesh is further refined close to the two additional streamwise boundaries, in order to properly solve for the boundary layers. In addition, the corners of the domain are now singular due to the incompatibility between the velocity imposed at the side boundaries and the fixed stream-wise walls. This is naturally smoothed by viscosity but is nevertheless an inevitable source of vorticity. In figure 19 (top) we report the results of a DNS confined simulation which reproduces both the control parameters (F r,Re) and the aspect ratio of the reference experiment illustrated in figure 15. For a direct comparison we report again the results of the reference experiment (bottom) already shown in figure 15 with the only difference that bulk oscillations are now filtered from the spatio temporal diagram of u at x = y = 0 and the origin of time axis is shifted forward to t = T ν /3 which is when we estimate that the PC flow is well established. One observes a striking good agreement between our DNS results and experimental result on both the spatial and temporal shapes of the selected mode. Computing the temporal spectrum we find that the temporal frequency as predicted by DNS is T = 50 ± 5 σ -1 which is compatible with the experimental one, while the consistency of spatial wavelengths is evident because in both cases an integer number of velocity maxima fit in the vertical mid plane y = 0. Finally we observe that as a whole the transition from the initial noise to the final non linear pattern takes almost the same time in DNS and the experiment. As a summary, figure 19 indicates that we correctly isolated the crucial factor which possibly alters the selection of the unstable mode, that is the stream-wise confinement. Incidentally we also report that additional DNS show that the form of the late non-linear stage is quite sensitive to initial and boundary conditions. For example slightly varying the box dimensions or the amount of initial noise the spatial shape of the selected mode is different and travelling waves or standing waves patterns can be alternatively present. To better investigate the role of stream-wise boundaries we perform additional DNS with a non stratified (F r = ∞) PC flow in exactly the same confined geometry and same Re number as the one we just described (figure 19). Results are reported in figure 20. We recall that u perturbations do not grow when considering periodic boundary conditions at F r = ∞, as expected by the fact that the unstratified PC flow is linearly stable. One observes that a strong vertical shear appears in a shorter time compared to both the stratified DNS and experiment in figure 19 while the vertical length scale is larger. A similar pattern is also observed in unstratified experiment as soon as Re 300. The vertical kinetic energy (black dashed line in figure 21) grows in an exponential way suggesting that a linear instability may act at a first stage. Conversely non linear mechanisms must act at larger time because the vertical kinetic energy rapidly saturates and the spatio-temporal diagram of figure 20 shows the presence of small scale features. A detailed investigation of the nature of this instability of a confined and unstratified PC flow is beyond the scope of the present work. Formally we recognise the scenario of a locally perturbed PC flow already described in literature (see Barkley & Tuckerman 1999, and citations inside), where some of the (stable) modes of the ideal PC flow become unstable. In particular the shape of the flow at the corners suggests that it may locally destabilise via centrifugal instability which will be studied in further investigations.

We conclude that the presence of boundaries may destabilize the PC flow with or without the presence of stratification. Now, the careful reader will agree that even if DNS fully justify differences between the linear analysis and the observed experimental pattern we are left with a cumbersome question regarding the origin of the perturbation pattern observed in the confined and stratified configuration: does this pattern coincide with a boundary-induced modification of the linear instability, or rather with the pure hydrodynamic boundary-induced instability modified by the stratification? For example one may claim that the patterns observed in figure 19 are a convoluted combination of the two showed in figure 11 and 20. The key to this answer resides in the same approach we followed with experiments, that is to detect if and when, our flow abruptly changes when varying the control parameter σ. We then consider further DNS which copy the parameters of another experiment of our reference series N 1 where F r and Re were 30% smaller than the unstable case described in figure 15. At the same time we repeat unstratified simulations at such a lower value of the Re number. We observe that the new unstratified case is almost unchanged while the subcritical stratified case shows a dramatic change. In this case perturbations are significant only close to the boundaries, and no instability develops in the bulk. In figure 21 we report the time evolution of the vertical kinetic energy for all the stream-wise confined DNS we have discussed above, together with the one performed in a periodic domain. One observes that in the unstratified case (dashed line), perturbations rapidly grow and saturate at the same value independently of the Re number. If we add stratification (solid lines), perturbation grows (black line) and saturates to the same value as periodic simulation if Re number is beyond the threshold predicted by the linear analysis. Conversely perturbations are damped (dashed line) when the Re number is below the threshold.

We have now enough elements to conclude that what we observe both in experiments and DNS (figure 19) is a signature of the linear instability of a PC flow, vertically stratified in density. Besides we observe that stream-wise boundaries are source of instability for a flow and likely affect the features, or possibly just the selection, of the unstable mode which shapes in the end the observed pattern. Incidentally we report that additional DNS were performed in a domain which is larger but still confined in the stream-wise direction, in order to explore when finite size effects become negligible and suitably design a larger experiment. Surprisingly we find that for a doubled size domain, the growth rate decreases and almost matches that of periodic simulations but the instability disappears (we observe a pattern similar to the subcritical stratified case) when further increasing the stream-wise domain (i.e. 4 times larger). With respect to this trend the long computation time demanded to consider even larger domains prevents us to be conclusive and further studies will be necessary.

The transition from the confined to the periodic case happens in a discontinuous way which needs to be further investigated.

At this stage we speculate that stream-wise boundaries may both introduce some forcing and inhibit the instability, perturbing the normal form and thus the resonance of the waves supported by the flow. Horizontal aspect ratio possibly controls the mutual importance of these two effects in a non trivial, non monotonous way, thus explaining the observed scenario.

Conclusions

We performed the linear stability analysis of the plane Couette Flow for a stably stratified fluid with a constant density gradient orthogonal to the shear. The domain has rigid closed boundaries in the direction of the shear, and open periodic boundaries in both vertical and stream-wise directions. Unstable stationary modes are found at strikingly moderate Reynolds number Re 700 and for a Froude number close to 1 for non vanishing horizontal and vertical wavenumbers with k x /k z ∼ 0.2. We then explore the stability of the flow in the (Re, F r) space. In the region we consider, the most unstable mode is always stationary and the growth rate remains relatively small, while the range of unstable F r numbers increases when increasing the Re number. Moreover the flow is unstable only to three dimensional perturbations, i.e. only for k x ,k z = 0. This result constitutes a fundamental difference with homogeneous shear flows for which the Squire theorem prescribes that the most unstable mode should be twodimensional. In the presence of stratification, hyperbolic tangent profile [START_REF] Deloncle | Three-dimensional stability of a horizontally sheared flow in a stably stratified fluid[END_REF] and Kolmogorov flow (Lucas et al. 2017) are also dominated by two-dimensional perturbations, while both boundary layer profile [START_REF] Chen | Instability of a boundary layer flow on a vertical wall in a stably stratified fluid[END_REF]) and strato-rotational instability [START_REF] Yavneh | Non-axisymmetric instability of centrifugally stable stratified Taylor-Couette flow[END_REF] similarly appear for three-dimensional perturbations. Curiously the study of [START_REF] Bakas | Gravity waves in a horizontal shear flow. part ii: Interaction between gravity waves and potential vorticity perturbations[END_REF] included the same modal approach but no unstable case was reported. Even more recently [START_REF] Chen | Instability of a boundary layer flow on a vertical wall in a stably stratified fluid[END_REF] showed that a linear instability does occur in parallel flows free from inflection points, when a vertical stratification is added, but only the plane Poiseuille flow and the viscous boundary layer were considered. Remarkably our instability also appears close to F r = 1, confirming the necessary coexistence of shear and stratification. The critical Reynolds number for the stratified plane Couette flow turns to be at least two times smaller than for the boundary layer and slightly larger than for the Poiseuille flow.

Looking at the most unstable mode of the linear problem, vertical velocity and density perturbations develop close to the boundaries which suggests that a crucial role may be played by lateral boundaries. Nonetheless, a comparable horizontal motion dominates in the mid vertical plane, and shows a vertically modulated pattern which is reminiscent of the deep equatorial currents, and staircase density layering in the Earth ocean [START_REF] Dunkerton | On the inertial stability of the equatorial middle atmosphere[END_REF][START_REF] Dengler | Equatorial deep jets and abyssal mixing in the indian ocean[END_REF]d'Orgeville et al. 2004).

A mechanism was proposed to explain the onset of the instability as the one suggested by [START_REF] Satomura | An investigation of shear instability in a shallow water[END_REF] adapted to the case of internal gravity waves, instead of shallow water waves. In this picture internal gravity waves are trapped close to the boundaries and Doppler shifted, thus allowing two counter propagating waves to become stationary and mutually resonant. The shape of the unstable region in the wave number space and the appearance of discrete additional resonances are also fully captured by the model, thus supporting its relevance. An analogous mechanism was also invoked at the origin of Strato-Rotational instability both in the plane Couette [START_REF] Kushner | Coupled Kelvin-wave and mirage-wave instabilities in semigeostrophic dynamics[END_REF][START_REF] Vanneste | Unbalanced instabilities of rapidly rotating stratified shear flows[END_REF]) and the Taylor-Couette (Yavneh et al. 2001;[START_REF] Park | The stably stratified TaylorCouette flow is always unstable except for solid-body rotation[END_REF] geometries.

The linear stability analysis was confirmed by DNS which fairly reproduce the spatial pattern and the growth rate. We report that no instability is observed when none of the unstable modes can properly fit in the domain. This confirms that the instability sharply selects the spatial pattern of the perturbation.

We analysed the experimental flow produced by a shearing device immersed in a tank filled with salty water linearly stratified in density. We report that when the F r number is close to 1 and for Re 1000 velocity perturbations are observed to grow in an exponential way. Remarkably we observe that perturbations start to grow only when the plane Couette profile is almost completely established. This was confirmed by ad-hoc versions of the linear problem for a collection of transient profiles and by DNS which mimic the transient flow of experiments. We conclude that the observed instability is crucially associated to the shape of the shear, namely the plane Couette profile.

Then we explored the stability of the flow in the (Re, F r) space varying the control parameter σ, that is equivalent to move along F r/Re = const lines, for a few different values of F r/Re. For each series of experiments we observe that an abrupt increase in the perturbation amplitude occurs, when σ is bigger than a threshold value σ c . When adding experimental data to the stability diagram predicted by linear theory we find that the threshold contour indicated by experiments qualitatively matches the margin of the linearly unstable region. Also, close to the threshold, the velocity perturbation shows a well organized pattern and is almost horizontal, which is in agreement with the solution of the linear problem. Nonetheless the unstable mode slowly oscillates in time and appears at higher wave numbers than the most unstable (stationary) mode indicated by the linear analysis. These two elements suggest that the mode selected in our experiment is not the most unstable of those predicted by linear analysis, or that these latter are possibly not the same when considering the finite size experimental apparatus.

In any case we claim that the origin of the discrepancy relies on the critically low value of the horizontal aspect ratio of our experimental domain, which is necessarily bounded in the stream-wise direction. The relevance of this hypothesis has been tested with complementary DNS where no slip rigid boundaries are now implemented also in the stream-wise direction. Remarkably, when copying the aspect ratio of our experiments we minutely reproduce the perturbation pattern observed in experiments. More generally DNS show that stream-wise confinement affects the stability of the flow irrespectively of the F r number (i.e. also without stratification), which questions the link between the instability observed in experiment and that predicted by linear analysis. We performed then DNS of a subcritical stratified experiment (i.e. Re and F r below the critical value) and show that the instability disappears. We then acknowledge the unstable pattern observed in both experiments and DNS as a true signature of the linear instability of a plane Couette flow vertically stratified in density.

Future studies are planned to closer investigate which is the critical aspect ratio to recover quantitatively the results of linear theory and periodic DNS. To this aim new DNS will be performed in a larger domain, which will possibly indicate how to correspondingly design a new set-up.

Quantitative measurements of the density field will be also performed in future experiments to quantify the density layering whose evidences were already available in our highest Reynolds experiments in the form of regularly spaced optical distortion. Such measurements will possibly add cues to the comprehension of the diapycnal mixing in the presence of horizontal layering as recently studied with experiments [START_REF] Woods | Non-invasive turbulent mixing across a density interface in a turbulent Taylor-Couette flow[END_REF][START_REF] Oglethorpe | Spontaneous layering in stratified turbulent TaylorCouette flow[END_REF]) and numerical simulations (Lucas & Caulfield 2017) in the case of the Taylor-Couette and Kolmogorov flows respectively.

Figure 1 .

 1 Figure1. Sketch of the analysed flow in a Cartesian reference x, ŷ, ẑ. The base flow is aligned with the stream-wise direction x, the constant shear is aligned with the cross-stream direction ŷ while density stratification and gravity are aligned with the vertical direction ẑ. We highlight in grey no slip lateral boundaries. Open periodic boundaries are not coloured.

Figure 2 .Figure 3 .

 23 Figure 2. Left: growth rate Im(ω) of the most unstable mode in the space (kx, kz) at F r = 1 and Re = 1000. Center: oscillation frequency Re(ω) of the most unstable mode in the space (kx, kz) at F r = 1 and Re = 1000. Notes that unstable modes (Im(ω) > 0) are stationary (Re(ω) = 0). Right: full spectrum at the most unstable mode kx = 0.815 and kz = 4.937, F r = 1, and Re = 1000. Crosses refer to Ny = 129 collocation points, circles to Ny = 257 and triangles to Ny = 513. The inset at the bottom right coincides with the area delimited by the red rectangle.

Figure 4 .

 4 Figure 4. Growth rate Im(ω) of the most unstable mode in the space (Re, F r). Each point is obtained taking the maximum value of Im(ω) over a collection of runs at fixed (Re, F r) and variable wave numbers (kx, kz). White dashed contours correspond to Im(ω) = 0.01,0.02,0.03,0.04,0.05,0.06. Black dashed lines correspond to Re = 10000 and F r = 0.4. White circles correspond to the points of the diagram analysed in figure 2 and 3.

Figure 5 .Figure 6 .

 56 Figure 5. Left: Im(ω) (growth rate) of the most unstable mode as a function of the Reynolds number at F r = 0.4. The dashed line correspond to the inviscid solution as obtained solving the eigenvalue problem. Right: solutions of the eigenvalue problem at F r = 0.4 with fixed kx = 1.29 and kz = 8.53. Different symbols correspond to different Re numbers. The inset corresponds to the thin rectangular region indicated by the red dashed line in the main graph. We highlight in red the two lowest Re numbers. One sees that starting from the third one (Re = 50000) the value of Im(ω) saturates to an asymptotic value.

Figure 7 .Figure 8 .

 78 Figure 7. Eigenfunctions of the most unstable mode at F r = 1, Re = 10000 and wave numbers kx = 0.767, kz = 4.937. Velocity fields and buoyancy are rescaled dividing by the maximum value of the perturbation u. Solid lines refer to the absolute value while dashed and dashed dotted lines refer to the real and imaginary part respectively.

F r = 1 FFigure 9 .

 19 Figure 9. Top: velocity C (n) ± as a function of kx as given by equation (3.1) at kz = 10.77 and F r = 1 and for the first n = 1 to 20 prograde (solid lines) and retrograde (dashed lines) confined internal gravity waves in the absence of mean flow (left), and the Doppler shifted velocities C (n) ± ∓ U d , where we set arbitrarily U d = 0.6 (right). The red lines indicate the limit of the dispersion relation for n → ∞, the crossing is situated at kx = 1/F rU d . We indicate with letters the resonances (crossing points) corresponding to different modes. Note that resonances (b) and (c) come in symmetric pairs. The red circle indicates the proximity of a degenerate crossing point. The single dashed-dotted line corresponds to C (1) + , i.e. a non transported prograde n = 1 mode. Bottom: Same as top right but at F r = 0.2, kz = 18 (left) and F r = 5, kz = 3 (right).

Figure 10 .

 10 Figure 10. Left: eigenvalues in the complex space for the reference case Re = 966, F r = 0.82, kx = 0.96, kz = 5.16 and Sc = ∞ (circles), Sc = 700 (crosses), Sc = 7 (squares), Sc = 1 (diamonds and stars). Right: zoom on the region contoured by the dashed red line in the left diagram, Sc = 7 symbols are reported in red.

Figure 11 .

 11 Figure11. Left: at the top, vertical kinetic energy (w 2 ) 1/2 as a function of time at Re = 966, F r = 0.82 and Sc = 7. The thick line refers to the DNS simulation while the thin line refers to the growth of the most unstable mode as indicated by linear stability analysis. At the bottom, horizontal velocity perturbation u at x = y = 0 as a function of t and z for the same DNS. One observes that a stationary pattern appears close to t = 500. Right: instantaneous 3D map of the buoyancy perturbation once the flow has become unstable (DNS). As predicted by the linear analysis the selected mode is mainly modulated in the vertical direction but still not stream-wise invariant (i.e. kx = 0). Note also that perturbations are concentrated near the boundaries y = ±1.

Figure 12 .

 12 Figure 12. Left: sketch of the experimental shearing device seen from the side (top) and from above (bottom). The two green shaded area correspond to two laser sheets which enlight the mid vertical plane (i.e. y = 0) and the horizontal mid plane (i.e. z = 0). Two cameras allow to image the flow in the enlighted areas. Right: schematic of confined DNS experiments. Two rigid lateral walls entrain the fluid at constant velocity, and two rigid walls confine the flow in the stream-wise direction. Vertical boundary conditions are periodic.

Figure 13 .

 13 Figure 13. Left: top view of the z = 0 plane. On the top we show a superposition of 40 images (i.e. ∼ 1.3 s) as captured by the camera. Only light contrast was exalted to show PIV particle trajectories. On the bottom the velocity field (u, v) as reconstructed via the PIV algorithm. The velocity plot is obtained averaging over 40 PIV fields (i.e. ∼ 1.3 s). Right: horizontal velocity as a function of the cross-stream direction y at four different times Tν /10 (crosses), Tν /5 (squares), Tν /3 (circles) and Tν (triangles). Symbols refer to experimental observations. Each profile corresponds to the time average of the horizontal velocity fields over 0.03Tν . Solid lines refer to the value of the expression (5.1) expected for the two infinite walls problem (increasing time from clear gray to black). The dashed line refers to the asymptotic t = ∞ solution.

Figure 14 .

 14 Figure 14. Top: horizontal velocity perturbation u for a moderate Re experiment (Re = 530, F r = 0.45) in the mid vertical plan y = 0. On the left, u as a function of the time t and the vertical coordinate z (x = 0). On the right, u as a function of the horizontal coordinate x and the time t (z = 0). Bottom left: horizontal velocity perturbation u for the same experiment in the mid horizontal plan z = 0, as a function of time t and y (x = -1). All the quantities are non dimensional. Bottom right: global frequency f box as a function of the imposed shear σ for a collection of experiments with different stratifications N (different symbols). f box was computed as the maximum in the temporal Fourier transform of ū(t), which is the vertical average of u at x = y = 0. In the inset we report the corresponding global period in unit of revolution time of the belt as a function of σ.

Figure 15 .

 15 Figure 15. Top: Horizontal velocity perturbation u at x = y = 0 as a function of time. Colorbar is set to ±15% of the wall speed. Here σ = 1.15 s -1 , F r = 0.82, Re = 969. All quantities are dimensionless. At t = 103.5 the imposed shear switched from 0.34 s -1 to 1.15 s -1. Bottom left: Evolution of the mean horizontal perturbation (u 2 ) 1/2 as a function of time for the same experiment (black line) and for a stable experiment where both F r and Re are diminished by a fraction 1/8 (red line). The inset shows log (u 2 ) 1/2 for the unstable case. At each time u 2 is obtained averaging over a short interval ∼ 4 , taking the square and finally averaging over the vertical direction. Bottom right: Snapshot of the horizontal velocity perturbation u in the plane y = 0, here t ∼ 900 .

Figure 16 .

 16 Figure 16. Mean horizontal perturbation (u 2 ) 1/2 as a function of the imposed shear σ. The value of u 2 is obtained averaging u over the vertical direction z and time. Close and open symbols are assessed to be stable and unstable cases respectively. Squares refer to large gap experiments (L) and circles to narrow gap experiments (N). Different experiments refer to different ratio F r/N as visible in figure 17.

Figure 17 .

 17 Figure 17. Growth rate of the most unstable mode in the space (Re, F r) exactly as in figure 4. Now we superpose the results of experiments (symbols) with the same code as in figure 16. One sees that transition from stable (closed symbols) to unstable (open symbols) cases happen close to the marginal contour Im(ω) = 0. The vertical dashed line corresponds to Re = 230 at which the non stratified flow (i.e. pure water) becomes unstable.

Figure 18 .

 18 Figure18. Solutions of the modified eigenvalue problem for a base flow given by expression (5.1). Here we choose the same Re, F r, and Ui (i.e. the belt velocity at t = 0) as the reference experiment. The eigenvalues problem is solved for the combination (kx, kz) which is the most unstable according to linear analysis and at seven different transient times. Different symbols correspond to different times. The dashed line marks the limit for instability.

Figure 19 .

 19 Figure 19. Left: Spatio-temporal diagram of the perturbation u at the center line x = 0, y = 0 for the reference case Re = 969, F r = 0.82 for confined DNS (top) and the reference experiment (bottom). Spurious bulk oscillations are filtered from the experimental data. Right: perturbation u in the plane y = 0 once the flow has become unstable for confined DNS (top) and the experiment (bottom). The red dashed rectangle indicates the area accessible to experimental measurements.

Figure 20 .

 20 Figure 20. Spatio temporal diagram of the perturbation u at the center line x = 0, y = 0 from DNS. Here Re = 966, F r = ∞ (i.e. unstratified) and the box geometry is the same as in the experiment.

  Figure21. Vertical kinetic energy density (w 2 ) 1/2 as a function of time, for five different DNS experiments. Black lines correspond to supercritical simulations (Re = 966) performed at F r = 0.82 (solid line), F r = ∞ (dashed line), and F r = 0.82 with periodic boundary conditions in the stream-wide direction (dash-dotted line). Red lines refer to subcritical simulations (Re = 629) performed at F r = 0.53 (solid line) and F r = ∞ (dashed line). The thin black line corresponds to the growth of the most unstable mode at Re = 966, F r = 0.82 as predicted by the linear theory. Horizontal dashed lines highlight the saturation level of supercritical unstratified and stratified DNS respectively.

Acknowledgments

This work has been carried out thanks to the support of the A*MIDEX grant (ANR-11-IDEX-0001-02) funded by the French Government "Investissements d'Avenir" program.