Mikael Cugnet

Florian Gallois

Angel Kirchev

Mikaël Cugnet

Denys Dutykh

Understanding self

What is self-discharge?

Self-discharge is a set of processes that decreases the performance of electrochemical power sources without ow of current through an external circuit.

Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion

What is self-discharge?

Self-discharge is a set of processes that decreases the performance of electrochemical power sources without ow of current through an external circuit.

Batteries that are prone to self-discharge exhibit a short shelf life and are thus disfavored in most applications.

Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion

What is self-discharge?

Self-discharge is a set of processes that decreases the performance of electrochemical power sources without ow of current through an external circuit.

Batteries that are prone to self-discharge exhibit a short shelf life and are thus disfavored in most applications. Self-discharge is more critical for primary batteries than secondary batteries:

It may have dramatic consequences for systems that cannot be powered.

It requires excessive charging that may be detrimental to the batteries. It demonstrates that batteries remain active even though they are unused:

∀ δt > 0, di dt ≤ 0, i = 0 ⇒ u (t + δt) < u (t)
It may have dramatic consequences for systems that cannot be powered.

It requires excessive charging that may be detrimental to the batteries. It demonstrates that batteries remain active even though they are unused:

∀ δt > 0, di dt ≤ 0, i = 0 ⇒ u (t + δt) < u (t)
This feature is neglected in most battery models for simplicity! Anodic: Anodic:

H 2 O ⇌ 1 2 O 2 + 2 H + + 2 e -(oxygen evolution) Pb + H 2 O ⇌ PbO + 2 H + + 2 e -(grid corrosion)
H 2 O ⇌ 1 2 O 2 + 2 H + + 2 e -(oxygen evolution) Pb + H 2 O ⇌ PbO + 2 H + + 2 e -(grid corrosion) Introduction Self-discharge of lead-acid cells
Modeling self-discharge of a lead-acid cell Conclusion

Redox shuttle process

Insuciently puried water or leaching impurities from cell components

Impurities, e.g. Fe 2+ /Fe 3+ , lead to discharging of the cell by means of a so-called redox shuttle process:

Positive electrode

PbO2 + H2SO4 + 2 H + + 2 Fe 2+ ⇌ PbSO4 + 2 H2O + 2 Fe 3+
Negative electrode

Pb + H2SO4 + 2 Fe 3+ ⇌ PbSO4 + 2 H + + 2 Fe 2+
This process continues as long as the electrode potentials allow the oxidation and reduction of ions.

The redox shuttle ions themselves are not consumed. Impurities, e.g. Fe 2+ /Fe 3+ , lead to discharging of the cell by means of a so-called redox shuttle process:

Positive electrode

PbO2 + H2SO4 + 2 H + + 2 Fe 2+ ⇌ PbSO4 + 2 H2O + 2 Fe 3+
Negative electrode

Pb + H2SO4 + 2 Fe 3+ ⇌ PbSO4 + 2 H + + 2 Fe 2+
This process continues as long as the electrode potentials allow the oxidation and reduction of ions.

The redox shuttle ions themselves are not consumed. The degradation of the separator must be avoided by proper choice of materials and designs. The degradation of the separator must be avoided by proper choice of materials and designs. Positive electrode as working electrode Positive electrode as working electrode The current density of charge transfer is described by the Butler-Volmer equation:

PbO 2 + H 2 SO 4 + 2 H + + 2 e -⇌ PbSO 4 + 2 H 2 O 2 Ag + H 2 SO 4 ⇌ Ag 2 SO 4 + 2 H + + 2 e -
i = i 0 exp (1 -α) 2 F R T • ∆ϕ -∆ϕ 0 -exp -α 2 F R T • ∆ϕ -∆ϕ 0
with the exchange current density

i 0 = n F k 0 γ ‡ a 1-α H 2 SO 4 ,ref a α H 2 SO 4
and the equilibrium potential at the reference temperature from the Nernst equation: The Arrhenius equation gives the dependence of the rate constant on

∆ϕ 0 Pb (T • ) = µ 0 PbSO 4 -µ 0 Ag 2 SO 4 2 F - R T • 2 F ln a H 2 SO 4 a H 2 SO 4 ,ref
T k 0 = k • 0 exp - ∆ Ḡ ‡ R T
Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion

Working electrode

Impact of temperature on current density and equilibrium voltage

The Arrhenius equation gives the dependence of the rate constant on

T k 0 = k • 0 exp - ∆ Ḡ ‡ R T
The Taylor polynomial of degree 2 for the equilibrium potential according to

T • ∆ϕ 0 (T ) = ∆ϕ 0 (T • ) + (T -T • ) ∂∆ϕ 0 ∂T p (T • ) + 1 2 (T -T • ) 2 ∂ 2 ∆ϕ 0 ∂T 2 p (T • )
allows to estimate the equilibrium potential for any temperature value T ̸ = T • from the entropy and heat capacity at constant pressure (isobaric process) The porous electrodes of lead-acid batteries are made of two phases: a solid phase (Pb at the negative electrode and PbO 2 at the positive electrode), a liquid phase (electrolyte). This approach is simplistic as the electrodes are also made of binders (electrical insulators) and conductive additives, but allows to consider that the whole volume of the electrode is conductive with two dierent types of conduction according to the phases: electronic in the solid phase (exclusively metallic), symbolized by the current i s , ionic in the liquid phase (electrolyte), symbolized by the current i l .

∆ϕ 0 (T ) = ∆ Ḡ 0 (T • ) n F -(T -T • ) ∆ S0 (T • ) n F -(T -T • ) 2 ∆ C 0 p (T • ) Introduction Self-discharge
The porous electrodes of lead-acid batteries are made of two phases: a solid phase (Pb at the negative electrode and PbO 2 at the positive electrode), a liquid phase (electrolyte). This approach is simplistic as the electrodes are also made of binders (electrical insulators) and conductive additives, but allows to consider that the whole volume of the electrode is conductive with two dierent types of conduction according to the phases: electronic in the solid phase (exclusively metallic), symbolized by the current i s , ionic in the liquid phase (electrolyte), symbolized by the current i l .

Electroneutrality (the double layer being negligible in volume)

Charges leaving the porous matrix of the electrode must pass into the solution and reciprocally:

∇ • i s + ∇ • i l = 0
The current per unit volume of the electrode ∇ • i l , transferred from the matrix to the solution, is involved either in faradaic reactions

∇ • i l = A n j=1 i j
where A is the specic interfacial area (m -1 ), n j=1 i j the sum of the current densities of all faradaic processes j involved (A m -2 ) Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion

Working porous electrode Kinetics and electrostatics are two sides of the same coin

The current per unit volume of the electrode ∇ • i l , transferred from the matrix to the solution, is involved either in faradaic reactions or in double-layer charging:

∇ • i l = A n j=1 i j + d dt (A q dl )
where A is the specic interfacial area (m -1 ), n j=1 i j the sum of the current densities of all faradaic processes j involved (A m -2 ) and q dl the surface charge density on the electrode side of the double layer (C m -2 ). In open circuit, no current ows, neither electronically (∇ • i s = 0), nor ionically:

Introduction

∇ • i l = A n j=1 i j + d dt (A q dl ) = 0
Assuming two faradaic processes: discharge of the electrode active material with the primary current i p hydrolysis of water contained in electrolyte with the gassing current i g Reducing self-discharge to sulfuric acid consumption works for the negative:

Negative electrode Anodic: Pb + H 2 SO 4 ⇌ PbSO 4 + 2 H + + 2 e -(discharge) Cathodic: 2 H + + 2 e -⇌ H 2 (hydrogen evolution)
but is not sucient for the positive, producing and consuming water:

Positive electrode Cathodic: PbO 2 + H 2 SO 4 + 2 H + + 2 e -⇌ PbSO 4 + 2 H 2 O (discharge) Anodic: H 2 O ⇌ 1
Self-discharge means sulfuric acid consumption and dilution

Reducing self-discharge to sulfuric acid consumption works for the negative:

Negative electrode Anodic: Pb + H 2 SO 4 ⇌ PbSO 4 + 2 H + + 2 e -(discharge) Cathodic: 2 H + + 2 e -⇌ H 2 (hydrogen evolution)
but is not sucient for the positive, producing and consuming water: A system of 7 ordinary dierential equations (2 for each electrode and 3 for the electrolytes)

Positive electrode Cathodic: PbO 2 + H 2 SO 4 + 2 H + + 2 e -⇌ PbSO 4 + 2 H 2 O (discharge) Anodic: H 2 O ⇌ 1 2 O 2 + 2 H + + 2 e -(oxygen evolution) ∀ i ∈ {H 2 O, H 2 SO 4 } , j ∈ {p, g} , k ∈ {pos, neg} , dn i dt = - j,k s i,j,k N k S k n j,k F i j,k
Electrode k

d∆ϕ k dt = - 1 C dl,k j i j,k dQ k dt = -s p,k N k S k i p,k SOC k = Q k max (Q k ) Electrode k d∆ϕ k dt = - 1 C dl,k j i j,k dQ k dt = -s p,k N k S k i p,k SOC k = Q k max (Q k )
The cell SOC corresponds to the minimum value of each electrode SOC k :

SOC = min max min (SOC k ) , 0 , 1 Introduction Self-discharge of lead-acid cells
Modeling self-discharge of a lead-acid cell Conclusion

Model equations

A system of 7 ordinary dierential equations (2 for each electrode and 3 for the electrolytes)

Electrode k

d∆ϕ k dt = - 1 C dl,k j i j,k dQ k dt = -s p,k N k S k i p,k SOC k = Q k max (Q k )
The cell SOC corresponds to the minimum value of each electrode SOC k :

SOC = min max min (SOC k ) , 0 , 1 Bulk electrolyte dn H 2 O dt = - j,k s H 2 O,j,k N k S k n j,k F i j,k dn H 2 SO 4 dt = - j,k s H 2 SO 4 ,j,k N k S k n j,k F i j,k m = 1 000 n H 2 SO 4 M H 2 O n H 2 O
A system of 7 ordinary dierential equations (2 for each electrode and 3 for the electrolytes)

Electrode k

d∆ϕ k dt = - 1 C dl,k j i j,k dQ k dt = -s p,k N k S k i p,k SOC k = Q k max (Q k )
The cell SOC corresponds to the minimum value of each electrode SOC k :

SOC = min max min (SOC k ) , 0 , 1 Bulk electrolyte dn H 2 O dt = - j,k s H 2 O,j,k N k S k n j,k F i j,k dn H 2 SO 4 dt = - j,k s H 2 SO 4 ,j,k N k S k n j,k F i j,k m = 1 000 n H 2 SO 4 M H 2 O n H 2 O
Reference electrolyte

dm ref dt = D ref S ref (m -m ref )
Positive electrode voltage Negative electrode voltage The positive electrode polarization is much higher than the negative one. Relative errors of density and voltages are very low (±0.3 %).

Diusion of species from the bulk is very slow but non negligible.

