Convection in solid planetary interiors and implications
for their evolution
S. Labrosse
LGL-TPE, Université Claude Bernard Lyon 1, Ens de Lyon, CNRSFDEPS, Nov. 2019

Subsolidus convection

I Solid surface planets and planetary objects (icy satellites, dwarf planets) show signs ofdeformation in the solid state, whether active or in their past.
I In many cases: thermal convection.
I Very large viscosity =⇒ slow motion. The bottleneck for the thermal evolution of planetary
objects with solid surface.
I In many cases, a liquid layer exists below (metallic core, water ocean) whose dynamics iscontrolled by the rate of heat extraction by convection in the solid.
I Convection can also happen in solid shells or spheres deep inside planetary objects: inner core, HP
ice layers of Titan, Ganymede.
=⇒ Important to understand convection in planetary mantles.

Part I
Fundamentals of Rayleigh-Bénard convection

Outline

Fundamentals of Rayleigh-Bénard convection
Historical backgroundbalance equations and the Boussinesq approximation
Linear stability analysis
Behaviour beyond the onset
High Rayleigh number dynamics and scaling of heat transfer

Experiments by Bénard
Bénard (1900a,b, 1901) conducted the first systematic experiments on flow driven by a destabilisingtemperature difference.

CL"'I44I'U;U~
la p~J~et/o~po;<A<m/~ue t/~oset'ig. 6. – .4/?pare~de courants de convection dans une nappe mince de 7j'ou~cnon i'olalil à 100" el pour assurer 7'ui!)/b7'7H/(c/?&r/a/!c descu~'econditions dans toute J'étendue de la(échelle 1/4'.

)'ig.i3.–)tt~À'<7'YuJe/t<7't'crcp~sf':La(/~r/.s~&~
ct.uA'f/7'e;-e;aj-r(~u~'t.rcest (-<aL7~'e.
(t'hutographier'n grandeurnaturelle.
aceti:lesnoires sont les coupures Liquide:sperm
Les ce!fu)es,bien lignes
deI Organisation of the flow in nearly perfect hexagonal cells (analogy to plant cells).

Rayleigh’s theory
I Rayleigh (1916) proposed the first theory for the linear stability of a steady conductive state in a
gravity field. He showed that a minimum temperature gradient is necessary for the onset of
convection, that depends on several physical parameters.
I Block (1956) showed that the flow in Bénard’s experiments is not driven by gravity but by
temperature–dependence of surface tension, the Marangoni effect. Pearson (1958) developed thecorresponding theory.
I Rayleigh–Bénard convection is still used to denote convection driven by thetemperature–dependence of density in a gravity field while Bénard’s setup is calledBénard–Marangoni.

Outline

Fundamentals of Rayleigh-Bénard convection
Historical backgroundbalance equations and the Boussinesq approximation
Linear stability analysis
Behaviour beyond the onset
High Rayleigh number dynamics and scaling of heat transfer

General principles

In order to pose a fluid dynamical problem, we write:
I Conservations equations: mass, momentum, energy.
I Well established, universal although several level of approximations are possible.
I Boundary conditions (BC): classical ones (Dirichlet, Neumann, Robin) or more exotic (phasechange BC).
I Constitutive equations: Fourier’s law, rheology, equation of state.
I Can be quite complex.
I Generally poorly constrained for the planetary interiors.

Conservation equations
udS

I Consider a fixed control volume.
I The balance equation for a quantity with massdensity f is written:
∂
∂t

Z

Z
ρf dV = −VZ
J f · dS +

S

σf dV
V

where the flux J f and the production σfexpress basic laws of physics.

V

I Use Gauss’ theorem:⇒∂ρf= −∇ · J f + σf∂tConservation equations
udS

I Consider a fixed control volume.
I The balance equation for a quantity with massdensity f is written:
∂
∂t

Z

Z
ρf dV = −VZ
J f · dS +

S

σf dV
V

where the flux J f and the production σfexpress basic laws of physics.

V

I Use Gauss’ theorem:⇒∂ρf= −∇ · J f + σf∂tConservation equations
udS

I Consider a fixed control volume.
I The balance equation for a quantity with massdensity f is written:
∂
∂t

Z

Z
ρf dV = −VZ
J f · dS +

S

σf dV
V

where the flux J f and the production σfexpress basic laws of physics.

V

I Use Gauss’ theorem:⇒∂ρf= −∇ · J f + σf∂tMass conservation

I No production: ⇒ σf = 0
I Convective flow only: J f = ρu
∂
∂t

Z

Z
ρ dV = −

V

S

⇒

Z
ρu · dS ⇒∂ρ
dV = −∂tV∂ρDρ+ u · ∇ρ ≡= −ρ∇ · u∂tDt

I Incompressible flow: ∇ · u = 0.
I NoteZD∂≡+ u · ∇.
Dt∂t

∇ · (ρu) dVVMomentumρDu= −∇P + ∇ · τ + ρg
Dt

Local expression of Newton’s 2nd law with
I forces applied to the surface: pressure P and deviatoric stress τ .
I body forces: gravity ρg.
I Other terms need to be added when dealing with other planetary layers: Coriolis acceleration,
Lorentz force.

Energy conservation

First principle of thermodynamics leads to:ρDe= −∇ · q − P∇ · u + τ : ∇u + ρh
DtIncludesI Viscous dissipation: τ : ∇u
I Radiogenic or tidal heat production: ρhEntropy balance

I Internal energy e is developed as function of two state variables s and ρ (add composition if
necessary).
De
Ds
P Dρ
de = T dS − P dV →
=T
+ 2
Dt
Dt
ρ Dt
I Combine the equation for internal energy:
ρTDs= −∇ · q + τ : ∇u + ρh
Dt

I In the generic form of a conservation equation:
ρ

Ds
q
−1
τ : ∇u + ρH
= −∇ ·+ 2 q · ∇T +Dt
TT
T
| {z } |{z}echangeproduction≥0Equation for the temperature

Depending on the choice of state variable, (T , P) or (T , ρ):
ρCp

DP
DT
= −∇ · q + αT
+ τ : ∇u + ρh
Dt
Dt

DT
= −∇ · q + αTKT ∇ · u + τ : ∇u + ρh
Dt
In the case of incompressibility (Boussinesq approximation, see below), the two equations become
identical.
ρCV

Mechanical boundary condition: dynamic topography
Ricard et al. (2014)
I Free surface BC applied on the topography z = h(x, y, t):σ(x, y, h(x, y, t), t) · n̂ = 0.

(1)

I Hydrostatic balance nearly holds⇒ τ = O(ρgh)

andh = O(α∆T d).

⇒ develop linearly eq. (1):σ(x, y, 0, t) · n̂ + ρg · n̂h(x, y, t) = 0.

(2)

I The slope of the topography is O(α∆Td/L) with L the horizontal wavelength of convection. It issmall enough that equation (2) can be simplified toσzz (x, y, 0, t) = −ρgh(x, y, t),(3)τzx (x, y, 0, t) = τzx (x, y, 0, t) = 0.

(4)

Mechanical boundary condition: dynamic topography
Ricard et al. (2014)
I Free surface BC applied on the topography z = h(x, y, t):σ(x, y, h(x, y, t), t) · n̂ = 0.

(1)

I Hydrostatic balance nearly holds⇒ τ = O(ρgh)

andh = O(α∆T d).

⇒ develop linearly eq. (1):σ(x, y, 0, t) · n̂ + ρg · n̂h(x, y, t) = 0.

(2)

I The slope of the topography is O(α∆Td/L) with L the horizontal wavelength of convection. It issmall enough that equation (2) can be simplified toσzz (x, y, 0, t) = −ρgh(x, y, t),(3)τzx (x, y, 0, t) = τzx (x, y, 0, t) = 0.

(4)

Mechanical boundary condition: dynamic topography
Ricard et al. (2014)
I Free surface BC applied on the topography z = h(x, y, t):σ(x, y, h(x, y, t), t) · n̂ = 0.

(1)

I Hydrostatic balance nearly holds⇒ τ = O(ρgh)

andh = O(α∆T d).

⇒ develop linearly eq. (1):σ(x, y, 0, t) · n̂ + ρg · n̂h(x, y, t) = 0.

(2)

I The slope of the topography is O(α∆Td/L) with L the horizontal wavelength of convection. It issmall enough that equation (2) can be simplified toσzz (x, y, 0, t) = −ρgh(x, y, t),(3)τzx (x, y, 0, t) = τzx (x, y, 0, t) = 0.

(4)

Mechanical boundary condition: dynamic topography
Ricard et al. (2014)
I Free surface BC applied on the topography z = h(x, y, t):σ(x, y, h(x, y, t), t) · n̂ = 0.

(1)

I Hydrostatic balance nearly holds⇒ τ = O(ρgh)

andh = O(α∆T d).

⇒ develop linearly eq. (1):σ(x, y, 0, t) · n̂ + ρg · n̂h(x, y, t) = 0.

(2)

I The slope of the topography is O(α∆Td/L) with L the horizontal wavelength of convection. It issmall enough that equation (2) can be simplified toσzz (x, y, 0, t) = −ρgh(x, y, t),(3)τzx (x, y, 0, t) = τzx (x, y, 0, t) = 0.

(4)

Equation for the topography
Ricard et al. (2014)
I Kinematic definition of the topography, with u h and wh the horizontal and vertical components ofvelocity of the surface:∂h+ u h · ∇h = wh .
∂t
I Neglect ∇h and assume wh = uz (x, y, 0, t)⇒∂h= uz (x, y, 0, t)∂t

I Difficulty: time scale for the evolution of the topography, τη = η/ρgd = O(3 × 103 yr)(post–glacial rebound) very short compared to mantle dynamics time scale.
⇒ Most numerical studies consider instantaneous adjustment of topography and use the free–slipBCs:uz (x, y, 0, t) = 0(5)τzx (x, y, 0, t) = τzx (x, y, 0, t) = 0.

(6)

I Topography is removed from the problem but can be computed as an output by
σzz (x, y, 0, t) = −ρgh(x, y, t).
I In experiments, non–deformable container implies no–slip BC: u = 0.

Equation for the topography
Ricard et al. (2014)
I Kinematic definition of the topography, with u h and wh the horizontal and vertical components ofvelocity of the surface:∂h+ u h · ∇h = wh .
∂t
I Neglect ∇h and assume wh = uz (x, y, 0, t)⇒∂h= uz (x, y, 0, t)∂t

I Difficulty: time scale for the evolution of the topography, τη = η/ρgd = O(3 × 103 yr)(post–glacial rebound) very short compared to mantle dynamics time scale.
⇒ Most numerical studies consider instantaneous adjustment of topography and use the free–slipBCs:uz (x, y, 0, t) = 0(5)τzx (x, y, 0, t) = τzx (x, y, 0, t) = 0.

(6)

I Topography is removed from the problem but can be computed as an output by
σzz (x, y, 0, t) = −ρgh(x, y, t).
I In experiments, non–deformable container implies no–slip BC: u = 0.

Equation for the topography
Ricard et al. (2014)
I Kinematic definition of the topography, with u h and wh the horizontal and vertical components ofvelocity of the surface:∂h+ u h · ∇h = wh .
∂t
I Neglect ∇h and assume wh = uz (x, y, 0, t)⇒∂h= uz (x, y, 0, t)∂t

I Difficulty: time scale for the evolution of the topography, τη = η/ρgd = O(3 × 103 yr)(post–glacial rebound) very short compared to mantle dynamics time scale.
⇒ Most numerical studies consider instantaneous adjustment of topography and use the free–slipBCs:uz (x, y, 0, t) = 0(5)τzx (x, y, 0, t) = τzx (x, y, 0, t) = 0.

(6)

I Topography is removed from the problem but can be computed as an output by
σzz (x, y, 0, t) = −ρgh(x, y, t).
I In experiments, non–deformable container implies no–slip BC: u = 0.

Equation for the topography
Ricard et al. (2014)
I Kinematic definition of the topography, with u h and wh the horizontal and vertical components ofvelocity of the surface:∂h+ u h · ∇h = wh .
∂t
I Neglect ∇h and assume wh = uz (x, y, 0, t)⇒∂h= uz (x, y, 0, t)∂t

I Difficulty: time scale for the evolution of the topography, τη = η/ρgd = O(3 × 103 yr)(post–glacial rebound) very short compared to mantle dynamics time scale.
⇒ Most numerical studies consider instantaneous adjustment of topography and use the free–slipBCs:uz (x, y, 0, t) = 0(5)τzx (x, y, 0, t) = τzx (x, y, 0, t) = 0.

(6)

I Topography is removed from the problem but can be computed as an output by
σzz (x, y, 0, t) = −ρgh(x, y, t).
I In experiments, non–deformable container implies no–slip BC: u = 0.

Equation for the topography
Ricard et al. (2014)
I Kinematic definition of the topography, with u h and wh the horizontal and vertical components ofvelocity of the surface:∂h+ u h · ∇h = wh .
∂t
I Neglect ∇h and assume wh = uz (x, y, 0, t)⇒∂h= uz (x, y, 0, t)∂t

I Difficulty: time scale for the evolution of the topography, τη = η/ρgd = O(3 × 103 yr)(post–glacial rebound) very short compared to mantle dynamics time scale.
⇒ Most numerical studies consider instantaneous adjustment of topography and use the free–slipBCs:uz (x, y, 0, t) = 0(5)τzx (x, y, 0, t) = τzx (x, y, 0, t) = 0.

(6)

I Topography is removed from the problem but can be computed as an output by
σzz (x, y, 0, t) = −ρgh(x, y, t).
I In experiments, non–deformable container implies no–slip BC: u = 0.

Equation for the topography
Ricard et al. (2014)
I Kinematic definition of the topography, with u h and wh the horizontal and vertical components ofvelocity of the surface:∂h+ u h · ∇h = wh .
∂t
I Neglect ∇h and assume wh = uz (x, y, 0, t)⇒∂h= uz (x, y, 0, t)∂t

I Difficulty: time scale for the evolution of the topography, τη = η/ρgd = O(3 × 103 yr)(post–glacial rebound) very short compared to mantle dynamics time scale.
⇒ Most numerical studies consider instantaneous adjustment of topography and use the free–slipBCs:uz (x, y, 0, t) = 0(5)τzx (x, y, 0, t) = τzx (x, y, 0, t) = 0.

(6)

I Topography is removed from the problem but can be computed as an output by
σzz (x, y, 0, t) = −ρgh(x, y, t).
I In experiments, non–deformable container implies no–slip BC: u = 0.

Thermal boundary conditions

I Solids in contact with low viscosity fluids above and/or below that can be considered as well
mixed: uniform temperature.
I Experiments: fluid in contact with a lid. Continuity of temperature and heat flux. In
dimensionless form, it can be written as a Robin BC:
Biθ +

∂θ
=0
∂z

with θ the temperature anomaly and Bi the Biot number.
I Bi → ∞: fixed temperature (Dirichlet BC)
I Bi → 0: fixed flux (Neumann BC)
Reality is often in–between. May apply to the effect of continents on mantle convection (Grigné
et al., 2007a,b).

Constitutive equations 1: Fourier’s law

q = −k∇T
I Second principle:−1
q · ∇T = k
T2



∇T
T

2
≥0⇒k>0I Valid for a very wide range of materials and temperature gradients.
I For crystals, usually anisotropic:
q = −k · ∇T ⇔ qi = −kij ∂j T
This is probably the case in the Earth’s mantle where seismic anisotropy is measured.
I Second principle : eigenvalues of k > 0Constitutive equations 1: Fourier’s law

q = −k∇T
I Second principle:−1
q · ∇T = k
T2



∇T
T

2
≥0⇒k>0I Valid for a very wide range of materials and temperature gradients.
I For crystals, usually anisotropic:
q = −k · ∇T ⇔ qi = −kij ∂j T
This is probably the case in the Earth’s mantle where seismic anisotropy is measured.
I Second principle : eigenvalues of k > 0Constitutive equations 2: Rheology I
I Total stress σ has to be related to the strain rate tensor, e = 21 (∇u + ∇u T ) ≡ (∂j ui + ∂i uj )/2.
Isolating the thermodynamic pressure, P:
σ = −PI + F(e).
I Define the P = σkk /3, the average pressure, and τ , the deviatoric stress as
σ = −PI + τ .
I Newtonian rheology: F is a linear function, i.e.
σij = −Pδij + cijkl ekl .
I For an isotropic fluid (which is not the case of Earth’s mantle, Pouilloux et al., 2007)
σij = −Pδij + λekk δij + 2ηeij .

Constitutive equations 2: Rheology I
I Total stress σ has to be related to the strain rate tensor, e = 21 (∇u + ∇u T ) ≡ (∂j ui + ∂i uj )/2.
Isolating the thermodynamic pressure, P:
σ = −PI + F(e).
I Define the P = σkk /3, the average pressure, and τ , the deviatoric stress as
σ = −PI + τ .
I Newtonian rheology: F is a linear function, i.e.
σij = −Pδij + cijkl ekl .
I For an isotropic fluid (which is not the case of Earth’s mantle, Pouilloux et al., 2007)
σij = −Pδij + λekk δij + 2ηeij .

Constitutive equations 2: Rheology I
I Total stress σ has to be related to the strain rate tensor, e = 21 (∇u + ∇u T ) ≡ (∂j ui + ∂i uj )/2.
Isolating the thermodynamic pressure, P:
σ = −PI + F(e).
I Define the P = σkk /3, the average pressure, and τ , the deviatoric stress as
σ = −PI + τ .
I Newtonian rheology: F is a linear function, i.e.
σij = −Pδij + cijkl ekl .
I For an isotropic fluid (which is not the case of Earth’s mantle, Pouilloux et al., 2007)
σij = −Pδij + λekk δij + 2ηeij .

Constitutive equations 2: Rheology I
I Total stress σ has to be related to the strain rate tensor, e = 21 (∇u + ∇u T ) ≡ (∂j ui + ∂i uj )/2.
Isolating the thermodynamic pressure, P:
σ = −PI + F(e).
I Define the P = σkk /3, the average pressure, and τ , the deviatoric stress as
σ = −PI + τ .
I Newtonian rheology: F is a linear function, i.e.
σij = −Pδij + cijkl ekl .
I For an isotropic fluid (which is not the case of Earth’s mantle, Pouilloux et al., 2007)
σij = −Pδij + λekk δij + 2ηeij .

Constitutive equations 2: Rheology II
I Considering the trace of tensors and the rest, one getsτij = η

∂ui
∂uj
2
+
− δij ∇ · u
∂xj
∂xi
3


P=P−

λ+

2
η
3


,

(7)


∇ · u.

(8)

I The bulk viscosity, ζ = λ + 2η/3 is difficult to measure and usually assumed zero (Stokes
hypothesis), which givesσij = −Pδij + η

∂uj
2
∂ui
+
− δij ∇ · u
∂xj
∂xi
3


.

Constitutive equations 2: Rheology II
I Considering the trace of tensors and the rest, one getsτij = η

∂ui
∂uj
2
+
− δij ∇ · u
∂xj
∂xi
3


P=P−

λ+

2
η
3


,

(7)


∇ · u.

(8)

I The bulk viscosity, ζ = λ + 2η/3 is difficult to measure and usually assumed zero (Stokes
hypothesis), which givesσij = −Pδij + η

∂uj
2
∂ui
+
− δij ∇ · u
∂xj
∂xi
3


.

Equation of state

I Origin of motion: change of density (ρ) with temperature (T ).
⇒ Thermal expansion coefficient:
α=−

1
ρ



∂ρ
∂T


.
P

I Minimal (linear) equation:ρ = ρ0 1 − α (T − T0 ) .
I Effect of pressureρ = ρ0P − P01 − α (T − T0 ) +KTImportant but not leading order since pressure variation is dominated by the hydrostatic, i.e. in
the direction of g. Not considered at first!
I Effect of composition: needs additional parameters such as the FeO mass fraction for the Earthmantle.

Equation of state

I Origin of motion: change of density (ρ) with temperature (T ).
⇒ Thermal expansion coefficient:
α=−

1
ρ



∂ρ
∂T


.
P

I Minimal (linear) equation:ρ = ρ0 1 − α (T − T0 ) .
I Effect of pressureρ = ρ0P − P01 − α (T − T0 ) +KTImportant but not leading order since pressure variation is dominated by the hydrostatic, i.e. in
the direction of g. Not considered at first!
I Effect of composition: needs additional parameters such as the FeO mass fraction for the Earthmantle.

Equation of state

I Origin of motion: change of density (ρ) with temperature (T ).
⇒ Thermal expansion coefficient:
α=−

1
ρ



∂ρ
∂T


.
P

I Minimal (linear) equation:ρ = ρ0 1 − α (T − T0 ) .
I Effect of pressureρ = ρ0P − P01 − α (T − T0 ) +KTImportant but not leading order since pressure variation is dominated by the hydrostatic, i.e. in
the direction of g. Not considered at first!
I Effect of composition: needs additional parameters such as the FeO mass fraction for the Earthmantle.

The Oberbeck–Boussinesq approximation

I Boussinesq (1903) and Oberbeck (1879) propose to simplify the full equations by setting thedensity constant in all terms but the buoyancy term.
I At the same level of approximation (more on that later) the dissipation is negligible andCp = Cv ≡ C .
I The minimal set of equations for convection are (neglecting internal heating for now)
∇ · u = 0,Du
ρ0= −∇P + ρg + η∇2 u,DtDT= κ∇2 T ,Dt
ρ = ρ0 1 − α(T − T0 ) ,I and boundary conditions.

(9)
(10)
(11)
(12)

Outline

Fundamentals of Rayleigh-Bénard convection
Historical backgroundbalance equations and the Boussinesq approximation
Linear stability analysis
Behaviour beyond the onset
High Rayleigh number dynamics and scaling of heat transfer

Generalities

I The problem (the equations) always admit several solutions, notably a motionless steady
conduction solution⇒ What controls the onset of motion? The (in–)stability of the steady conduction solution.
I What forms do the solutions take with motion?

Dimensional analysis: the Rayleigh number

T=T0

I Buoyancy: ρgα∆T ∼ ρv/τc ∼ ρd/τc2 .
d

T=T0+ΔT⇒ Convective time: τc2 = d/gα∆T .
I Diffusive time: τd = d 2 /κ.
I Viscous time: τv = d 2 /ν = ρd 2 /ηI Convection if τv τd /τc2 >> 1
Ra ≡α∆Tgd 3
> Rc ∼ 103
κν

Gross estimate for Earth’s mantle: Ra ∼ 108

Dimensional analysis: the Rayleigh number

T=T0

I Buoyancy: ρgα∆T ∼ ρv/τc ∼ ρd/τc2 .
d

T=T0+ΔT⇒ Convective time: τc2 = d/gα∆T .
I Diffusive time: τd = d 2 /κ.
I Viscous time: τv = d 2 /ν = ρd 2 /ηI Convection if τv τd /τc2 >> 1
Ra ≡α∆Tgd 3
> Rc ∼ 103
κν

Gross estimate for Earth’s mantle: Ra ∼ 108

Dimensional analysis: the Rayleigh number

T=T0

I Buoyancy: ρgα∆T ∼ ρv/τc ∼ ρd/τc2 .
d

T=T0+ΔT⇒ Convective time: τc2 = d/gα∆T .
I Diffusive time: τd = d 2 /κ.
I Viscous time: τv = d 2 /ν = ρd 2 /ηI Convection if τv τd /τc2 >> 1
Ra ≡α∆Tgd 3
> Rc ∼ 103
κν

Gross estimate for Earth’s mantle: Ra ∼ 108

Dimensional analysis: the Rayleigh number

T=T0

I Buoyancy: ρgα∆T ∼ ρv/τc ∼ ρd/τc2 .
d

T=T0+ΔT⇒ Convective time: τc2 = d/gα∆T .
I Diffusive time: τd = d 2 /κ.
I Viscous time: τv = d 2 /ν = ρd 2 /ηI Convection if τv τd /τc2 >> 1
Ra ≡α∆Tgd 3
> Rc ∼ 103
κν

Gross estimate for Earth’s mantle: Ra ∼ 108

Perturbation equations
The system of equation admits a motionless (u = 0) steady (∂t = 0) conduction solution:∇P = ρg(13)

2∇ T =0

(14)









ρ = ρ0 1 − α(T − T0 )(15)

T (d) = T0 and T (0) = T0 + ∆Tz
⇒ Tc = T0 + ∆T − ∆Td(16)

ρc = ρ0 1 − α∆T 1 −

z
d



⇒ Pc = ...

(17)
(18)

Write equations for the perturbations of the steady conduction solution, θ = T − Tc , p = P − Pc :∇·u =0Duρ0= −∇p − ρ0 αθg + η∇2 u
Dt
Dθ
∆T
=
uz + κ∇2 θ
Dt
d

(19)
(20)
(21)

Dimensionless equations
I There are several ways of doing it but I choose here
x 0, y0 =

pd 2
z
1
θ
κt
x, y 0; z = + ; θ0 =; t 0 = 2 ; p0 =
d
d
2
∆T
d
κη

I We get, after dropping the 0 s:∇·u =01 Du= −∇p + ∇2 u + RaθẑPr Dt
Dθ= uz + ∇2 θ
Dt

(22)
(23)
(24)

I withρ0 gα∆Td 3
the Rayleigh number
κη
η
Pr =
the Prandtl number
ρ0 κ

Ra =

I and boundary conditions at z = ±1/2 (free–slip for now)
θ = 0; uz = 0; ∂z ux = ∂z uy = 0 ⇒ ∂z2 uz = 0.

(25)
(26)

Dimensionless equations
I There are several ways of doing it but I choose here
x 0, y0 =

pd 2
z
1
θ
κt
x, y 0; z = + ; θ0 =; t 0 = 2 ; p0 =
d
d
2
∆T
d
κη

I We get, after dropping the 0 s:∇·u =01 Du= −∇p + ∇2 u + RaθẑPr Dt
Dθ= uz + ∇2 θ
Dt

(22)
(23)
(24)

I withρ0 gα∆Td 3
the Rayleigh number
κη
η
Pr =
the Prandtl number
ρ0 κ

Ra =

I and boundary conditions at z = ±1/2 (free–slip for now)
θ = 0; uz = 0; ∂z ux = ∂z uy = 0 ⇒ ∂z2 uz = 0.

(25)
(26)

Dimensionless equations
I There are several ways of doing it but I choose here
x 0, y0 =

pd 2
z
1
θ
κt
x, y 0; z = + ; θ0 =; t 0 = 2 ; p0 =
d
d
2
∆T
d
κη

I We get, after dropping the 0 s:∇·u =01 Du= −∇p + ∇2 u + RaθẑPr Dt
Dθ= uz + ∇2 θ
Dt

(22)
(23)
(24)

I withρ0 gα∆Td 3
the Rayleigh number
κη
η
Pr =
the Prandtl number
ρ0 κ

Ra =

I and boundary conditions at z = ±1/2 (free–slip for now)
θ = 0; uz = 0; ∂z ux = ∂z uy = 0 ⇒ ∂z2 uz = 0.

(25)
(26)

The Prandtl number

I Characteristics of the working fluid
I Liquid water: Pr ∼ 7
I Earth’s mantle: Pr ∼ 1025
I Water ice: Pr ∼ 1017
⇒ Inertia term negligible for convection in solids!
I Kinetic energy of Earth’s mantle (mass 4 × 1024 kg), assuming a mean velocity of 3 cm/yr is∼ 2 × 106 J. Similar to a car driving at 100 km/hr.
⇒ The Prandtl number is taken as infinite.

Mode decomposition for the linear problem I
I Considering infinitely small perturbations of the conduction solution, the problem can be linearised:∇·u =01 ∂u= −∇p + ∇2 u + RaθẑPr ∂t∂θ= uz + ∇2 θ
∂t

(27)
(28)
(29)

I The perturbation can be developed in time–dependent Fourier modes and, for a linear problem,
each mode can be analysed independently. The problem is independent of the horizontal
orientation and we choose:(θ, p, ux , uz ) = (Θ(z), P(z), U (z), W (z))eσt eikx .
I If <(σ) > 0 the instability grows.
I The conduction solution is stable if all modes of perturbation have <(σ) < 0Mode decomposition for the linear problem IId
I Denoting D ≡ dz

ikU + DW = 0,

(30)

 i

(31)

h

2

Pr −ikP + D −kh2

U = σU ,i

Pr − DP + D2 −k 2 W + RaΘ = σW ,

2

W + D −k2

Θ = σΘ

(32)(33)

I This is a generalised eigenvalue problem of the form:
L · X = σR · X(34)

with X = (P; U ; W ; Θ)T the global vertical mode and R a diagonal matrix with 0 or 1 on the
diagonal.
I with boundary conditions for the free–slip case


W

±

1
2




= 0; DU

±

1
2




= 0; Θ ±

1
2


=0

The linear operator is self–adjoint
I Define the dot product of two modes asZ
hX 2 |X 1 i =ei(k1 −k2 )x

Z 1

2P¯2 P1 + U¯2 U1 + W̄2 W 1 + RaΘ̄2 Θ1 dz dx





1
−2

with Ū the complex–conjugate of U , etc.
I Then a series of integrations by part, using the boundary conditions, allows you to show (Schlüter
et al., 1965) thathX 2 |LX 1 i = hLX 2 |X 1 imeaning that the Linear problem is self–adjoint. The R operator is also self–adjoint.
⇒ All the eigenvalues are real:hX|LXi = hX|σRXi = σ hX|RXi
hLX|Xi = hσRX|Xi = σ̄ hRX|Xi⇒ σ = <(σ)
I The instability is characterised by the change of sign of σ. We just need to search for σ = 0(neutral stability).

Solution for free–slip BCs

I When both boundaries are free–slip,
W = cos(πz) provides the solution.
I Neutral stability:Rac =π2 + k 2

3

k2

I Minimum valueRc =I First unstable mode has wavelength
λc =⇒ rolls√
2 wider than they are tall.

√2π= 2 2.
kc

27π 4
π' 657 for kc = √42

Numerical technique for more general BCs
I Rigid BC: uz = ux = 0. Since ∂x ux + ∂z uz = 0, this implies that ∂z uz = 0 at the boundary.
I For at least one rigid BC, solution has to be computed numerically, for example by series
expansion of the vertical mode.
I A simple way: use differentiation matrices with Chebyshev colocation pseudo-spectral method(e.g. Guo et al., 2012).
I A function F (z) defined on [−1, 1] is discretised on Chebyshev–Gauss–Lobatto nodal points
zi = cos(iπ/N ), for i = 0..N , to get vector F ≡ (Fi )i=0..N .
I The n th derivative at the same points is simply obtained using the corresponding differentiation
matrix D(n) :F (n) = D(n) · F
I The linear operator L is written as a block matrix where the D operator is replaced by thedifferentiation matrix.
I Lines corresponding to the boundaries are replaced by corresponding BCs, with a 0 on the
diagonal of R, or simply removed for Dirichlet BCs.
I Use an eigenvector–eigenvalue solver (e.g. in Python) to find the eigenvalue as function of Ra.
Differentiation matrices are available in Python at https://github.com/labrosse/dmsuite.
I Procedure:
I For a given wavenumber k, find the value Ra(k) that makes σ = 0.
I Find the value of k that minimizes Ra(k)Numerical technique for more general BCs
I Rigid BC: uz = ux = 0. Since ∂x ux + ∂z uz = 0, this implies that ∂z uz = 0 at the boundary.
I For at least one rigid BC, solution has to be computed numerically, for example by series
expansion of the vertical mode.
I A simple way: use differentiation matrices with Chebyshev colocation pseudo-spectral method(e.g. Guo et al., 2012).
I A function F (z) defined on [−1, 1] is discretised on Chebyshev–Gauss–Lobatto nodal points
zi = cos(iπ/N ), for i = 0..N , to get vector F ≡ (Fi )i=0..N .
I The n th derivative at the same points is simply obtained using the corresponding differentiation
matrix D(n) :F (n) = D(n) · F
I The linear operator L is written as a block matrix where the D operator is replaced by thedifferentiation matrix.
I Lines corresponding to the boundaries are replaced by corresponding BCs, with a 0 on the
diagonal of R, or simply removed for Dirichlet BCs.
I Use an eigenvector–eigenvalue solver (e.g. in Python) to find the eigenvalue as function of Ra.
Differentiation matrices are available in Python at https://github.com/labrosse/dmsuite.
I Procedure:
I For a given wavenumber k, find the value Ra(k) that makes σ = 0.
I Find the value of k that minimizes Ra(k)Numerical technique for more general BCs
I Rigid BC: uz = ux = 0. Since ∂x ux + ∂z uz = 0, this implies that ∂z uz = 0 at the boundary.
I For at least one rigid BC, solution has to be computed numerically, for example by series
expansion of the vertical mode.
I A simple way: use differentiation matrices with Chebyshev colocation pseudo-spectral method(e.g. Guo et al., 2012).
I A function F (z) defined on [−1, 1] is discretised on Chebyshev–Gauss–Lobatto nodal points
zi = cos(iπ/N ), for i = 0..N , to get vector F ≡ (Fi )i=0..N .
I The n th derivative at the same points is simply obtained using the corresponding differentiation
matrix D(n) :F (n) = D(n) · F
I The linear operator L is written as a block matrix where the D operator is replaced by thedifferentiation matrix.
I Lines corresponding to the boundaries are replaced by corresponding BCs, with a 0 on the
diagonal of R, or simply removed for Dirichlet BCs.
I Use an eigenvector–eigenvalue solver (e.g. in Python) to find the eigenvalue as function of Ra.
Differentiation matrices are available in Python at https://github.com/labrosse/dmsuite.
I Procedure:
I For a given wavenumber k, find the value Ra(k) that makes σ = 0.
I Find the value of k that minimizes Ra(k)Numerical technique for more general BCs
I Rigid BC: uz = ux = 0. Since ∂x ux + ∂z uz = 0, this implies that ∂z uz = 0 at the boundary.
I For at least one rigid BC, solution has to be computed numerically, for example by series
expansion of the vertical mode.
I A simple way: use differentiation matrices with Chebyshev colocation pseudo-spectral method(e.g. Guo et al., 2012).
I A function F (z) defined on [−1, 1] is discretised on Chebyshev–Gauss–Lobatto nodal points
zi = cos(iπ/N ), for i = 0..N , to get vector F ≡ (Fi )i=0..N .
I The n th derivative at the same points is simply obtained using the corresponding differentiation
matrix D(n) :F (n) = D(n) · F
I The linear operator L is written as a block matrix where the D operator is replaced by thedifferentiation matrix.
I Lines corresponding to the boundaries are replaced by corresponding BCs, with a 0 on the
diagonal of R, or simply removed for Dirichlet BCs.
I Use an eigenvector–eigenvalue solver (e.g. in Python) to find the eigenvalue as function of Ra.
Differentiation matrices are available in Python at https://github.com/labrosse/dmsuite.
I Procedure:
I For a given wavenumber k, find the value Ra(k) that makes σ = 0.
I Find the value of k that minimizes Ra(k)Numerical technique for more general BCs
I Rigid BC: uz = ux = 0. Since ∂x ux + ∂z uz = 0, this implies that ∂z uz = 0 at the boundary.
I For at least one rigid BC, solution has to be computed numerically, for example by series
expansion of the vertical mode.
I A simple way: use differentiation matrices with Chebyshev colocation pseudo-spectral method(e.g. Guo et al., 2012).
I A function F (z) defined on [−1, 1] is discretised on Chebyshev–Gauss–Lobatto nodal points
zi = cos(iπ/N ), for i = 0..N , to get vector F ≡ (Fi )i=0..N .
I The n th derivative at the same points is simply obtained using the corresponding differentiation
matrix D(n) :F (n) = D(n) · F
I The linear operator L is written as a block matrix where the D operator is replaced by thedifferentiation matrix.
I Lines corresponding to the boundaries are replaced by corresponding BCs, with a 0 on the
diagonal of R, or simply removed for Dirichlet BCs.
I Use an eigenvector–eigenvalue solver (e.g. in Python) to find the eigenvalue as function of Ra.
Differentiation matrices are available in Python at https://github.com/labrosse/dmsuite.
I Procedure:
I For a given wavenumber k, find the value Ra(k) that makes σ = 0.
I Find the value of k that minimizes Ra(k)Numerical technique for more general BCs
I Rigid BC: uz = ux = 0. Since ∂x ux + ∂z uz = 0, this implies that ∂z uz = 0 at the boundary.
I For at least one rigid BC, solution has to be computed numerically, for example by series
expansion of the vertical mode.
I A simple way: use differentiation matrices with Chebyshev colocation pseudo-spectral method(e.g. Guo et al., 2012).
I A function F (z) defined on [−1, 1] is discretised on Chebyshev–Gauss–Lobatto nodal points
zi = cos(iπ/N ), for i = 0..N , to get vector F ≡ (Fi )i=0..N .
I The n th derivative at the same points is simply obtained using the corresponding differentiation
matrix D(n) :F (n) = D(n) · F
I The linear operator L is written as a block matrix where the D operator is replaced by thedifferentiation matrix.
I Lines corresponding to the boundaries are replaced by corresponding BCs, with a 0 on the
diagonal of R, or simply removed for Dirichlet BCs.
I Use an eigenvector–eigenvalue solver (e.g. in Python) to find the eigenvalue as function of Ra.
Differentiation matrices are available in Python at https://github.com/labrosse/dmsuite.
I Procedure:
I For a given wavenumber k, find the value Ra(k) that makes σ = 0.
I Find the value of k that minimizes Ra(k)Numerical technique for more general BCs
I Rigid BC: uz = ux = 0. Since ∂x ux + ∂z uz = 0, this implies that ∂z uz = 0 at the boundary.
I For at least one rigid BC, solution has to be computed numerically, for example by series
expansion of the vertical mode.
I A simple way: use differentiation matrices with Chebyshev colocation pseudo-spectral method(e.g. Guo et al., 2012).
I A function F (z) defined on [−1, 1] is discretised on Chebyshev–Gauss–Lobatto nodal points
zi = cos(iπ/N ), for i = 0..N , to get vector F ≡ (Fi )i=0..N .
I The n th derivative at the same points is simply obtained using the corresponding differentiation
matrix D(n) :F (n) = D(n) · F
I The linear operator L is written as a block matrix where the D operator is replaced by thedifferentiation matrix.
I Lines corresponding to the boundaries are replaced by corresponding BCs, with a 0 on the
diagonal of R, or simply removed for Dirichlet BCs.
I Use an eigenvector–eigenvalue solver (e.g. in Python) to find the eigenvalue as function of Ra.
Differentiation matrices are available in Python at https://github.com/labrosse/dmsuite.
I Procedure:
I For a given wavenumber k, find the value Ra(k) that makes σ = 0.
I Find the value of k that minimizes Ra(k)Numerical technique for more general BCs
I Rigid BC: uz = ux = 0. Since ∂x ux + ∂z uz = 0, this implies that ∂z uz = 0 at the boundary.
I For at least one rigid BC, solution has to be computed numerically, for example by series
expansion of the vertical mode.
I A simple way: use differentiation matrices with Chebyshev colocation pseudo-spectral method(e.g. Guo et al., 2012).
I A function F (z) defined on [−1, 1] is discretised on Chebyshev–Gauss–Lobatto nodal points
zi = cos(iπ/N ), for i = 0..N , to get vector F ≡ (Fi )i=0..N .
I The n th derivative at the same points is simply obtained using the corresponding differentiation
matrix D(n) :F (n) = D(n) · F
I The linear operator L is written as a block matrix where the D operator is replaced by thedifferentiation matrix.
I Lines corresponding to the boundaries are replaced by corresponding BCs, with a 0 on the
diagonal of R, or simply removed for Dirichlet BCs.
I Use an eigenvector–eigenvalue solver (e.g. in Python) to find the eigenvalue as function of Ra.
Differentiation matrices are available in Python at https://github.com/labrosse/dmsuite.
I Procedure:
I For a given wavenumber k, find the value Ra(k) that makes σ = 0.
I Find the value of k that minimizes Ra(k)Numerical technique for more general BCs
I Rigid BC: uz = ux = 0. Since ∂x ux + ∂z uz = 0, this implies that ∂z uz = 0 at the boundary.
I For at least one rigid BC, solution has to be computed numerically, for example by series
expansion of the vertical mode.
I A simple way: use differentiation matrices with Chebyshev colocation pseudo-spectral method(e.g. Guo et al., 2012).
I A function F (z) defined on [−1, 1] is discretised on Chebyshev–Gauss–Lobatto nodal points
zi = cos(iπ/N ), for i = 0..N , to get vector F ≡ (Fi )i=0..N .
I The n th derivative at the same points is simply obtained using the corresponding differentiation
matrix D(n) :F (n) = D(n) · F
I The linear operator L is written as a block matrix where the D operator is replaced by thedifferentiation matrix.
I Lines corresponding to the boundaries are replaced by corresponding BCs, with a 0 on the
diagonal of R, or simply removed for Dirichlet BCs.
I Use an eigenvector–eigenvalue solver (e.g. in Python) to find the eigenvalue as function of Ra.
Differentiation matrices are available in Python at https://github.com/labrosse/dmsuite.
I Procedure:
I For a given wavenumber k, find the value Ra(k) that makes σ = 0.
I Find the value of k that minimizes Ra(k)The linear operator with free-slip BCs



0:N

0:N

0:N

1:N −1

0

ikI

D

0

 0


−Pr ikI
L =

 0

 −Pr D

0

D



Pr D

(2)

0
2

−k I



0

D
0
0

0



Pr D

(2)



− k2I

I



0:N

 0
0


 1:N −1
0


0
 N

 1:N −1
Pr RaI
D(2) − k 2 I

1:N −1applies to vector X = (P, U , W , Θ)T as
L · X = σR · X.
I For infinite Pr , divide the velocity equations by Pr and make the corresponding diagonal values of
R zero.
I Filter out infinite eigenvalues that are due to the zeros on the diagonal of R.

Effect of mechanical BCs on the linear stability
6000

Rigid–Rigid BCs
Ramin = 1707.76; k = 3.12

Neutral Ra

5000

Free–Rigid BCs
Ramin = 1100.65; k = 2.68
Free–Free BCs
Ramin = 657.52; k = 2.23

4000

3000

2000

1000

1

2

3

4

Wavenumber, k

5

6

7

Outline

Fundamentals of Rayleigh-Bénard convection
Historical backgroundbalance equations and the Boussinesq approximation
Linear stability analysis
Behaviour beyond the onset
High Rayleigh number dynamics and scaling of heat transfer

Weakly non-linear analysis I

I A method pioneered by Malkus and Veronis (1958) for free-slip BCs and Schlüter et al. (1965) for
rigid BCs. See also Manneville (2004).
I The conservation equations are split in their linear and non-linear parts:
L(∂t , ∂x , ∂z , Ra)X = N(X, X),(35)

with X = (p; u; w; θ)T . The linear operator is further developed around the critical Rayleighnumber as
L = Lc − (Ra − Rac )M.
(36)
The solution X and the Rayleigh number are developed as
X = X 1 + 2 X 2 + 3 X 3 + ...
2

Ra = Rac + Ra1 +  Ra2 + ...

(37)
(38)

Weakly non-linear analysis II
and we get a set of equations for the increasing order of :
Lc X 1 = 0,

(39)

Lc X 2 = N(X 1 , X 1 ) + Ra1 MX 1 ,(40)

Lc X 3 = N(X 1 , X 2 ) + N(X 2 , X 1 ) + Ra1 MX 2 + Ra2 MX 1 .

(41)

n−1

Lc X n =X

n−1

N(X l , X n−l ) +l=1

X

Ral MX n−l

(42)

l=1

I Solvability condition for each degree obtained by taking the dot product by X c and using theHermitian property ⇒ values of Rai .
I We can prove (Labrosse et al., 2018) that Ra2n+1 = 0 so that, to leading order,r

Ra − Rac,
Ra2
Ra − Rac
Nu = 1 + A
,
Rac
=

(43)
(44)

with Nu = qd/k∆T the Nusselt number, the dimensionless heat flux.
I In the case of free-slip boundary conditions, the coefficient A = 2 can be determined analytically.

Stability of finite amplitude solutions
I Schlüter et al. (1965) showed that only rolls are stable finite amplitude solutions close to the onsetof convection.
I Busse (1967) showed that a finite range of wavenumber leads to stable roll solution.

(Busse 1967; Busse & Whitehead, 1971)

Example calculation close to onset

I Ra = 800, Pr = ∞.
I Aspect ratio = 32 × 32 × 1.
I Initial condition: conductive solution plus
random noise.
I Pattern dynamics with long–distanceinteractions between defects.
I Steady–state: Rolls at π/4 angle so that a√natural number of 2 2 wavelength fit.
I Method of solution: finite differences andmultigrid (Sotin and Labrosse, 1999).

Experimental test of Busse (1967)’s theory
Busse and Whitehead (1971)(Busse & Whitehead, 1971)

Experimental test of Busse (1967)’s theory
Busse and Whitehead (1971)(Busse & Whitehead, 1971)

Experimental test of Busse (1967)’s theory
Busse and Whitehead (1971)

Fig. 10: Cross-roll instability at R=3000, d=5 mm, α=2π/1.64.Time intervals between
subsequent photographs are 10, 4, 3, 7 and 28 min, respectively.

Experimental test of Busse (1967)’s theory
Busse and Whitehead (1971)

Fig. 11: Zigzag instability at R=3600, d=5 mm, α=2π/2.8.Time intervals between
subsequent photographs are 9, 10, 10, 26 and 72 min, respectively.

Experimental test of Busse (1967)’s theory
Busse and Whitehead (1971)

Fig. 15: Pinching instability at R=18x103, d=1 cm, α1=2π/2.55, α1=2π/1.7.
Time intervals between the photographs is 35 min. (Busse & Whitehead, 1971)

Experiments vs. theory of Busse balloon

I Very good match of observations to thetheory!

Origin of the hexagonal flow

I Many experiments (starting with Bénard’s) lead to hexagonal patterns.
I Hexagonal patterns are non–symmetrical with respect to z → −z tranformation, whereas rolls are.

1

1

0.8

0.8

0.8 1

60

1.2

1.2

1.2

0.6

0.6

20

40

60

80

100

120

140

160

180

200

220

240

20

40

0.6

80

100

120

I Hexagonal flow is obtained for asymmetrical conditions such as provided by depth– or
temperature–dependent properties (e.g. η(T )) or volumetric heat generation (figure).

Outline

Fundamentals of Rayleigh-Bénard convection
Historical backgroundbalance equations and the Boussinesq approximation
Linear stability analysis
Behaviour beyond the onset
High Rayleigh number dynamics and scaling of heat transfer

Example calculation at high Ra
I Ra = 107 , Pr = ∞, aspect ratio 4 × 4 × 1
I Initial conditions: T = 1/2 and exponential variation in thin layers to match BCs plus small
random noise.
I Two iso-temperature surface represented.

Regime diagram in the (Ra, Pr ) space
Krishnamurti (1973)

Ra = 107

Ra = 105(Sotin & Labrosse, 1999)(Krishnamurti, 1973)

Temperature profiles
0.0
0.1

Ra=107
Ra=106
Ra=105Boundary layer

0.2

Depth

0.3
0.4
0.5

I Efficient mixing in the bulk of the domain ⇒uniform temperature.
I Matching the boundary conditions ⇒boundary layers.
I Increasing the Rayleigh number makes the
thickness of boundary layers decrease.

Well mixed interior

0.6
0.7
0.8
0.9
1.0

Boundary layer

0.0

0.2

0.4

0.6

Mean temperature

0.8

1.0

Advection and conduction profiles
I Integrate the energy balance equation between the top boundary and any depth z, averaged over
time:


∂T
1
∂T
qtop ≡ −z==−(z) + uz (T (z) − T ).
∂z
2∂z
I Increase of velocity with Ra makes the advection increase ⇒ thickness of boundary layers
decreases to match the heat flow.
Min temperature
Avg. Temperature
Max temperature

Ra=1061.0

0.8

Radius0.6
Diffusion
Advection
0.4

0.2

0.0
0.00

0.25

0.50
0.75
Temperature

1.00

0

5

10
Heat flux

15

20

Simple dimensional argument for the heat flow

I Dimensionless heat flow Nu = qd/k∆T = f (Ra) = ARaβ to be valid over large range of Ravalues.
I At very large Ra, boundary layers and the resulting plumes get very small.
I The dynamics of convection and the resulting heat flow should become independent of the total
thickness:βk∆T gα∆Td 3
1
⇒β= .
q=A
d
κν
3

Boundary layer instabilities

t=0.00010, q=25.04

t=0.00020, q=18.68

t=0.00040, q=13.74

t=0.00052 q=13.11

t=0.00058 q=15.43

t=0.00066 q=24.67

Howard’s model I
Howard (1964) proposed a simple model for convection at high Ra:
I Based on observations of boundary layers (BL) instabilities.
I Boundary layer grows by diffusion.
I BL becomes unstable when its Rayleigh number reaches a critical value.
I Cycle restarts.
I Ergodicity assumptions: horizontal average on an infinite layer equals time average on one cycle.
I Diffusive part of the cycle takes much longer than destabilisation part ⇒ dominates the time
average.


Z tc
1
∆T
0.5 − z
√
⇒T =
erf
dt
tc 0
2
2 κt
for the top boundary layer and something equivalent at the bottom.
I Heat flux
Z tc
1
k∆T
k∆T
k∆T
√
q=
dt = √=.
tc 0 2 πκtδcπκtcHoward’s model I
Howard (1964) proposed a simple model for convection at high Ra:
I Based on observations of boundary layers (BL) instabilities.
I Boundary layer grows by diffusion.
I BL becomes unstable when its Rayleigh number reaches a critical value.
I Cycle restarts.
I Ergodicity assumptions: horizontal average on an infinite layer equals time average on one cycle.
I Diffusive part of the cycle takes much longer than destabilisation part ⇒ dominates the time
average.


Z tc
1
∆T
0.5 − z
√
⇒T =
erf
dt
tc 0
2
2 κt
for the top boundary layer and something equivalent at the bottom.
I Heat flux
Z tc
1
k∆T
k∆T
k∆T
√
q=
dt = √=.
tc 0 2 πκtδcπκtcHoward’s model I
Howard (1964) proposed a simple model for convection at high Ra:
I Based on observations of boundary layers (BL) instabilities.
I Boundary layer grows by diffusion.
I BL becomes unstable when its Rayleigh number reaches a critical value.
I Cycle restarts.
I Ergodicity assumptions: horizontal average on an infinite layer equals time average on one cycle.
I Diffusive part of the cycle takes much longer than destabilisation part ⇒ dominates the time
average.


Z tc
1
∆T
0.5 − z
√
⇒T =
erf
dt
tc 0
2
2 κt
for the top boundary layer and something equivalent at the bottom.
I Heat flux
Z tc
1
k∆T
k∆T
k∆T
√
q=
dt = √=.
tc 0 2 πκtδcπκtcHoward’s model I
Howard (1964) proposed a simple model for convection at high Ra:
I Based on observations of boundary layers (BL) instabilities.
I Boundary layer grows by diffusion.
I BL becomes unstable when its Rayleigh number reaches a critical value.
I Cycle restarts.
I Ergodicity assumptions: horizontal average on an infinite layer equals time average on one cycle.
I Diffusive part of the cycle takes much longer than destabilisation part ⇒ dominates the time
average.


Z tc
1
∆T
0.5 − z
√
⇒T =
erf
dt
tc 0
2
2 κt
for the top boundary layer and something equivalent at the bottom.
I Heat flux
Z tc
1
k∆T
k∆T
k∆T
√
q=
dt = √=.
tc 0 2 πκtδcπκtcHoward’s model I
Howard (1964) proposed a simple model for convection at high Ra:
I Based on observations of boundary layers (BL) instabilities.
I Boundary layer grows by diffusion.
I BL becomes unstable when its Rayleigh number reaches a critical value.
I Cycle restarts.
I Ergodicity assumptions: horizontal average on an infinite layer equals time average on one cycle.
I Diffusive part of the cycle takes much longer than destabilisation part ⇒ dominates the time
average.


Z tc
1
∆T
0.5 − z
√
⇒T =
erf
dt
tc 0
2
2 κt
for the top boundary layer and something equivalent at the bottom.
I Heat flux
Z tc
1
k∆T
k∆T
k∆T
√
q=
dt = √=.
tc 0 2 πκtδcπκtcHoward’s model I
Howard (1964) proposed a simple model for convection at high Ra:
I Based on observations of boundary layers (BL) instabilities.
I Boundary layer grows by diffusion.
I BL becomes unstable when its Rayleigh number reaches a critical value.
I Cycle restarts.
I Ergodicity assumptions: horizontal average on an infinite layer equals time average on one cycle.
I Diffusive part of the cycle takes much longer than destabilisation part ⇒ dominates the time
average.


Z tc
1
∆T
0.5 − z
√
⇒T =
erf
dt
tc 0
2
2 κt
for the top boundary layer and something equivalent at the bottom.
I Heat flux
Z tc
1
k∆T
k∆T
k∆T
√
q=
dt = √=.
tc 0 2 πκtδcπκtcHoward’s model I
Howard (1964) proposed a simple model for convection at high Ra:
I Based on observations of boundary layers (BL) instabilities.
I Boundary layer grows by diffusion.
I BL becomes unstable when its Rayleigh number reaches a critical value.
I Cycle restarts.
I Ergodicity assumptions: horizontal average on an infinite layer equals time average on one cycle.
I Diffusive part of the cycle takes much longer than destabilisation part ⇒ dominates the time
average.


Z tc
1
∆T
0.5 − z
√
⇒T =
erf
dt
tc 0
2
2 κt
for the top boundary layer and something equivalent at the bottom.
I Heat flux
Z tc
1
k∆T
k∆T
k∆T
√
q=
dt = √=.
tc 0 2 πκtδcπκtcHoward’s model II

I The thickness δc is determined by the instability of the boundary layerRagα∆T δc3
=
Raδ =
2κν
2



I Heat flux then scales as
q=

k∆T
d



Ra
2Rδc

δc
d

3
= Rδc .

1/3.

I As noted by Howard (1964) Rδc is different from the critical value for the stability of the wholelayer since it concerns the destabilisation of a curved profile.
I This theory provides a scaling relation for the fluctuation time of the boundary layer:
tc =

δc2
d2
=
πκ
πκ



2Rδc
Ra

2/3
.

Howard’s model II

I The thickness δc is determined by the instability of the boundary layerRagα∆T δc3
=
Raδ =
2κν
2



I Heat flux then scales as
q=

k∆T
d



Ra
2Rδc

δc
d

3
= Rδc .

1/3.

I As noted by Howard (1964) Rδc is different from the critical value for the stability of the wholelayer since it concerns the destabilisation of a curved profile.
I This theory provides a scaling relation for the fluctuation time of the boundary layer:
tc =

δc2
d2
=
πκ
πκ



2Rδc
Ra

2/3
.

Howard’s model II

I The thickness δc is determined by the instability of the boundary layerRagα∆T δc3
=
Raδ =
2κν
2



I Heat flux then scales as
q=

k∆T
d



Ra
2Rδc

δc
d

3
= Rδc .

1/3.

I As noted by Howard (1964) Rδc is different from the critical value for the stability of the wholelayer since it concerns the destabilisation of a curved profile.
I This theory provides a scaling relation for the fluctuation time of the boundary layer:
tc =

δc2
d2
=
πκ
πκ



2Rδc
Ra

2/3
.

Howard’s model II

I The thickness δc is determined by the instability of the boundary layerRagα∆T δc3
=
Raδ =
2κν
2



I Heat flux then scales as
q=

k∆T
d



Ra
2Rδc

δc
d

3
= Rδc .

1/3.

I As noted by Howard (1964) Rδc is different from the critical value for the stability of the wholelayer since it concerns the destabilisation of a curved profile.
I This theory provides a scaling relation for the fluctuation time of the boundary layer:
tc =

δc2
d2
=
πκ
πκ



2Rδc
Ra

2/3
.

Experiments at very high Ra
Niemela et al. (2000)(Niemela et al, 2000)

I Working fluid: cryogenic helium.
I Pr ∼ 1

β=2/7

e
op

Sl

9

30
0.

β=

β=1/3

β=0.309

I 1 m–high tank.

mean field

I Exponent β close to but different
from 1/3.

References I
Bénard, H. (1900a). Les tourbillons celullaires dans une nappe liquide. Deuxième partie : procédésmécaniques et optiques d’examen. Lois numériques des phénomènes. 12: 1309–1328.
Bénard, H. (1900b). Les tourbillons celullaires dans une nappe liquide. Première partie : descriptiongénérale des phénomènes. 12: 1261–1271.
Bénard, H. (1901). Les tourbillons celullaires dans une nappe liquide transportant la chaleur parconvection en régime permanent. Ann. Chim. Phys. 23: 62–144.
Block, M. J. (1956). Surface Tension as the Cause of Bénard Cells and Surface Deformation in a
Liquid Film. Nature 178: 650–651.
Boussinesq, J. (1903). Théorie Analytique de la Chaleur. Vol. 2. pp. 157–161. Gauthier-Villars.
Busse, F. H. (1967). On the stability of two-dimensional convection in a layer heated from below. J.
Math. and Phys. 46: 140–149.
Busse, F. H. and J. A. Whitehead (1971). Instabilities of convection rolls in a high Prandtl number
fluid. J. Fluid Mech. 47: 305–320.
Grigné, C., S. Labrosse, and P. J. Tackley (2007a). Convection under a lid of finite conductivity in wide
aspect ratio models: effect of continents on the wavelength of mantle flow. J. Geophys. Res. 112:
B08403.

References II
Grigné, C., S. Labrosse, and P. J. Tackley (2007b). Convection under a lid of finite conductivity: Heatflux scaling and application to continents. J. Geophys. Res. 112: B08402.
Guo, W., G. Labrosse, and R. Narayanan (2012). The Application of the Chebyshev-Spectral Methodin Transport Phenomena. Berlin: Springer-Verlag.
Howard, L. N. (1964). “Convection at High Rayleigh number”. Proceedings of the EleventhInternational Congress of Applied Mechanics. Ed. by H. Gortler. New York: Springer-Verlag,
1109–1115.
Krishnamurti, R. (1973). Some further studies on the transition to turbulent convection. J. Fluid
Mech. 60: 285–303.
Labrosse, S., A. Morison, R. Deguen, and T. Alboussière (2018). Rayleigh-Bénard convection in acreeping solid with a phase change at either or both horizontal boundaries. J. Fluid Mech. 846: 5–36.
Malkus, W. and G. Veronis (1958). Finite amplitude cellular convection. J. Fluid Mech. 4: 225–260.
Manneville, P. (2004). Instabilities, Chaos and Turbulence - An introduction to nonlinear dynamics andcomplex systems. London: Imperial College Press.
Niemela, J. J., L. Skrbek, K. R. Sreenivasan, and R. J. Donnelly (2000). Turbulent convection at very
high Rayleigh numbers. Nature 404: 837–841.

References III
Oberbeck, A. (1879). Über die Wärmeleitung des Flüssigkeiten bei Berücksichtigung des Strömungeninfolge von Temperaturdifferenzen. Ann. Phys. Chem. 7: 271–292.
Pearson (1958). On convection cells induced by surface tension. J. Fluid Mech. 4: 489–500.
Pouilloux, L., E. Kaminski, and S. Labrosse (2007). Anisotropic rheology of a cubic medium andimplications for geological materials. Geophysical Journal International 170: 876–885.
Rayleigh, L. (1916). On convection currents in a horizontal layer of fluid, when the higher temperature
is on the under side. Phil. Mag. 32: 529–546.
Ricard, Y., S. Labrosse, and F. Dubuffet (2014). Lifting the cover of the cauldron: Convection in hot
planets. Geochem. Geophys. Geosyst. 15: 4617–4630.
Schlüter, A., D. Lortz, and F. Busse (1965). On the Stability of Steady Finite Amplitude Convection.
J. Fluid Mech. 23: 129–144.
Sotin, C. and S. Labrosse (1999). Three-dimensional Thermal convection of an isoviscous,
infinite–Prandtl–number fluid heated from within and from below: applications to heat transfer inplanetary mantles. Phys. Earth Planet. Inter. 112: 171–190.