The model shows how the acid is diluted by water production.

Electrolyte density Cell voltage

Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion

Model extrapolation

Simulating the impact of a one-year selfdischarge at 25 • C on electrolyte density, cell voltage and SOC

Electrolyte density Cell voltage

Cell SOC Thermodynamics and kinetics make the electrolyte density, cell voltage, and SOC decrease exponentially, showing that the lower the SOC, the lower the selfdischarge phenomenon.

Measurements from the electrolyte reserve located above the electrodes may lead to lower density values than the actual ones inside the pores. (19151998)

  of lead-acid cells Modeling self-discharge of a lead-acid cell Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Introduction What is self-discharge? Why self-discharge? Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion

  Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Redox shuttle process Insuciently puried water or leaching impurities from cell components

  Number and dimensions of positive and negative electrodes Number and dimensions of separators Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Self-discharge in practice Tear-down of an individual electrochemical cell Number and dimensions of positive and negative electrodes Number and dimensions of separators Weight of each electrode (still soaked with electrolyte) Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Self-discharge in practice Tear-down of an individual electrochemical cell Number and dimensions of positive and negative electrodes Number and dimensions of separators Weight of each electrode (still soaked with electrolyte) Weight of each electrode (after drying to remove electrolyte) Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Self-discharge in practice Tear-down of an individual electrochemical cell Number and dimensions of positive and negative electrodes Number and dimensions of separators Weight of each electrode (still soaked with electrolyte) Weight of each electrode (after drying to remove electrolyte) Weight of each grid (after removing active materials) Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Reference electrode Study of each electrochemical reaction with the silver-silver sulfate reference electrode To measure the potential of an individual working electrode, you must combine it with a reference electrode to form a galvanic cell. Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Reference electrode Study of each electrochemical reaction with the silver-silver sulfate reference electrode To measure the potential of an individual working electrode, you must combine it with a reference electrode to form a galvanic cell.

  Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Reference electrode Study of each electrochemical reaction with the silver-silver sulfate reference electrode To measure the potential of an individual working electrode, you must combine it with a reference electrode to form a galvanic cell.

  IntroductionSelf-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell ConclusionWorking electrodeImpact of temperature on current density and equilibrium voltage

  of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Working porous electrode Electroneutrality and conservation of charge The porous electrodes of lead-acid batteries are made of two phases: a solid phase (Pb at the negative electrode and PbO 2 at the positive electrode), a liquid phase (electrolyte). Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Working porous electrode Electroneutrality and conservation of charge

  Self-discharge of lead-acid cellsModeling self-discharge of a lead-means no charge transferredIn open circuit, no current ows, neither electronically (∇ • i s = 0), nor ionically:

  and a constant double layer capacity per unit area of the electrode dC dl dt = 0 d∆ϕ dt = -i p + i g C dl Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Model equations Self-discharge means sulfuric acid consumption and dilution Reducing self-discharge to sulfuric acid consumption works for the negative: Negative electrode Anodic: Pb + H 2 SO 4 ⇌ PbSO 4 + 2 H + + 2 e -(discharge) Cathodic: 2 H + + 2 e -⇌ H 2 (hydrogen evolution) Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Model equations Self-discharge means sulfuric acid consumption and dilution

  Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Model validation with experimental data Electrode voltage and SOC during 49 days at 25 • C, 34 days at 35 • C, and 16 days at 45 • C Positive electrode voltage Negative electrode voltage State of charge The positive electrode polarization is much higher than the negative one. The positive electrode voltage decreases faster than the negative one. Bulk electrolyte density Bulk and ref molalities Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Model validation with experimental data Bulk and reference electrolytes during 49 days at 25 • C, 34 days at 35 • C, and 16 days at 45 • C Relative errors of density and voltages are very low (±0.3 %). Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Model validation with experimental data Bulk and reference electrolytes during 49 days at 25 • C, 34 days at 35 • C, and 16 days at 45 • C Relative errors of density and voltages are very low (±0.3 %).Diusion of species from the bulk is very slow but non negligible.

  Self-discharge of lead-acid cellsModeling self-discharge of a lead-acid cell ConclusionModel improvementsElectrolyte density measurement issues, long-term selfdischarge at 40 • C, and sealed technologies Measurements from the electrolyte reserve located above the electrodes may lead to lower density values than the actual ones inside the pores.Measurements with a densimeter have to be done in accordance with the allowed temperature range (measurements above 40 • C may be inaccurate). Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion Model improvements Electrolyte density measurement issues, long-term selfdischarge at 40 • C, and sealed technologies Measurements from the electrolyte reserve located above the electrodes may lead to lower density values than the actual ones inside the pores. Measurements with a densimeter have to be done in accordance with the allowed temperature range (measurements above 40 • C may be inaccurate). A long-term self-discharge at 40 • C has been started recently to: experimentally validate the model beyond one-month, potentially add the impact of corrosion and loss of water by evaporation, if necessary. Electrolyte density measurement issues, long-term selfdischarge at 40 • C, and sealed technologies Measurements from the electrolyte reserve located above the electrodes may lead to lower density values than the actual ones inside the pores. Measurements with a densimeter have to be done in accordance with the allowed temperature range (measurements above 40 • C may be inaccurate). A long-term self-discharge at 40 • C has been started recently to: experimentally validate the model beyond one-month, potentially add the impact of corrosion and loss of water by evaporation, if necessary. The model should be applied on sealed lead-acid cells as well. Introduction Self-discharge of lead-acid cells Modeling self-discharge of a lead-acid cell Conclusion What we have learned What should be done next If you don't know precisely what happens inside one cell, when there is no current going through it, do you sincerely believe you could drop that step and jump directly to the modeling of battery (dis)charge with more success? If you don't know precisely what happens inside one cell, when there is no current going through it, do you sincerely believe you could drop that step and jump directly to the modeling of battery (dis)charge with more success? "It is better to do the right problem the wrong way than the wrong problem the right way." R. Hamming Richard Hamming

  Impact of the reference electrode on the Butler-Volmer model of the negative electrode

PbO 2 + H 2 SO 4 + 2 H + + 2 e -⇌ PbSO 4 + 2 H 2 O 2 Ag + H 2 SO 4 ⇌ Ag 2 SO 4 + 2 H + + 2 e - Negative electrode as working electrode Pb + H 2 SO 4 ⇌ PbSO 4 + 2 H + + 2 e - Ag 2 SO 4 + 2 H + + 2 e -⇌ 2 Ag + H 2 SO 4

  Electrode voltage and SOC during 49 days at 25 • C, 34 days at 35 • C, and 16 days at 45 • C Electrode voltage and SOC during 49 days at 25 • C, 34 days at 35 • C, and 16 days at 45 • C

	Introduction	Self-discharge of lead-acid cells Self-discharge of lead-acid cells	Modeling self-discharge of a lead-acid cell Modeling self-discharge of a lead-acid cell	Conclusion Conclusion
	Model validation with experimental data Model validation with experimental data
	Positive electrode voltage Positive electrode voltage	Negative electrode voltage Negative electrode voltage	State of charge State of charge

n F T

O 2 + 2 H + + 2 e -(oxygen evolution)

Self-discharge has an obvious impact on operational performance and safety.

Self-discharge in practice Tear-down of an individual electrochemical cell Number and dimensions of positive and negative electrodes Number and dimensions of separators Weight of each electrode (still soaked with electrolyte) Weight of each electrode (after drying to remove electrolyte) Weight of each grid (after removing active materials) Estimation of the total volume of electrolyte in the cell Acknowledgments Acknowledgments Co-workers on this subject Acknowledgments Co-workers on this subject Colleagues for their fruitfull discussion and managers for support Acknowledgments Co-workers on this subject Colleagues for their fruitfull discussion and managers for support Craig Wohlers for giving me the opportunity to present my work and for organizing this great event of the battery community. Acknowledgments Co-workers on this subject Colleagues for their fruitfull discussion and managers for support Craig Wohlers for giving me the opportunity to present my work and for organizing this great event of the battery community. My wife and my three kids

Kinetics and electrostatics are two sides of the same coin

The current per unit volume of the electrode ∇ • i l , transferred from the matrix to the solution, is involved either in faradaic reactions or in double-layer charging:

where A is the specic interfacial area (m -1 ), n j=1 i j the sum of the current densities of all faradaic processes j involved (A m -2 ) and q dl the surface charge density on the electrode side of the double layer (C m -2 ).

where C dl is the double layer capacity per unit area of the electrode (F m -2 ). The positive electrode polarization is much higher than the negative one.

The positive electrode voltage decreases faster than the negative one.

The positive electrode SOC decreases slower than the negative one.


P Le Quéré
email: patrick.le-quere@limsi.fr

Natural convection in air-filled differentially heated isoflux cavities: scalings and transition to unsteadiness, a long story made short

Keywords: natural convection, isoflux cavity, instability, unsteadiness, time integration, eigenmodes computation Nomenclature

Natural convection in air-filled differentially heated isoflux cavities: Scalings and transition to unsteadiness, a long story made short

(x, z) Horizontal and vertical coordinates β

Introduction

Buoyancy-induced flows abound in nature or in industrial configurations and the study of their many specific configurations has been and is still the object of intense research. Among the many configurations, the so-called isothermal differentially heated cavity has become one of the most classical problems of the heat transfer literature. The determination of the flow structure, corresponding heat transfer, loss of stability, transition to turbulence have been investigated for more than 60 years now. One specific problem which has resisted years of investigation is the transition to unsteadiness in an air filled cavity differentially heated (and cooled) at uniform heat flux at its vertical boundaries. This is the subject of the present paper.

After a historical journey in the field, I will first address the scaling issue, showing that it is not a good idea to perform the time integration of the unsteady equations in the usual way, that is starting from a value obtained at a lower value of Ra as initial condition, unless one assumes a scaling that removes the dependence of the stratification on the Rayleigh number. Second I will determine the critical value of the Rayleigh number in a square cavity through both time integration and computation of the leading eigenvalues of the spectrum of steady solution. I will discuss the physics of various groups of eigenmodes, and perform a grid independence test that shows the requirements needed in terms of spatial resolution. I will also discuss the spatial structure of the most unstable eigenmode and fluctuating solution. I will then determine upper and lower bounds of the critical Rayleigh number for values of the aspect ratio in the range 0.2 to 8 and show that the traveling waves in the boundary layers have similar characteristics than those resulting from blending the linear stability analysis of the buoyancy layer with the analytical solution. I will finally show how the leading eigenvalues evolve with aspect ratio.

A historical perspective

For one of the simplest configurations i.e. the flow adjacent to a heated vertical plate in an infinite ambient medium, literature abounds with analysis of analytical solutions for various types of thermal boundary conditions solutions on the plate or in the ambient medium (see e.g. [START_REF] Gebhart | Buoyancy-induced flows and transport[END_REF]).

Whereas early studies considered the case of uniform wall temperature, the case of uniform heat flux was apparently first considered in [START_REF] Sparrow | Laminar free convection from a vertical plate with uniform surface heat flux[END_REF]. Let us first recall the basic scalings of the isoflux natural convection boundary layer. Consider a fluid of kinematic viscosity ν, thermal conductivity λ, and thermal diffusivity κ along a vertical plate of height H heated at a uniform flux φ in an isothermal environment at T ∞ .

The governing equations under both the boundary layer and Boussinesq approximations read :

∂u ∂x + ∂w ∂z = 0 (1)
u ∂w ∂x + w ∂w ∂z = ν ∂ 2 w ∂x 2 + g β (T -T ∞ ) (2)
u ∂T ∂x + w ∂T ∂z = κ ∂ 2 T ∂x 2 (3)
where u and w are the velocity components in horizontal and vertical directions, respectively, and T the temperature. Classical scaling analysis [START_REF] Bejan | Convection Heat Transfer[END_REF] results in the following characteristic scales for the velocity magnitude W and boundary layer thickness δ T :

δ T ∝ HRa -1/5 (4) W ∝ κ H Ra 2/5 ( 5
)
where Ra is the flux Rayleigh number (Ra = gβφH 4 νκλ ). The temperature difference across the boundary layer scales as

∆T ∝ φH λ Ra -1/5 (6)
and the entrainment velocity U like κ H Ra 1/5 . In terms of local vertical height z, the thermal boundary layer thickness δ T thus increases like z 1/5 , the vertical velocity like z 3/5 , ∆T like z 1/5 and the entrainment velocity decreases like z -1/5 . These evolutions can be compared with those characterizing the flow along an isothermal plate for which the boundary layer thickness δ T increases like z 1/4 , the vertical velocity like z 1/2 , and the entrainment velocity decreases like z -1/4 . Compared to the isothermal case, the vertical velocity in the isoflux boundary layer thus increases more rapidly and the boundary layer thickness less rapidly with altitude, which is due to the constant input of buoyancy flux.

A sub category of thermally buoyancy-induced flows are those which occur in enclosures when the walls are maintained at different temperatures. We will concentrate here on those induced by horizontal thermal gradients, idealized as the so-called differentially heated cavity, a configuration consisting of a rectangular cavity whose vertical opposite walls are maintained at different constant temperatures, the horizontal ones being either adiabatic or perfectly conducting. The first attempt to elucidate its flow structure and related heat transfer is presumably due to Batchelor [START_REF] Batchelor | Heat transfer by free convection across a closed cavity between vertical boundaries at different temperatures[END_REF]. As precisely described in [START_REF] Ostrach | Natural convection in enclosures[END_REF], Batchelor hypothesized that the core of the cavity would be isothermal, from which he derived that the flow in the core would be rotating like a solid body. Later experiments ( [START_REF] Eckert | Natural convection in an air layer enclosed between two vertical plates with different temperatures[END_REF], [START_REF] Elder | Laminar free convection in a vertical slot[END_REF]) showed that, for sufficiently large values of the Rayleigh number, the core was not isothermal but thermally stratified instead. This led to the identification of the boundary layer regime whose first analytical investigation was performed by Gill [START_REF] Gill | The boundary-layer regime for convection in a rectangular cavity[END_REF] for large values of the Prandtl number. The free constant that appears in the general solution was set requiring that the vertical velocity be zero at both horizontal walls. This was revisited by Quon some time later [START_REF] Quon | Free convection in an enclosure revisited[END_REF] and improved by Bejan [START_REF] Bejan | Note on Gill's solution for free convection in a vertical enclosure[END_REF] who required that the constant be determined requiring the vertical heat flux be zero, leading to a much better agreement with the numerical solutions available.

Kimura and Bejan [START_REF] Kimura | The boundary layer natural convection regime in a rectangular cavity with uniform heat flux from the side[END_REF] investigated the isoflux differentially heated cavity along the lines of what Gill had done for the isothermal cavity [START_REF] Gill | The boundary-layer regime for convection in a rectangular cavity[END_REF] using the alternative proposed by Bejan [START_REF] Bejan | Note on Gill's solution for free convection in a vertical enclosure[END_REF] to relate the vertical stratification in the core to the cavity characteristics. The general solution depends on two parameters p and q, which are determined from symmetry properties of the solution with respect to the cavity center. They found that the temperature of the walls rises linearly with height at the same rate as the temperature in the core which is motionless. In the reference quantities [START_REF] Batchelor | Heat transfer by free convection across a closed cavity between vertical boundaries at different temperatures[END_REF][START_REF] Ostrach | Natural convection in enclosures[END_REF][START_REF] Eckert | Natural convection in an air layer enclosed between two vertical plates with different temperatures[END_REF] the vertical velocity and temperature profiles across the boundary layer read ;

w = 32 q 3 exp(- q 4 x) sin( q 4 x) θ = 4 q exp(- q 4 x) cos( q 4 x) (7)
where the key parameter q is linked to the cavity characteristics :

q 9 = 8192 H L Ra 1/5 H (8)
This means that the dimensional vertical velocity behaves as Ra2/5 × Ra -1/15 i.e. like Ra 1/3 and the boundary layer thickness scales like Ra -1/5 × Ra -1/45 i.e. like Ra -2/9 . Likewise the temperature difference across the boundary layer scales like φH λ Ra -2/9 . This is precisely the velocity scaling proposed by Kouroudis et al [START_REF] Kouroudis | Direct numerical simulation of natural convection in a square cavity with uniform heat fluxes at the vertical sides: Flow structure and transition[END_REF] 2 , who derived their scales using directly the thermal balance between upward convection and downward conduction in the scaling analysis. This results in a slightly different scaling for temperature φH λ Ra -1/9 . They also provided scalings for the horizontal velocity and thickness of the horizontal boundary layer and the top and bottom walls ( U ∝ W Ra -1/18 and δ ∝ HRa -1/6 , respectively). They also addressed the transition to unsteadiness to which we will come back later on.

The flow structure in the isoflux cavity is thus made of two boundary layers whose thickness is independent of height, a noticeable feature that is due to the fact that the equilibrium temperature at the walls increases linearly with height as does the stratification in the core, both at the same rate. A detailed comparison between the isothermal and isoflux cavity for values of the aspect ratio ranging from 1/8 to 8 and two Prandtl values corresponding to air and water was performed by [START_REF] Turan | Influences of boundary conditions on laminar natural convection in rectangular enclosures with differentially heated side walls[END_REF] with particular emphasis on heat transfer. As shown by Prandtl [START_REF] Prandtl | Essentials of Fluid Dynamics[END_REF] in his study of katabatic winds, linear and identical variation of the wall and ambient temperatures with height is indeed a (sufficient) condition for the existence of a one-dimensional solution of the Navier-Stokes equations, named the buoyancy layer which, in boundary layer units, reads :

w(x) = exp(-x) sin(x) θ(x) = exp(-x) cos(x) (9)
The linear stability of this solution was studied by Gill and Davey [START_REF] Gill | Instabilities of a buoyancy-driven system[END_REF] who found a critical Reynolds number of approximately 100 for a fluid of Prandtl number equal to that of air (0.72). This stability analysis was further extended by Iyer [START_REF] Iyer | Instabilities in buoyancy-driven boundary-layer flows in a stably stratifed medium[END_REF] to the case of boundary layers at an angle to the vertical. Iyer and Kelly [START_REF] Iyer | Supercritical solutions for the buoyancy boundary layer[END_REF] investigated the nature of the instability, carrying out non-linear expansions of the 2D solution at second order. They found only supercritical solutions to occur. This analysis was also carried out by True and Bruun Nielsen [START_REF] True | On bifurcation of stable periodic flows in an Ekman layer and in a convection boundary layer[END_REF], who found that the instability was supercritical for values of the Prandtl number smaller than 9.9 and subcritical for larger Prandtl values.

Desrayaud [START_REF] Desrayaud | Stability of the flow near a heat-flux plate and comparison with numerical simulation in a square cavity[END_REF] extended the linear stability analysis to the case of homogeneous Neumann boundary condition for the temperature fluctuation. He showed that changing this boundary condition has a substantial effect on the critical parameters for small values of the Prandtl number. In particular for a Prandtl number of 0.71, the critical Reynolds number drops from 101 to 44, corresponding to a critical wave number and wave speed of 0.387 and 0.319 respectively. Mc Bain et al [START_REF] Mcbain | Instability of the buoyancy layer on an evenly heated vertical wall[END_REF] further extended the stability analysis which they corroborated through direct numerical simulations of the buoyancy layer along a vertical plate. They used a finite difference algorithm assuming periodicity in the streamwise direction for two Prandtl numbers of 7 and 0.1 and Reynolds numbers of 9 and 135, respectively. They observed the development of traveling waves in good agreement with the linear stability results and confirmed the supercritical nature of the bifurcation.

A related configuration is the flow between two vertical walls of infinite vertical extension whose temperatures increase linearly with height. In this configuration, the governing equations admit a family of 1D analytical solutions that depend on the ratio of the imposed stratification to the temperature difference across the slot, known as the stratification parameter. For sufficiently large value of this parameter, and assuming zero net vertical flux, the solution reduces to two buoyancy layers, one upward and one downward, whose expression was given by Elder [START_REF] Elder | Laminar free convection in a vertical slot[END_REF]. The linear stability of the whole family of solutions was investigated by Bergholz [START_REF] Bergholz | Instability of steady natural convection in a vertical fluid layer[END_REF]. Weakly nonlinear effects as well as numerical simulations for large supercritical values of Ra showing the development of soliton-like structures were investigated in [START_REF] Podvin | Nonlinear dynamics between two differentially heated vertical plates in the presence of stratification[END_REF].

Owing to its well posedness natural convection flow in differentially heated cavities has become one of the most popular test cases for testing the efficiency and accuracy of numerical methods to solve the incompressible Navier-Stokes equations. It was proposed as a benchmark problem in the 1980's [START_REF] Vahl | Natural convection in a square cavity: a comparison exercise[END_REF] for the case of air (P r = 0.71) for values of the Rayleigh number limited to 10 6 . Using Richardson extrapolation G. de Vahl Davis provided the first accurate solutions for values of Rayleigh number in the range 10 3 -10 6 [START_REF] Vahl | Natural convection of air in a square cavity : a benchmark numerical solution[END_REF], although his solution for Ra = 10 6 was later found to be a few percents off ( [START_REF] Quéré | Computation of natural convection in twodimensional cavities with Chebyshev polynomials[END_REF][START_REF] Hortmann | Finite volume multigrid prediction of laminar natural convection: Bench-mark solutions[END_REF]). Continuous improvements in numerical algorithms, in particular Chebyshev spectral methods which proved ideally suited for this configuration ( [START_REF] Haldenwang | Unsteady numerical simulation by Chebyshev spectral methods of natural convection at high Rayleigh number[END_REF], [START_REF] Quéré | Computation of natural convection in twodimensional cavities with Chebyshev polynomials[END_REF]) and available computer resources have increased the benchmark solutions up to values of Rayleigh number of 10 8 [START_REF] Chen | Accurate solutions to the square thermally driven cavity at high Rayleigh number[END_REF] close to the end of existence of stable steady state solutions. The accuracy of the solution for Ra = 10 8 has been recognized many times since (see e.g. [START_REF] Janssen | Influence of Prandtl number on instability mechanisms and transition in a differentially heated square cavity[END_REF], [START_REF] Nobile | Simulation of time-dependent flow in cavities with the additive-correction multigrid method, part II: applications[END_REF], [START_REF] Gjesdal | Spectral element benchmark simulations of natural convection in two-dimensional cavities[END_REF], [START_REF] Contrino | Lattice-Boltzmann simulations of the thermally driven 2D square cavity at high Rayleigh numbers[END_REF])

These solutions have allowed to study the flow structure in detail, and in particular that the boundary layer regime in an air-filled cavity is characterized by a quasi linearly stratified core with a vertical gradient very close to 1 in units of ∆T H (see eg [START_REF] Henkes | Scaling of the laminar natural convection flow in a heated square cavity[END_REF], [START_REF] Ravi | On the high-Rayleigh-number structure of steady laminar natural-convection flow in a square enclosure[END_REF]). One notable feature is that this stratified core is thus able to support internal gravity waves, as was first shown by Imberger and Patterson [START_REF] Patterson | On the existence of an oscillatory approach to steady natural convection in cavities[END_REF] in their study of the transients for the flow establishment following a sudden imposition of the temperature difference across of the cavity. They showed that these waves originate from the intrusion layers emanating from the boundary layers piling-up against the opposite walls and that it is the time needed to damp these waves that characterizes the time to achieve steady state.

Predicting the transition to unsteadiness in differentially heated cavities and the subsequent route to chaos has drawn considerable interest and efforts in the mid 80's and early 90's ( [START_REF] Paolucci | Transition to chaos in a differentially heated vertical cavity[END_REF], [START_REF] Paolucci | Direct simulation of two dimensional turbulent natural transition in an enclosed cavity[END_REF], [START_REF] Henkes | On the stability of the natural convection flow in a square cavity heated from the side[END_REF], [START_REF] Quéré | From onset of unsteadiness to chaos in a differentially heated square cavity[END_REF], [START_REF] Janssen | Influence of Prandtl number on instability mechanisms and transition in a differentially heated square cavity[END_REF], [START_REF] Xin | Direct numerical simulations of two-dimensional chaotic natural convection in a differentially heated cavity of aspect ratio 4[END_REF]). Methodologies to determine the transition to unsteadiness have proceeded along two main streams initially quite closely related to the type of numerical approximation and methods to solve the governing equations. The use of finite element approximation was generally related to solving directly the steady state equations using inversion of the discretized equations coupling directly the momentum and continuity equations. The availability of the jacobian in matrix form thus allowed for computing its spectrum directly providing the leading eigenvalues and eigenvectors characterizing the stability of the solution. One of the early works along this line was carried out by Winters [START_REF] Winters | Hopf bifurcation in the double glazing problem with conducting boundaries[END_REF] who determined the loss of stability of the square differentially heated cavity with perfectly conducting top and bottom walls. He computed a critical Rayleigh number close to 2.1×10 6 , in good agreement with the value determined experimentally by Briggs and Jones [START_REF] Briggs | Two-dimensional periodic natural convection in a rectangular enclosure of aspect ratio one[END_REF]. As problems got harder and necessitated the use of larger spatial resolutions, direct methods faced numerical difficulties in particular memory requirements which grow very rapidly. This promoted the development of matrix-free methods either to solve the steady state equations either stable or unstable through Newton's iterations or to get the leading eigenvalues responsible for the loss of stability of the corresponding solutions ( [START_REF] Edwards | Krylov methods for the incompressible Navier-Stokes equations[END_REF], [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF], [START_REF] Loiseau | Time-stepping and Krylov methods for large-scale instability problems[END_REF]).

The second methodology consists in integrating the time dependent equations and observing the asymptotic nature of the solution, either steady or unsteady. This methodology places additional constraints on the numerical algorithm which has to be accurate enough from a time-space view point. This was generally done by progressively increasing the Rayleigh number starting from an initial condition corresponding to the steady solution found for the immediately smaller value of the Rayleigh number. Its efficiency relies on some a priori knowledge of the time needed to reach to new asymptotic mean solution, and from this standpoint the fact that in the isothermal differentially heated cavity the thermal stratification remains quasi-constant as the Rayleigh number increases results in the fact that the transient reduces to a "simple" adjustment of the thickness of the boundary layers in a first stage, followed by the time needed to damp the internal waves generated during the transient. Using this approach Le Quéré and Alziary de Roquefort [START_REF] Gelfgat | Transition to unsteady natural convection of air in vertical differentially heated cavities: influence of thermal boundary conditions on the horizontal walls[END_REF], using a second order time stepping method coupled to Chebyshev spatial approximation, were able to determine bounds on the Rayleigh number corresponding to the transition to unsteadiness for both adiabatic and perfectly conducting horizontal walls for aspect ratios in the range 1 to 10 for perfectly conducting horizontal walls and 2 to 10 for adiabatic walls. It was shown that depending on the aspect ratio and thermal boundary conditions on the top and bottom walls, three different instability mechanisms are responsible for the transition to unsteadiness. For large values of the aspect ratio (say larger than 4), the instability is of boundary layer type independently of the thermal boundary conditions on the horizontal walls. For smaller values of the aspect ratio and adiabatic top and bottom walls the instability originates at the base of the detached region that develops along the ceiling after the boundary layer impinges on the ceiling and starts to flow horizontally. The crossover of the two mechanisms correspond to a co-dimension 2 bifurcation. Using a Galerkin Chebyshev approximation with a resolution of (32 × 50) the critical value of the aspect ratio corresponding to the crossover was later determined by Yahata [START_REF] Yahata | Stability analysis of natural convection in vertical cavities with lateral heating[END_REF] to be equal to 3.41 3 . Xin and Le Quéré [START_REF] Xin | Natural convection flows in air-filled differentially heated cavities with adiabatic horizontal walls[END_REF] showed that for values of the aspect ratio in the range 3-4 the branch of steady solutions has two turning points, which they attributed to a sudden adjustment of the size of the detached structure to the cavity width. For small values of the aspect ratio and perfectly conducting walls, the instability is of thermal type, as the fluid is more and more heated from below as it is drawn into the ascending boundary layer. It was also found that different instability modes with different number of structures co-exist, with different symmetry properties, and that the solution can saturate on either one of these modes. Hysteresis effects exist when increasing or decreasing the Rayleigh number ( [START_REF] Jones | Periodic two-dimensional cavity flow: Effect of linear horizontal thermal boundary condition[END_REF], [START_REF] Quéré | Onset of unsteadiness routes to chaos and simulations of chaotic flows in cavities heated from the side: a review of present status[END_REF], [START_REF] Xin | Natural convection flows in air-filled differentially heated cavities with adiabatic horizontal walls[END_REF]).

The results were later extended to square cavities with adiabatic horizontal walls ([50], [START_REF] Quéré | From onset of unsteadiness to chaos in a differentially heated square cavity[END_REF]) and the critical Rayleigh number corresponding to the transition to unsteadiness was found to occur at 1.82×10 8 . It was also found that the instability mechanism is linked to the base of the detached region and that the first instability mode does not possess the centrosymmetry property of the base flow. The accuracy of this critical value has been confirmed many times since, either through direct numerical integration ( [START_REF] Gjesdal | Spectral element benchmark simulations of natural convection in two-dimensional cavities[END_REF], [START_REF] Norris | Prandtl number dependence of the onset of instability for differentially heated cavity flow[END_REF]) or through direct search of the most unstable critical eigenvalue ( [START_REF] Yahata | Stability analysis of natural convection in vertical cavities with lateral heating[END_REF], [START_REF] Xin | Natural convection flows in air-filled differentially heated cavities with adiabatic horizontal walls[END_REF], [START_REF] Yu | Stability of convective flows in cavities: solution of benchmark problems by a low-order finite volume method[END_REF]). Another noteworthy result reported in [START_REF] Quéré | From onset of unsteadiness to chaos in a differentially heated square cavity[END_REF] was that the fluctuating temperature field displayed inclined stripes in the stratified core which were attributed to internal gravity waves due to the fact that the primary instability mechanism is characterized by a low frequency, smaller that the cut-off Brunt-Väisälä frequency related to the stratification in the core. The stripes inclination's angle was found to closely match the analytical relationship between the angle and the ratio of the main frequency to the Brunt-Väisälä frequency (see e.g. [START_REF] Tritton | Physical Fluid Dynamics[END_REF], pp. 208-211). This was later confirmed by Xin and Le Quéré [START_REF] Xin | Natural convection flows in air-filled differentially heated cavities with adiabatic horizontal walls[END_REF] who showed that the internal waves were intrinsically part of the spatial structure of the most unstable eigenmode, which is consistent with the fact that the transition corresponds to a Hopf bifurcation (see also [START_REF] Lappa | Thermal Convection: Patterns, Evolution and Stability[END_REF], pp. 238-239). Critical values have progressively become available for a wide range of aspect ratios and Prandtl numbers (see e.g. [START_REF] Roux | Numerical simulation of oscillatory convection in low-Pr fluids : a GAMM Workshop[END_REF], [START_REF] Henry | Two-and three-dimensional numerical simulations of the transition to oscillatory convection in low-Prandtl-number fluids[END_REF], [START_REF] Burroughs | Linear stability of flow in a differentially heated cavity via large scale eigenvalue calculations[END_REF], [START_REF] Yu | Stability of convective flows in cavities: solution of benchmark problems by a low-order finite volume method[END_REF], [START_REF] Christon | Computational predictability of timedependent natural convection flows in enclosures (including a benchmark solution)[END_REF], [START_REF] Xin | Natural convection flows in air-filled differentially heated cavities with adiabatic horizontal walls[END_REF]).

Let us now come back to the isoflux cavity after this seemingly long, but hopefully not totally useless, digression. On the basis of the results provided by Kimura and Bejan [START_REF] Kimura | The boundary layer natural convection regime in a rectangular cavity with uniform heat flux from the side[END_REF] and of the stability analysis ( [START_REF] Desrayaud | Stability of the flow near a heat-flux plate and comparison with numerical simulation in a square cavity[END_REF], [START_REF] Mcbain | Instability of the buoyancy layer on an evenly heated vertical wall[END_REF]), one can derive the value of the critical Rayleigh number which would correspond to the loss of stability of the boundary layer and hence presumably correspond to the transition to unsteadiness in a differentially heated isoflux cavity. This seems particularly relevant for this configuration since, as the boundary layer is of constant thickness, one could legitimately expect that the transposition of the stability results for the buoyancy layer to the cavity would be particularly meaningful. In order to rescale the velocity and boundary layer thickness to arrive at the expressions ( 9) the reference velocity should therefore be taken as 32 q 3 κ H Ra 2/5 and the scaling for the boundary layer thickness as 4 q HRa -1/5 . Using relation [START_REF] Gill | The boundary-layer regime for convection in a rectangular cavity[END_REF], one arrives at the relationship between the Reynolds number characterizing the buoyancy layer and the cavity characteristics, first derived by Desrayaud [START_REF] Desrayaud | Stability of the flow near a heat-flux plate and comparison with numerical simulation in a square cavity[END_REF],

Ra = (Re × P r) 9 2048 ( H L ) 4 . ( 10
)
Assuming a critical value of Re c = 44 [START_REF] Desrayaud | Stability of the flow near a heat-flux plate and comparison with numerical simulation in a square cavity[END_REF] results in a critical Rayleigh number approximately equal to 1.4×10 10 a very large value compared to that found in a cavity with isothermal walls. This large value places very severe constraints on the numerical requirements needed to address the transition to unsteadiness. First of all, at the corresponding critical value of Ra, the boundary layer thickness is less than one hundredth of the cavity width. From [START_REF] Desrayaud | Stability of the flow near a heat-flux plate and comparison with numerical simulation in a square cavity[END_REF], the wavelength in boundary layer units is equal to 2π 0.387 ≈ 16.24 which corresponds to a wavelength equal to 0.13 × H, meaning that any numerical algorithm should have at least several hundreds of grid points in each direction.

Second and may be more importantly, contrarily to the isothermal cavity for which the stratification in the core does not change as the Rayleigh number is increased, this is not the case for the isoflux cavity. It thus turns out that, integrating the equations from a previous value of Ra, the time needed to reach the new asymptotic solution is very long, on the order of the diffusive time scale on the cavity height ie H 2 κ , that is Ra 2/5 time units in the natural scales [START_REF] Batchelor | Heat transfer by free convection across a closed cavity between vertical boundaries at different temperatures[END_REF][START_REF] Ostrach | Natural convection in enclosures[END_REF], all the more so due to the fact that there is no entrainment in or ejection from the boundary layer to help speed up this diffusive process. Using time steps smaller than the CFL limit for stability and accuracy reasons, one thus needs several hundreds of thousands of time steps, not to say millions, to reach the asymptotic solution.

Several hundreds of grid points in each direction, hundreds of thousands of time steps placed this configuration in the category of hard problems in the 90's. As reported in [START_REF] Quéré | Onset of unsteadiness routes to chaos and simulations of chaotic flows in cavities heated from the side: a review of present status[END_REF], early attempts to predict the transition to unsteadiness in the isoflux cavity have proven quite inconclusive. Integration of the governing equations in the range 1 to 2×10 10 did not show a hint of transition to unsteadiness casting some doubts on the pertinence of the analogy. Increasing progressively the value of Ra, Le Quéré [START_REF] Quéré | Onset of unsteadiness routes to chaos and simulations of chaotic flows in cavities heated from the side: a review of present status[END_REF] was eventually able to approximately locate the transition to take place in between 2×10 12 and 4×10 12 , a value more than two orders of magnitude larger that the alleged critical value. This order of magnitude corresponding to the transition to unsteadiness was recently confirmed by Kouroudis et al [START_REF] Kouroudis | Direct numerical simulation of natural convection in a square cavity with uniform heat fluxes at the vertical sides: Flow structure and transition[END_REF] who performed time integrations of the unsteady equations with a finite element approximation with biquadratic elements for velocity and temperature and bilinear approximation for pressure. They used spatial resolution of 301 × 601 up to 401 × 601 nodes for the higher Ra number of 4 × 10 14 and a Prandtl number of 1. They report that they found the transition to occur at a critical Ra around 3 × 10 12 . They also checked that the scaling they have proposed is valid over the entire range of Ra numbers investigated, that is up to 4 × 10 14 , two orders of magnitude above the onset of unsteadiness.

The large difference between between the critical values resulting from linear stability theory of the boundary layer solution and the bifurcation to unsteadiness in the closed cavity, either isothermal or isoflux, has been logically ascribed to the difference between the concepts of convective and absolute instability [START_REF] Quéré | Onset of unsteadiness routes to chaos and simulations of chaotic flows in cavities heated from the side: a review of present status[END_REF], [START_REF] Armfield | A direct boundary-layer stability analysis of steady-state cavity convection flow[END_REF]. It has motivated several studies of wave properties [START_REF] Armfield | Wave properties of natural-convection boundary layers[END_REF] or convective-absolute stability analysis of various similarity solutions of boundary layers adjacent to vertical walls, [START_REF] Krizhevsky | Convective and absolute instabilities of a buoyancy-induced flow in a thermally stratified medium[END_REF], [START_REF] Tao | Spatio-temporal instability of natural convection boundary layer in thermally stratified medium[END_REF], [START_REF] Tao | Absolute and convective instabilities of natural convection flow in boundary-layer regime[END_REF], [START_REF] Tao | Nonlinear global instability in buoyancy-driven boundary-layer flows[END_REF]. Although Tao [START_REF] Tao | Nonlinear global instability in buoyancy-driven boundary-layer flows[END_REF] has been able to confirm quantitatively and relate the results from a convective-absolute linear stability theory to the existence of a global instability with sustained oscillations in the numerical simulation of a two dimensional flow developing along a vertical plate in a stratified environment, none of these studies have succeeded in thoroughly explaining the relationship between the instability studies of the buoyancy layer and the transition to unsteadiness in the isoflux cavity.

As was said above, one difficulty when dealing with time integration of the unsteady equations lies in the fact that the stratification changes with Rayleigh number. There are two ways to get around the stratification issue. One can think of building an initial condition as close as possible from the asymptotic solution at a given value of Ra. This is quite simple owing to particular structure of the flow, 1D except in the ends, and linear stratification everywhere. To this end one needs the stratification whose dimensional value was reported by Kimura and Bejan [START_REF] Kimura | The boundary layer natural convection regime in a rectangular cavity with uniform heat flux from the side[END_REF] in the conclusion of their paper as follows :

( dT dz
) core = (8192) For the sake of completeness computation of the stratification is reported in annex A. The magnitude of this vertical stratification can be appreciated by comparing the resulting vertical temperature difference ∆T H to the temperature difference across the cavity ∆T L = 8 q φH λ Ra -1/5 . This ratio reads :

∆T H ∆T L = (2) 2/9 4 ( H L ) 5/9 Ra 1/9 (13)
In a square cavity at Ra = 10 9 , ∆T H is approximately 3 times larger that ∆T L , a stratification much larger than that found in the isothermal differentially heated cavity. This ratio increases with Ra and for Ra = 10 12 , it reaches 6.3, a very large stratification that obviously helps strongly stabilize the flow. Another way to circumvent the stratification issue is to use a scaling which removes the Ra dependence of the core stratification, that is, as equation ( 12) shows, to define a temperature scale such as ∆T = κν gβH 3 Ra 8/9 [START_REF] Prandtl | Essentials of Fluid Dynamics[END_REF] or, equivalently :

∆T = φH λ Ra -1/9 (15)
which is precisely the temperature scale proposed in [START_REF] Kouroudis | Direct numerical simulation of natural convection in a square cavity with uniform heat fluxes at the vertical sides: Flow structure and transition[END_REF].

Governing equations

Scalings

We now review and discuss the different scalings to address these questions. We assume that the flow is governed by the Navier-Stokes equations under the Boussinesq approximation.

∂u ∂x + ∂w ∂z = 0 (16)
∂u ∂t + u ∂u ∂x + w ∂u ∂z = - ∂P ∂x + ν ∇ 2 u ( 17
)
∂w ∂t + u ∂w ∂x + w ∂w ∂z = - ∂P ∂z + ν ∇ 2 w + gβ(T -T 0 ) (18)
∂θ ∂t + u ∂θ ∂x + w ∂θ ∂z = κ ∇ 2 θ ( 19
)
The equations are made dimensionless using a length scale L, a reference velocity Ṽ , and a temperature scale ∆T , yielding :

∂u ∂x + ∂w ∂z = 0 (20
( H L ) 4/9 ∆T = φH λ Ra -1/9
Table 1: Scaling I is the natural scaling along an isoflux boundary layer; scaling II results from Kimura and Bejan's analysis; scaling III corresponds to a constant stratification. This is the scaling used in [START_REF] Kouroudis | Direct numerical simulation of natural convection in a square cavity with uniform heat fluxes at the vertical sides: Flow structure and transition[END_REF] 5

In order to illustrate the effect of using these different scalings on the numerical computation of the solution at a given Ra number, the following test cases are performed with a numerical algorithm which we briefly describe.

Numerical algorithm

The dimensionless equations are integrated in unsteady form. The time stepping algorithm is based on the incremental prediction-projection method ( [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two-or three-dimensional cavity flows[END_REF], [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]). The prediction step consists in obtaining an intermediate velocity field V = (u , w ), assuming a known pressure field P n . The time discretization combines an implicit treatment of the linear viscous terms whereas the convective terms are treated by a second order extrapolation. This scheme reads:

3u -4u n + u n-1 2∆t + 2V n • ∇u n -V n-1 • ∇u n-1 = - ∂P n ∂x + P r P e ∇ 2 u (24) 3w -4w n + w n-1 2∆t + 2V n • ∇w n -V n-1 • ∇w n-1 = - ∂P n ∂z + P r P e ∇ 2 w + B θ n+1 (25)
which yields two independent Helmholtz problems for each of the components u and w . These Helmholtz problems are solved using an incremental ADI solver.

The second step consists in projecting V onto the space divergence free vector fields with null normal trace, yielding V n+1 . This is carried out by defining a scalar variable ψ such as:

V n+1 -V = ∆t∇ψ (26)
ψ is determined by taking the divergence of equation ( 26) which yields,

∇ 2 ψ = -(∇ • V )/∆t (27)
This elliptical equation for ψ is associated with boundary condition ∂ψ/∂n = 0. This equation is solved using a multigrid algorithm. The final velocity is obtained from:

V n+1 = V + ∆t∇ψ (28)
The variable ψ can be seen as a pressure correction and one has:

P n+1 = P n + 3/2 ψ (29)
The spatial discretization makes use of the classical staggered grid arrangement. The discretized equations are obtained by a finite volume formulation, which guarantees conservative properties.

Both the viscous and convective terms are treated by finite centered differences. We use a grid of N × M grid points with a cosine mesh in x and a uniform mesh in z.

Reaching steady state through time integration

4.1. Which scaling to use ?

Figure 1 presents the time traces of the temperature at a point located in the core for Ra = 10 9 starting from rest and uniform temperature. The three time traces, although computed with different time steps, are virtually indistinguishable which confirm the correctness of the scalings proposed in table 1 (if needed). They also show that the time needed to reach steady state is on the order of 1000 time units in the time unit H 2 κ Ra -2/5 . The high frequency oscillations observed up to a time of 200 are due to the internal waves. This figure confirms that in this case it is the diffusive time scale on the cavity height which governs the steady state, and not the time needed to damp the internal waves, contrarily to the isothermal cavity. As can also be seen in figure 2, the period of the internal waves decreases with increasing time, as a consequence of the fact that the core stratification increases with time. In order to further illustrate the effect of using these different scalings on the numerical computation of the solution at a given Ra number, the following test case is performed. The solution for Ra = 10 10 is computed for different scalings or initial conditions. Figure 3 presents the time evolution of the temperature at a given point for the following conditions. Time trace (a) corresponds to the integration of the governing equations using scaling I starting from an initial solution for Ra = 10 9 as initial condition. For time trace (b) the initial condition is the manufactured solution 6 corresponding to Ra = 10 10 using scaling I form of the governing equations. Time trace (c) corresponds to the use of scaling III starting from an initial solution for Ra = 10 9 as initial condition. The three time traces show that the asymptotic solution is steady, but one clearly sees that one can conclude much more rapidly on the asymptotic nature of the solution either starting from the manufactured initial condition or using scaling III. All the computations reported below have been performed using scaling I.

Internal gravity waves

As recalled in section 2, it was shown in [START_REF] Patterson | On the existence of an oscillatory approach to steady natural convection in cavities[END_REF] that the stratified core is able to sustain internal gravity waves. A characteristic angular frequency of these waves is the Brunt-Väisälä frequency µ based on the stratification in the core (µ = √ BS with the notations of Table 1). In linearly stratified enclosures Thorpe [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] showed that the spectrum of the standing waves is discrete and quantified according to

τ 2 n,m = µ 2 m 2 /L 2 m 2 /L 2 + n 2 /H 2 = µ 2 1 1 + n 2 m 2 1 AR 2 (30)
where µ is the basic angular Brunt-Väisälä frequency, and n and m are the number of half wavelengths in the vertical and horizontal directions respectively, showing that the basic Brunt-Väisälä frequency is an upper bound of the waves frequency. In shallow cavities (AR << 1), the spectrum thus extends from µ down to ≈ µ × AR as can be seen from figure 4 that presents the spectrum of the Nusselt number across the vertical mid-plane, known to be a good indicator of the internal wave activity, in a cavity of AR = 0.2 for Ra = 2.2 × 10 12 . The spectrum clearly shows different peaks corresponding to odd values of m, that match the analytical formula to better than 1 0 / 00 . That n should be taken equal to 1 is evidenced by figure 5 which shows that at a given abscissa the crests and troughs of the fluctuations are all in phase.

Transition to unsteadiness in a square cavity

Transition through time integration

Transition to unsteadiness was first determined in the brute force way, by setting the value of Ra to a specified value and observing the asymptotic response. To avoid the stratification issue we used the first strategy presented above, i.e., starting the time integration from the initial manufactured condition previously described at the given Ra number. For instance figure 6a presents the time evolution of the temperature at several points located in the boundary layer in a square cavity for Ra = 3 × 10 12 . It can be seen that, right after the beginning of time integration, a wave packet develops and travels upward as time increases but the solution eventually reaches a steady state. The situation is very similar to what was reported in [START_REF] Armfield | A direct boundary-layer stability analysis of steady-state cavity convection flow[END_REF] for the isothermal cavity. A close up of the evolution is displayed in figure 6b which shows that the front and rear of the wave packet progress linearly with time but not at the same speed, resulting in a spreading of the wave packet. This indicates that the boundary layer has become convectively unstable but the asymptotic solution is still steady.

When the Ra number is increased to Ra = 4.2 × 10 12 , the time traces of the temperatures in figure 7b and7a show that after the initial wave packet has left the boundary layer, the solution at the two uppermost points continues to display time oscillations, first characterized by several frequencies, which slowly give way to a time periodic solution, as confirmed by the power spectrum of the signal for the uppermost point displayed in figure 8.

Figure 9 presents a snapshot of the evolution of the temperature fluctuation ( θ = θ(x, z, t)θ(x, z)) in the upward boundary layer as a function of altitude. It shows that the amplitude of the fluctuation grows exponentially with z in the upper middle part of the cavity.

To confirm this we have computed the most unstable eigenmode of the steady state solution at this value of the Rayleigh number. The (unstable) steady solution was obtained using the procedure proposed by Tuckerman et al in a series of papers [START_REF] Tuckerman | Steady-state solving via Stokes preconditioning; recursion relations for elliptic operators[END_REF][START_REF] Mamun | Asymmetry and Hopf bifurcation in spherical Couette flow[END_REF][START_REF] Tuckerman | Bifurcation analysis for timesteppers[END_REF] that we implemented successfully in [START_REF] Gadoin | A general methodology for investigating flow instabilities in complex geometries: application to natural convection in enclosures[END_REF]. The unsteady linearized equations around that steady solution are then integrated in time starting from an initial random temperature field. The equations are inte- grated long enough in time to let the solution polarize along the most unstable eigenmode. As this may require very long time integration, the solution was periodically rescaled to prevent overflow as shown in figure 10. Figure 11 displays the evolution of the amplitude of the eigenmode along the same vertical line as in figure 9. It shows that the eigenmode possesses a spatial structure that evolves exponentially with z with an amplification of almost 10 orders of magnitude from the bottom to the top. This figure also helps realizing the difficulty to numerically approximate this problem accurately as the eigenmode shows there is room for at least 20 wavelengths all along the boundary layer. Assuming that one needs approximately 30 points per wavelength to accurately represent the time evolution of the solution results in the fact that one should have on the order of 500 to 1000 grid points in the vertical direction. It also shows that a uniform equidistant mesh in the vertical direction is probably the best way to approximate this solution.

Comparison of figures 9 and 11 shows that the fluctuating temperature differs from the eigenmode in the half lower part of the boundary layer where the fluctuation are on the order of 10 -4 whereas the amplitude of the eigenmode continues to decrease exponentially. This is due to small albeit unavoidable nonlinear effects that mask the shape of the eigenmode.

Grid independence was performed increasing the spatial resolution to (N, M ) = (512, 1024). The unsteady solution obtained for Ra = 4.2 × 10 12 with (N, M ) = (512, 512) was interpolated to the finer mesh and time integration was resumed. The time signal showed that the solution converged to a steady state. We had to increase the Rayleigh number to Ra = 4.5 × 10 12

Unsteadiness via spectrum computation 5.2.1. Methodology

One way to confirm the previous results is to compute the evolution of the spectrum of the Jacobian of the steady solution of the equations. As recalled in the introduction, computation of the spectrum requires being able to compute first the corresponding steady solution. Since the time needed to reach the steady state solution diverges as the Rayleigh number approaches the critical value computation of steady solutions require the use of iterative algorithms. Since most fixed point methods more or less mimic unsteady algorithms and therefore encounter the same convergence difficulties, there is no really other choice than to resort to Newton's iteration.

Given a non-linear equation

F (u) = 0
Newton's iteration consists in repeating iteratively until convergence

DF (u k ) δu k = -F (u k ) ( 31
)
u k+1 = u k + δu k
where DF (u k ) is the Jacobian of F on the current iterate u k . It requires solving equation [START_REF] Gjesdal | Spectral element benchmark simulations of natural convection in two-dimensional cavities[END_REF] which rises several difficulties. The Jacobian belongs to the class of large sparse matrices. Its order is on the order of 4N M as it is made of writing the discrete linearized equations linking the discrete values of (u, w, P, T ). Ranking each variable in lexicographic order, figure 13a displays the typical Jacobian pattern showing its sparsity pattern and typical bandwidth. These features make direct methods intractable and has motivated the development of alternative approaches such as those proposed by Tuckerman and co-workers ( [START_REF] Tuckerman | Steady-state solving via Stokes preconditioning; recursion relations for elliptic operators[END_REF], [START_REF] Edwards | Krylov methods for the incompressible Navier-Stokes equations[END_REF], [START_REF] Mamun | Asymmetry and Hopf bifurcation in spherical Couette flow[END_REF], [START_REF] Tuckerman | Bifurcation analysis for timesteppers[END_REF]). The basis of the method consists in using a time stepping algorithm which relies on a classical implicit discretization of the linear terms coupled to an explicit discretization of the non-linear terms as a means of getting iteratively the action of the Jacobian on a given vector, requiring the use of a matrix free method to solve [START_REF] Gjesdal | Spectral element benchmark simulations of natural convection in two-dimensional cavities[END_REF]. Although effective at moderate spatial resolutions, this solution procedure was found to suffer convergence difficulties for spatial resolutions required for the present configuration.

We have thus decided to resort to direct resolution of (31) via LU-factorisation. This procedure will also allow us to use the Arpack package to compute the leading eigenvalues in the shift and invert mode which requires being able to compute efficiently (DF -σI) -1 where σ is a complex shift.

As said above, dealing with the Jacobian in the form depicted in figure 13a is intractable. Two manipulations are needed, on the one hand to handle the sparsity and on the other to reduce the bandwidth. We have developed a hand-made renumbering algorithm which basically, without going in too much detail, consists in grouping together the unknowns pertaining to a given cell which reduces the bandwidth to approximately 4N (see figure 13b). Storage and manipulation is done under the CSR (Compressed Sparse Row) format, using routines from the Sparskit tool box ( [START_REF] Saad | SPARSKIT: a basic tool kit for sparse matrix computations -version 2[END_REF]). It is also well known that the Jacobian has to be regularized since its kernel is at least of dimension 1. This fundamental mode in incompressible flow is related to the fact that a constant pressure results in a null velocity field. The regularization thus consists in imposing a value of pressure at a given point replacing the corresponding divergence equation. In the case of the isoflux cavity, the dimension of the Jacobian's kernel is at least 2 since with Neumann boundary conditions at all 4 walls, the temperature field is determined up to a constant. It is thus necessary to drop one temperature equation and replace it by setting a given value at one point or any linear combination to make the Jacobian non singular. In order to maintain the centro-symmetry of the solution a linear combination linking the 4 points around the cavity center was used, which also requires starting the iterations with an initial guess that shares this property. The resulting regularized Jacobian is then solved (31) using the PARDISO package ( [START_REF] Alappat | A recursive algebraic coloring technique for hardware-efficient symmetric sparse matrix-vector multiplication[END_REF][START_REF] Bollhöfer | State-of-the-Art Sparse Direct Solvers[END_REF][START_REF] Bollhöfer | Large-scale sparse inverse covariance matrix estimation[END_REF]).

Putting this altogether enables us to compute steady solution of the governing equations starting increasing progressively the Rayleigh number value step by step. Usually 6 to 8 iterations are required to bring the residual of the steady equations down to machine accuracy. Increments of 10 are taken up Ra of 10 9 , increments of 2 to 3 for larger values.

In order to reach Ra values of 10 12 and above, we had to increase the spatial resolution up to (N, M ) = (512, 512). This corresponds to a Jacobian of order ≈ 10 6 . Our reordering algorithm results in a bandwidth of ≈ 2500. Figure 14 presents the corresponding temperature field and stream function for Ra = 10 12 . Its vertical stratification and temperature difference across the cavity are 9.88 and 1.591 which, not surprisingly, are in very close agreement with the analytical ones (9.99 and 1.593).

Spectrum computation

The spectrum of a given steady solution was computed using the Arpack library ( [START_REF] Lehoucq | ARPACK users guide: Solution of large scale eigenvalue problems by implicitly restarted Arnoldi methods[END_REF], [START_REF] Sorensen | Numerical methods for large eigenvalue problems[END_REF]). Since one anticipates that the loss of stability be due to a Hopf bifurcation, the eigenvalues of interest are those of maximum real part with an imaginary part corresponding to the frequency of the observed unsteadiness, making the algorithm genuinely complex. We have used the generalized eigenvalue problem

DF (u) x = λ M x ( 32
)
where M is the identity matrix with zeros on the diagonal elements corresponding to the enforcement of the incompressibility condition. The procedure is very similar to what has been done by others ([57], [START_REF] Yu | Stability of convective flows in cavities: solution of benchmark problems by a low-order finite volume method[END_REF]). We have found by trial and error that the shift-invert algorithm generally converges when requiring a few tenths eigenvalues. Getting a global view of the spectrum thus requires repeating the procedure with different shifts σ judiciously scattered in the complex plane. Figure 15 presents a global view of the spectrum for Ra = 10 12 for a spatial resolution (N, M ) = (512, 512). One can distinguish several groups of eigenvalues and corresponding eigenmodes. The eigenvalues of imaginary part larger than 7 correspond to boundary layer modes (fig. 16a). At this Ra value they are still very stable and we will come back to them later. The second group of modes is characterized by eigenvalues of imaginary part in the range 3 to 5. They correspond to short wavelength traveling waves along the ceiling and floor of the cavity connected with long wavelength waves along the vertical walls (fig. 16b). The third group of modes correspond to eigenvalues smaller that the cut-off Brunt-Väisälä frequency which is equal to 2.66. These are internal gravity modes, whose typical structure is displayed in fig. 16c. As can be seen they are characterized by inclined stripes. The eigenmode which is shown corresponds to an imaginary part equal to 2.22 and the stripes angle to the horizontal is very close to 57 • , which agrees with the theoretical inclination arcsin( 2. 22 2.66 ) to better than 1 • . The fourth group of modes are those of null imaginary part. They are characterized by horizontal stripes as shown in fig. 16d. From now on, we will concentrate on the eigenvalues of imaginary part larger than 8, taking for granted, unless proven wrong, that they are responsible for the transition to unsteadiness.

As can be seen the boundary layer modes lie on two parabola-like curves open to the left. They come in pairs with nearly identical imaginary parts. Those located on the rightmost parabola have the anti centro-symmetry property (ACS) of the base flow (ie θ(x, z) = -θ(-x, -z) assuming the origin is at the cavity center), whereas those on the leftmost parabola have the opposite symmetry (CS) (ie θ(x, z) = θ(-x, -z)). From this standpoint the situation is different from that in the isothermal cavity where both families of modes ACS and CS were found to lie on the same parabola [START_REF] Gadoin | A general methodology for investigating flow instabilities in complex geometries: application to natural convection in enclosures[END_REF]. As Ra increases, figure 17 shows that the eigenvalues of imaginary part around 8 move to the right until one of them will eventually cross the imaginary axis which is about to happen for Ra = 4 × 10 12 . As can be seen the imaginary part of most dangerous eigenvalues also increases with Ra, which is due to the fact that the boundary layer gets thinner with increasing Ra and therefore the wavelength decreases, resulting in an increased frequency.

Figure 18 shows part of the spectrum for Ra = 4 × 10 12 and 4.2 × 10 12 showing that for the latter value, one eigenvalue has just crossed the imaginary axis. From this spectrum it could be concluded that for a spatial resolution (N, M ) = (512, 512) the critical Rayleigh number corresponding to the transition to unsteadiness lies in between Ra = 4 × 10 12 and The spectrum for Ra = 4.2 × 10 12 also helps understand some of the characteristic features of time signals shown previously. The two most unstable eigenvalues are (1.52 × 10 -3 , 9.65) and (-2.85 × 10 -3 , 9.87). The difference between the two angular frequencies (= 0.22) correspond approximately to a time length of ≈ 28 which is precisely the period of the low frequency modulation that is seen at the beginning of the time signal shown in figure 7b. The low frequency modulation observed during the transients before reaching a time-periodic evolution is thus a beat frequency between the most unstable modes. On the other hand the small difference between the real parts ≈ 4.3 × 10 -3 explains the long time needed to separate the most unstable eigenmodes. Getting completely rid of the second most unstable mode while integrating the linearized equations would require on the order of log (10 15 ) 4.3×10 -3 ≈ 8000 time units. Grid refinement tests are shown in figure 19 that displays the evolution of the spectrum for Ra = 4.4 × 10 12 . The figure clearly shows that increasing spatial resolution has a stabilizing effect on the solution. Whereas the solution for (N, M ) = (512, 512) is unsteady with 3 eigenmodes unstable, that for (N, M ) = (512, 768) is just marginally unstable while that for (N, M ) = (512, 1024) is definitively stable. It is also seen that increasing spatial resolution has a substantial effect on the frequency of oscillation. One can conclude that a spatial resolution of (N, M ) = (512, 512) is not sufficient enough to accurately approximate this solution. Finally figure 20 shows that for (N, M ) = (512, 1024) transition to unsteadiness takes place slightly below Ra = 4.5 × 10 12 which is consistent with the results from time integration presented earlier. The imaginary part of the most unstable is equal to 9.735 which corresponds to a frequency of 1.549 which agrees to better than 1 0 / 00 with the dominant frequency of the time signal (1.550) in figure 12.

Figure 21 displays the amplitude of the θ-component of the most unstable eigennmode along a vertical line parallel to the heated wall showing the exponential increase of the mode amplitude as a function of altitude. This figure is in very good agreement with the most unstable eigenmode obtained through integration of the linearized equations shown in figure 11, validating both approaches.

Global results

Table 2 presents the results of many hours of computation corresponding to determinations of lower and upper bounds for the transition to unsteadiness using both time integration and Jacobian's spectrum computation, which are very complementary. As said above, computation of part of the spectrum is performed using ARPACK in the shift-and-invert mode requiring a few tens of eigenvalues. The outcome of the procedure thus yields eigenmodes corresponding to eigenvalues in the vicinity of the shift. Good guessed values for the shift, in particular its imaginary part, is provided by time integration. It also provides a good bracket of the lower and upper bounds, thus helping reduce both the number of solutions computed with Newton's method as well as the number of calls to ARPACK with different shifts.

The steady-unsteady transition curve given in figure 22, the value which is plotted is determined by linear interpolation. This figure shows that the critical values obtained by blending the linear stability analysis with the isoflux cavity characteristics are not relevant. The transition to unsteadiness occurs at a value of Ra more than two orders magnitude larger that the predicted value. Likewise the aspect ratio dependence is not followed. This shows that the isoflux cavity has its own stability characteristics which are not those of the vertical boundary layers it houses despite the fact that the boundary layers being of constant thickness and that, due to the large number of wavelengths, one could imagine that the enclosure effect would be negligible. However comparison of the rightmost columns shows that at criticality the characteristic values that characterize the traveling waves agree reasonably well with those inferred from from blending the linear stability results from [START_REF] Desrayaud | Stability of the flow near a heat-flux plate and comparison with numerical simulation in a square cavity[END_REF] and the analytical values from [START_REF] Kimura | The boundary layer natural convection regime in a rectangular cavity with uniform heat flux from the side[END_REF].

In particular for large values of the aspect ratio the critical wave-speeds agree to within a few percents.

Figure 23 shows comparisons of the amplitude of the eigenmode and of fluctuating temperature in a cavity of AR = 5 for Ra = 5.1 × 10 13 . These evolutions display similar trends as those already observed for the square cavity but with an increased number of wavelengths. The corresponding wavelength is equal to ≈ 0.188, showing again the need to have at least 1000 grid points in z, and that a uniform mesh is indeed the best way to approximate this solution.

Les us conclude this section by reporting some surprising behaviors of the eigenvalues. As said above for the square cavity the eigenvalues lie on two parabola open to the left, those corresponding to modes having the anti-centro-symmetry (ACS) of the base flow on the right most parabola and those having the opposite symmetry (CS) on the left one. The modes on both parabolas have almost exactly identical imaginary parts (see figures 17 and 18 for

Rayleigh

Aspect ratio [START_REF] Desrayaud | Stability of the flow near a heat-flux plate and comparison with numerical simulation in a square cavity[END_REF] and the analytical values from [START_REF] Kimura | The boundary layer natural convection regime in a rectangular cavity with uniform heat flux from the side[END_REF] instance), corresponding to the fact that they have the same spatial structure, same number of wavelengths, with opposite symmetry. This is no longer true for other values of the aspect ratio. For instance in tall cavities figure 24 shows that the eigenvalues lie on a single parabola and we have checked that the modes with increasing imaginary part have alternatively ACS and CS symmetry. The situation is analogous to what was found in isothermal cavities ( [START_REF] Gadoin | A general methodology for investigating flow instabilities in complex geometries: application to natural convection in enclosures[END_REF], [START_REF] Xin | Natural convection flows in air-filled differentially heated cavities with adiabatic horizontal walls[END_REF]) presumably due to a coupling of both boundary layers eigenmodes through the small horizontal extension of the boundary layer eigenmodes which is apparent in figure 16a. This characteristic (all eigenvalues on the same parabola) extends to values of the aspect ratio down to 3.

Intermediate plots of the eigenvalues location rearrangement can be seen in figure 25 for AR = 2 which shows that the single parabola found for large values of the aspect ratio has started to split into two parabolas, the eigenmodes with ACS lying on the rightmost one and those with CS on the leftmost. The imaginary parts of the CS modes remain approximately half way between the two neighboring ACS modes. With decreasing aspect ratio, the two parabolas move further apart (see figure 26 for AR = 1.2). The rearrangement of the imaginary parts takes place on a very narrow range of aspect ratios (1-1.1), as shown in figure 27, just after the two parabolas have started to move back towards one another. As the aspect ratio continues to decrease, the two parabolas get very close, giving rise to pairs of nearby eigenvalues, as can be seen from figure 28 for AR = 0.5. For smaller values of the aspect ratio this process continues, giving rise to eigenvalues of multiplicity two (to better than 10 -5 ), corresponding to the fact that there is no reason why two eigenmodes of identical spatial structure only distinct in their ACS or CS symmetry should have different stability properties. What remains unexplained is what is seen for nearly square cavities. Why modes with different symmetries have so different stability criteria ?

Conclusion

We have computed the transition to unsteadiness in air filled differentially heated isoflux cavities for values of the aspect ratio in the range 0.2-8. Transition to unsteadiness was investigated using three different approaches : through time integration which required revisiting scaling issues ; through time integration of the linearized equations, and finally from direct computation of the most dangerous part of the Jacobian of steady solutions, either stable or unstable. It was found that for a square cavity filled with air (P r = 0.71), transition occurs for a Rayleigh number very close to 4.5 × 10 12 a very large value which requires that the solution be computed with at least 1000 grid points in the vertical direction. We have discussed the nature of eigenmodes, showing that they belong to four different groups corresponding to different physical phenomenology. We have shown that the modes responsible for the transition are characterized by a wavelike structure of constant wavelength whose amplitude increases exponentially with height in the upward boundary layer along the heated wall (and symmetrically along the cooled wall). We have shown that these cavities possess their own stability characteristics, although, at criticality, the traveling waves in the boundary layers have similar characteristics than those resulting from blending the linear stability analysis of the buoyancy layer with the analytical solution. We have finally discussed the evolution of the spectrum pattern observed for different values of the aspect ratio and the way the leading eigenvalues corresponding to eigenmodes with alternate symmetry which lie on a single parabola for large values of the aspect ratio eventually merge in eigenvalues of multiplicity 2 at small aspect ratio.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: Evolution of the temperature at point (0.16, 0.62) for the 3 different scalings, Ra = 10 9 , AR = 1. Numerical values for scalings II and III have been converted to scaling I units

FrequencyFigure 4 :

 4 Figure 4: Density power spectrum of the evolution of the Nusselt number across the vertical mid-plane in the time evolution leading to steady state. Ra = 2.2 × 10 12 ; AR = 0.2. The basic Brunt-Väisälä frequency µ is equal to 0.307

Figure 5 :

 5 Figure 5: Isovalues of the fluctuating temperature field during the transient regime for Ra = 2.5×10 12 , AR = 0.2 ; isovalues ; ±5e -5, ±4e -5, ±3e -5, (N, M ) = (512, 256)

  Evolution of the temperature very close to the heated wall x ≈ 0.0001 and equidistant values in z = 0.25, 0.375, 0.5, 0.625, 0.75, 0.875 and 0.98.

  Close-up over the first 80 time units

Figure 6 :

 6 Figure 6: Time traces of temperature: Ra = 3 × 10 12 ; (N, M ) = (512, 512)

Figure 7 :FrequencyFigure 8 :Figure 9 : 2 TimeFigure 10 :

 789210 Figure 7: Time traces of temperature very close to the heated wall x ≈ 0.0001: Ra = 4.2 × 10 12 ; (N, M ) = (512, 512)

Figure 11 :

 11 Figure 11: Evolution of log(| θ|) with z ; Ra = 4.2 × 10 12 ; (N, M ) = (512, 512) ; full line corresponds to positive values, dashed line to negative ones ; one wavelength is thus made of two consecutive arches

- 3 Figure 12 :

 312 Figure 12: Time trace of temperature: Ra = 4.5 × 10 12 ; (N, M ) = (512, 1024)

  (a) (u,w,P,T) and corresponding equations ranked in natural lexicographic order ; the figures on the left correspond to the numbering of the equations; the 500 first lines are the equations for u, the lines from 500 to 1000 are those for w, the next 500 are those for the divergence and the last 500 are those for Θ (b) after reordering

Figure 13 :

 13 Figure 13: Pattern of Jacobian ; resolution N = 16, M = 32; Jacobian is of order ≈ 2000

Figure 14 :

 14 Figure 14: Solution for Ra = 10 12

Figure 15 :

 15 Figure 15: Global view of the most unstable part of the spectrum for Ra = 1 × 10 12 ; spatial resolution (N, M ) = (512, 512) ; (+) group I modes ; (×) group II; ( * ) group III; ( ) group IV. Only half of complex plane with positive imaginary part is shown.

Figure 16 :

 16 Figure 16: Real part of temperature component of eigenmodes. The interface line between yellow and dark blue is the zero amplitude contour. Amplitude is arbitrary

Figure 17 :

 17 Figure 17: Evolution of part of spectrum with Ra; Ra = 1 × 10 12 (+) ; Ra = 2 × 10 12 (×); Ra = 3 × 10 12 ( * ); Ra = 4 × 10 12 ( ); spatial resolution (N, M ) = (512, 512)

Figure 18 :Figure 19 :

 1819 Figure 18: Evolution of part of spectrum for Ra = 4 × 10 12 (×) and 4.2 × 10 12 (+); spatial resolution (N, M ) = (512, 512)

Figure 20 :

 20 Figure 20: Evolution of part of spectrum for Ra = 4.4 × 10 12 (+) , 4.5 × 10 12 (×) , 4.6 × 10 12 ( * ); spatial resolution : (N, M ) = (512, 1024)

Figure 21 :

 21 Figure 21: Evolution of log 10(|Θ|) in z ; spatial resolution (N, M ) = (512, 1024)

Figure 22 :

 22 Figure 22: Evolution of the critical Rayleigh number as function of aspect ratio

Figure 23 :

 23 Figure 23: Evolution of log 10(| Θ|) with z ; Ra = 5.1 × 10 13 ; AR = 5 ; (N, M ) = (512, 1024)

  012 -0.01 -0.008 -0.006 -0.004 -0

Figure 24 :Figure 25 :

 2425 Figure 24: Evolution of part of spectrum for AR = 8 ; Ra = 2 × 10 14 (+) ; Ra = 2.05 × 10 14 (×); Ra = 2.1 × 10 14 ( * ) ; (N, M ) = (256, 1536)

Figure 26 :Figure 27 :

 2627 Figure 26: Evolution of part of spectrum for AR = 1.2 ; Ra = 5 × 10 12 (+) ; Ra = 6 × 10 12 (×); (N, M ) = (512, 1024)

Figure 28 :

 28 Figure 28: Evolution of part of spectrum for AR = 0.5 ; Ra = 2.1 × 10 12 (+) ; Ra = 2.4 × 10 12 (×); Ra = 2.7 × 10 12 ( * ); (N, M ) = (1024, 512)

Table 1

 1 presents the values of parameters P e and B, of the temperature boundary condition and resulting stratification for different choices of scalings (it is assumed that L = H).

	)

Table 2 :

 2 Global values for 0.2 ≤ AR ≤ 8 ; σ i , λ and c are the angular frequency of the first unstable mode, its wavelength (measured graphically), and corresponding wave-speed, respectively ; σ th i , λ th and c th are the theoretical values obtained from blending the linear stability results from

	H L	N × M	Ra inf	Ra sup	σ i	λ	c = λσ i 2π	σ th i	λ th	c th
	0.2 1024 × 512 1.7 × 10 12 1.9 × 10 12 10.50 0.0630 0.105 15.37 0.054	0.132
	0.25 1024 × 512 1.9 × 10 12 2 × 10 12 10.12 0.0628 0.101 14.75 0.052	0.122
	0.5 1024 × 512 2.4 × 10 12 2.7 × 10 12 9.57	0.052	0.079 13.20 0.0453 0.095
	1	512 × 1024 4.4 × 10 12 4.5 × 10 12 9.73 0.0435 0.067 12.39 0.0368 0.0726
	2	512 × 1024 8.8 × 10 12 9 × 10 12 10.15 0.032	0.052 11.81 0.0293 0.0551
	3	512 × 1024 1.6 × 10 13 1.8 × 10 13 10.66 0.0261 0.0442 11.94 0.0242 0.0461
	4	512 × 1024 3 × 10 13 3.1 × 10 13 11.17 0.0217 0.0386 12.25 0.0206 0.0403
	5	512 × 1024 5 × 10 13 5.2 × 10 13 11.66 0.0188 0.0349 12.64 0.0179 0.0361
	6	512 × 1024 8.1 × 10 13 8.2 × 10 13 12.12 0.0162 0.0313 13.05 0.0159 0.0329
	8	256 × 1536 2.0 × 10 14 2.1 × 10 14 13.36 0.013 0.0276 14.14 0.0125 0.0281

Apparently the work by Kimura and Bejan was unknown to these authors, as it is not referenced, although it appears in the references of some papers they quote

In his study[START_REF] Yahata | Stability analysis of natural convection in vertical cavities with lateral heating[END_REF] Yahata wrongly labeled some parts of the instability boundary as due to internal gravity modes. In a stably stratified environment the internal gravity modes are always stable although they can be part of the oscillating motion if the primary instability mechanism has a frequency lower that the cut-off Brunt-Väisälä frequency as shown in[START_REF] Xin | Natural convection flows in air-filled differentially heated cavities with adiabatic horizontal walls[END_REF]

Since a good agreement between the analytical value and numerical simulations was reported in[START_REF] Kimura | The boundary layer natural convection regime in a rectangular cavity with uniform heat flux from the side[END_REF] (see their figure6), comparison of the two expressions[START_REF] Kimura | The boundary layer natural convection regime in a rectangular cavity with uniform heat flux from the side[END_REF] and[START_REF] Kouroudis | Direct numerical simulation of natural convection in a square cavity with uniform heat fluxes at the vertical sides: Flow structure and transition[END_REF] strongly suggests there was a misprint in the expression reported in the conclusion of[START_REF] Kimura | The boundary layer natural convection regime in a rectangular cavity with uniform heat flux from the side[END_REF]. The exponent should have been[START_REF] Batchelor | Heat transfer by free convection across a closed cavity between vertical boundaries at different temperatures[END_REF] 9 and not 1 9 , as confirmed by A. Bejan and S.Kimura (personal communication)

Note however that their dimensionless equations[START_REF] Quon | Free convection in an enclosure revisited[END_REF][START_REF] Bejan | Note on Gill's solution for free convection in a vertical enclosure[END_REF][START_REF] Kimura | The boundary layer natural convection regime in a rectangular cavity with uniform heat flux from the side[END_REF] are not consistent with these scalings

i.e. initializing a centro-symmetric solution made of two buoyancy layers[START_REF] Elder | Laminar free convection in a vertical slot[END_REF] along each wall with corresponding parameters[START_REF] Gill | The boundary-layer regime for convection in a rectangular cavity[END_REF][START_REF] Kouroudis | Direct numerical simulation of natural convection in a square cavity with uniform heat fluxes at the vertical sides: Flow structure and transition[END_REF]

What remains unexplained is the reason for which the critical value of the Rayleigh number for the buoyancy layer and the isoflux cavity differ so widely.

Appendix A

Computation of the stratification comes from balancing the vertical conduction of heat with the heat transported by convection in the boundary layers, that is :

where L is the cavity width. In the boundary layer regime, using the results from [ [START_REF] Kimura | The boundary layer natural convection regime in a rectangular cavity with uniform heat flux from the side[END_REF]], this expression reduces to :

x)dx [START_REF] Ravi | On the high-Rayleigh-number structure of steady laminar natural-convection flow in a square enclosure[END_REF] The integral in the expression above is exactly 1 2q , which results in :

and upon replacing q 5 from ( 8), the dimensional stratification finally reads :