Planetary interiors:
Magnetic fields, Convection and Dynamo Theory
2. Convection in planetary interiors
Chris Jones, Department of Applied Mathematics
University of Leeds UK

FDEPS Lecture 2, Kyoto, 28th November 2017

Section 2.
Convection in planetary interiors
2.1 Boussinesq Rayleigh-Bénard convection(2.1) Boussinesq Rayleigh-Bénard convection

2/85

Rayleigh-Bénard convection

Layer of liquid of depth d between two parallel horizontal plates
heated from below.
Hotter fluid expands, becomes less dense, and rises up. At the top,
fluid cools, becomes denser and sinks back to the bottom.
This continual motion carries heat from the bottom of the layer tothe top: heat transport by convection.
If the heating is small, the flow may be steady or vary slowly.
If the heating is strong, the flow is very turbulent.
(2.1) Boussinesq Rayleigh-Bénard convection

3/85

Convection occurs on many different scales

Left image is Bénard’s original experiment: layer depth about
1mm, convecting cells similar size. Fluid flow made visible byaluminium flakes which reflect the light.
Right image is a snapshot of the solar granulation. Each granule ishot fluid rising, and is typically 1000km across. This processtransports the heat through the outer layers.
The Bénard cells are fairly steady, the solar granules are very
turbulent.
(2.1) Boussinesq Rayleigh-Bénard convection

4/85

Boussinesq approximation 1.
The density change that drives the convection is small, so we canignore variations of density except when they are multiplied bygravity.
For a liquid, the density isρ = ρ0 (1 − αT ).

(2.1.1)α is small, so mass conservation is∇ · u = 0.

(2.1.2)

Boussinesq equation of motion is∂u1+ (u · ∇)u = − ∇p − g (1 − αT )ẑ + ν∇2 u.
∂t
ρ0(2.1.3)

which is the Navier-Stokes equation. Here p is the pressure, ν isthe kinematic viscosity (units m2 s−1 ).
(2.1) Boussinesq Rayleigh-Bénard convection

5/85

Boussinesq approximation 2.
The temperature equation is∂T+ u · ∇T = κ∇2 T .
∂t(2.1.4)

Here κ is the thermal diffusivity. It also has units m2 s−1 .
The thermal conductivity of the fluid is K = ρcp κ, where cp is the
specific heat at constant pressure.
This Boussinesq convection problem was discussed by Rayleigh in
1916. A nice treatment is given in Chandrasekhar, 1961.
The Boussinesq equations have a basic state conduction solution,
where the fluid is at rest, and heat is tranported only by thermal
conduction.
As the temperature difference across the layer is increased, at a
certain value the fluid starts to move.
(2.1) Boussinesq Rayleigh-Bénard convection

6/85

Basic State
The basic state solution has zero velocity and temperature fixed intime and a function of z only0=−1
∇p − g (1 − αT )ẑρ0(2.1.5)0 = κ∇2 T .

(2.1.6)andThe solution isT =T̄ = T0 + ∆T (1 − z/d),1
p = p̄ = p0 − ρ0 gz(1 − αT0 − α∆T ) − ρ0 g α∆Tz 2 .
2The heat flux carried by conduction is∆T
F = −K dT
dz = κρcp d Watts per square metre.
(2.1) Boussinesq Rayleigh-Bénard convection(2.1.7)(2.1.8)7/85

2.2 Linear theory of Rayleigh-Bénard convection(2.2) Linear theory of Rayleigh-Bénard convection

8/85

Stability of the basic state: onset of convection
We imagine the bottom plate temperature being slowly increased,
i.e. ∆T slowly increased. There comes a point when the layer
starts to convect.
Close to onset, the velocity will be small, so u · ∇u will be small
compared to ν∇2 u, so we can neglect it. The temperature will
only be slightly different from the basic state, so we write
T = T0 + ∆T (1 − z/d) + θ(2.2.1)

where θ is small compared to ∆T , and p = p̄ + p 0 .
The linearised dimensional equations are then1∂u
= − ∇p 0 + g αθẑ + ν∇2 u,∂t
ρ0(2.2.2)∇ · u = 0,(2.2.3)∂θ= uz ∆T /d + κ∇2 θ.
∂t(2.2.4)(2.2) Linear theory of Rayleigh-Bénard convection

9/85

Non-dimensional equations 1.
Note that the two terms which can be neglected close to onset are(u · ∇)u in the equation of motion, and u · ∇θ in the temperatureequation. They become important again as ∆T is increased.
We now write the equations in dimensionless form. The advantageof doing this is that the key dimensionless parameters appear, so
we know how big the parameter space of our problem is.
Sometimes though it is useful to have dimensional equations as
well, as then we can more easily see which terms are important ingeophysical applications.
There are several possible choices, all lead to the same set of
dimensionless parameters, but they appear in different places in theequations.

(2.2) Linear theory of Rayleigh-Bénard convection

10/85

Non-dimensional equations 2.
Our choice is
κ ∂
∂
= 2 ,
∂t
d ∂ t̃

∇=

1˜
∇,
d

u=

ν
ũ,
d

θ=

ν∆T
θ̃,
κ

p0 =

ρ0 ν 2 0p̃
d2(2.2.5)giving1 ∂ũ
˜ 2 ũ,= −∇p̃ 0 + Raθ̃ẑ + ∇Pr ∂ t̃∇ · ũ = 0,∂ θ̃˜ 2 θ̃.
= ũz + ∇∂ t̃

(2.2.6, 7)(2.2.8)Hereg α∆Td 3
ν, Pr = ,(2.2.9)κνκare the two dimensionless parameters called the Rayleigh number
and the Prandtl number.
Ra =(2.2) Linear theory of Rayleigh-Bénard convection

11/85

Non-dimensional equations 3.

At this stage it is conventional to drop the tildes in the equations.
This sometimes confuses students, because the untilded variables
might be dimensionless or might have dimensions. But it is too
tiresome to keep writing tildes, so you just have to be careful tocheck whether the variables are dimensionless or not. So thedimensionless linearised equations are
1 ∂u= −∇p 0 + Raθẑ + ∇2 u,
Pr ∂t(2.2.10)∇ · u = 0,(2.2.11)∂θ= uz + ∇2 θ.
∂t(2.2.121)(2.2) Linear theory of Rayleigh-Bénard convection

12/85

How big is the Rayleigh number?
The Rayleigh number isg α∆Td 3.
κν
In a Bénard experiment with water, in S.I units (metres, seconds,
kg, ◦ K) g ≈ 10, α ≈ 10−4 , ∆T ≈ 1, d ≈ 10−3 , κ ≈ 10−5 ,ν ≈ 10−4 , for water.
Ra =

10−12
≈ 1000,10−9which as we will see is around the critical value for the onset of
convection.
The Prandtl number ν/κ for water is usually taken as 6.8, for air it
is about 0.7, and it is small for liquid metals.
For solar material, ∆T and d are very large, so the Rayleighnumber is enormous, and it is very large in geophysicalapplications. So in practice Ra is massively supercritical. Is linear
theory relevant in planets?
Ra =(2.2) Linear theory of Rayleigh-Bénard convection

13/85

Boundary conditions

The temperature is fixed at z = 0 and the top boundary z = 1. So
we need θ = 0 at z = 0 and z = 1.
The fluid cannot penetrate the boundaries, so we must haveuz = 0 at z = 0 and z = 1.
The natural two remaining boundary conditions would be to set thehorizontal velocities ux and uy zero at z = 0 and z = 1. From thecontinuity equation (8) this gives duz /dz = 0 at z = 0 and z = 1.
However, Rayleigh noticed that if we take stress-free boundary
conditions, dux /dz = duy /dz = 0, this gives d 2 uz /dz 2 = 0 on theboundaries and then there is a simple exact solution, so people
often use stress-free boundary conditions (illustrative boundary
conditions) to get an initial idea of the behaviour.

(2.2) Linear theory of Rayleigh-Bénard convection

14/85

Solving the linear equations
We take the curl of the equation of motion (2.2.10). Recall thatζ = ∇ × u is the vorticity of the fluid, and
ζz =

∂uy
∂ux
−
.
∂x
∂y

(2.2.13)

The curl of (2.2.10) gives the vorticity equation, and the z
component gives1 ∂ζz= ∇ 2 ζz .
(2.2.14)Pr ∂t
The buoyancy force only generates horizontal vorticity, not verticalvorticity.
An important point is that ∇ × ∇p = 0, so taking the curl
eliminates the awkward pressure term.
(2.2.13) is a diffusion equation, so if the vertical vorticity is zero atthe boundaries, it diffuses away to zero everywhere, ζz = 0.
(2.2) Linear theory of Rayleigh-Bénard convection

15/85

The double curl equation
We take the curl of the vorticity equation, which is the double curlof the equation of motion (2.2.10).
∂u
1
∇×∇×
= Ra∇ × ∇ × θẑ + ∇ × ∇ × ∇2 uPr∂t

(2.2.15)

Using ∇ · u = 0 and the vector identity curl curl = grad div - del2 ,we get∂2∂2+.
∂x 2 ∂y 2(2.2.16)
The double curl procedure looks complicated at first sight, but it isoften the best way to solve the convection equations even if there
are nonlinear equations and additional affects like rotation,
magnetic fields or compressibility.
1 ∂ 2∇ uz = Ra∇2H θ + ∇4 uz ,Pr ∂t(2.2) Linear theory of Rayleigh-Bénard convection

where

∇2h =

16/85

The linear stability ODE
With the θ equation (2.2.11)∂θ= uz + ∇2 θ,∂t(2.2.17)

we now have just two scalar unknowns, uz and θ.
We look for solutions which are periodic in x and y , and havetime-dependence exp σt , uz = W (z)exp[i(kx x + ky y ) + σt],θ = Θ(z)exp[i(kx x + ky y ) + σt] and putting these into (2.2.15)
and (2.2.16), letting a2 = kx2 + ky2 ,

d2
2
−
a
W +Ra a2 W = 0.
dz 2(2.2.18)
This is a 6th order system of ODE’s in z, so we need 6 boundary
conditions.

d2
σ− a2 −
dz 2
Pr



d2
− a2 − σdz 2(2.2) Linear theory of Rayleigh-Bénard convection



17/85

Equation for the growth rate σWith stress free boundary conditions, W = 0 and D 2 W = 0,
D = d/dz, and then Θ = 0 is equivalent to D 4 uz = 0 since
 2

 2
d
d
σ
2
2
−
a
−
−
a
W = Ra a2 Θ.
(2.2.19)
dz 2
Pr
dz 2
So the solution is W (z) = sin nπz, for integer n, giving
σ 2 + (n2 π 2 + a2 )σ(1 + Pr ) + (n2 π 2 + a2 )2 Pr −a2 RaPr= 0.
(n2 π 2 + a2 )(2.2.20)

This gives a steady solution σ = 0 when
Ra = Rac =(n2 π 2 + a2 )3,
a2(2.2.21)and if Ra is greater than this critical value, there is a solution with
σ > 0 so the disturbance grows exponentially.
(2.2) Linear theory of Rayleigh-Bénard convection

18/85

Dispersion relation for Ra

If Ra < Rac , σ has negative real part, so the disturbance decaysaway with time.
The values of a and n which minimise the critical Rayleigh number
are the important critical values. n = 1 is always the minimum.
√
There is a minimum at ac = π/ 2, where Rac = 27π 4 /4 ≈ 657,
Below this value all disturbances decay, while if Ra > Rac a bandof wavenumbers grow exponentially.
(2.2) Linear theory of Rayleigh-Bénard convection

19/85

No-slip boundary caseIn the more realistic no-slip case the boundary conditions are
W = DW = Θ = 0 on z = 0, 1, and then sin πz isn’t a solution.
However, the 6th order ODE can be integrated numerically with
the no-slip boundary conditions, and the critical Rayleigh number
plot looks similar.
In this case, Rac ≈ 1707.76 and the minimising a ≈ 3.117, close toπ.
Lots of other cases are possible, e.g. one no-slip, one stress-freeboundary, or fixed flux boundaries dΘ/dz = 0.

(2.2) Linear theory of Rayleigh-Bénard convection

20/85

Pattern selection in convection

If Ra is justq above Rac , a band of wavenumbers neara = ac = kx2 + ky2 can grow.
Initially there will be small random perturbations due toimperfections in the experiment, which will be made up of allwavenumbers. Only the unstable ones will grow, the others will
decay, so the growing disturbance will contain only the unstablewavenumbers.
However, this does not give a unique pattern. There could be
two-dimensional rolls aligned along the x axis, with ky = ac andkx = 0, or aligned along the y axis, with kx = ac and ky = 0.
√
Could also have square cells with kx = ky = ac / 2. Even possible
to have hexagonal cells.

(2.2) Linear theory of Rayleigh-Bénard convection

21/85

Some observed convection patterns

Left: Two dimensional convection rolls. Middle: Square cells.
Right: hexagonal cells.
Which mode actually occurs depends on the nonlinear terms.
Although very small, over a very long time scale a preferredpattern generally appears.
Which pattern depends on the boundary conditions, and theRayleigh and Prandtl numbers.

(2.2) Linear theory of Rayleigh-Bénard convection

22/85

Stably stratified fluid

What happens if ∆T < 0, so top boundary is hotter than lower
boundary?
Ra < 0, and we have stably stratified fluid.
The limit Ra → −∞ corresponds to the small diffusion limit, and(2.2.20) becomesa2 RaPr.
(2.2.22)σ2 = 2 2n π + a2
Restoring the dimensions,a2 d 2,σ 2 = −g αβ 2 2n π + a2 d 2(2.2.23)

where β = −∆T /d, so β is now a positive temperature gradient.

(2.2) Linear theory of Rayleigh-Bénard convection

23/85

Internal gravity waves
Taking the square root of this negative quantity,
σ = ±i (g αβ)1/2 √adn 2 π 2 + a2 d 2.

(2.2.24)Since σ is now imaginary, the time dependence exp σt corresponds
to waves.
These waves are called internal gravity waves. The quantity
N = (g αβ)1/2
is called the buoyancy frequency. If there is some viscous andthermal diffusion, then the internal gravity waves are damped.
Physically, in a stably stratified layer, if a parcel of fluid moves
upwards, it finds itself in less dense surroundings, so it falls back.
The buoyancy is now a restoring force.
(2.2) Linear theory of Rayleigh-Bénard convection

24/85

2.3 Nonlinear Rayleigh-Bénard convection(2.3) Nonlinear Rayleigh-Bénard convection

25/85

Nonlinear convection

As the Rayleigh number isincreased, the horizontally
averaged temperature profilechanges.
Boundary layers appear, across which the temperature varies
rapidly. The interior temperature becomes nearly constant.
The heat transport is mainly by conduction in the boundary layers,and by convection in the bulk interior.
If the boundary layer has thickness δ, the heat flux transported is
F =K

∆T
∆T
>> Fbasic = K2δd(2.3.1)The dimensionless ratio F /Fbasic = Nu, the Nusselt number.
(2.3) Nonlinear Rayleigh-Bénard convection

26/85

Malkus-Howard convection theory
At large Rayleigh number, Nu ≈ d/2δ.
The Rayleigh number (2.2.9) isg α∆Td 3.
κνMalkus-Howard idea is that the thickness of the boundary layer
makes it close to critical, soRa =Rabl =

g α∆T δ 3
≈ 103 .
κν

(2.3.2)SoRa =g α∆T δ 3
g α∆Td 3≈
8Nu 3 ≈ 8 × 103 Nu 3 ,κνκν

(2.3.3)Nu ∼ Ra1/3 .

(2.3.4)

givingExperiments suggest that the 1/3 law is close, though not all
experiments agree.
(2.3) Nonlinear Rayleigh-Bénard convection

27/85

Averaged heat flux
The temperature equation (2.1.4) is∂T+ u · ∇T = κ∇2 T .
∂tAverage over a horizontal plane at z,
1
haih =Lx Ly

Z Lx /2 Z Ly /2
a dxdy ,−Lx /2

Lx → ∞, Ly → ∞. (2.3.5)−Ly /2∂ hT ih+ h∇ · (uT ) − T ∇ · uih = κ
∂t

 2

∂ T
∂2T
∂2T
+
+
,
∂x 2
∂y 2
∂z 2 hbut averages of x or y derivatives are zero, and ∇ · u = 0. So


 2 
∂ hT ih∂∂ T+(uz T ) = κ∂t∂z∂z 2 hh(2.3) Nonlinear Rayleigh-Bénard convection

28/85

Convective and conductive heat flux
We now take the time-average, so the ∂/∂t term goes,

 2 

∂ T
∂
(uz T ) = κ∂z∂z 2 hhtrue at every z level. Integrate from 0 to z 0 to get




∂T
∂T
−κ
h(uz T )iz 0 = κ∂z z 0∂z 0∂Th(uz T )iz 0 − κ= F,∂z z 0(2.3.6)(2.3.7)

because F is the heat flux conducted in at the bottom of the layer.
< (uz T ) >z 0 is the convected heat flux at level z 0 and
00−κ < ∂T∂z >z is the conducted heat flux at level z .
Note the convected heat flux is positive because hot fluid rises andcold fluid sinks.
(2.3) Nonlinear Rayleigh-Bénard convection

29/85

Dissipation integral 1.
Equation of motion (2.1.3) is∂u1+ (u · ∇)u = − ∇p 0 + g αTẑ + ν∇2 u,∂tρ0taking out the hydrostatic part of the pressure. We take the scalar
product of this with u and average over x, y , z, and t∂u∂(u2 /2)=, u·(u·∇)u = ∇·(u(u2 /2)), u·∇p 0 = ∇·(p 0 u),∂t∂t
and all these terms vanish when we average. For the last term


∂
∂ui
∂ui ∂ui∂ 2 ui=ui−,(2.3.8)
ui
2
∂xj
∂xj
∂xj ∂xj∂xj
u·givingZ d
gα
00

huz T iz 0 dz = ν(2.3) Nonlinear Rayleigh-Bénard convection

Z d
0∂ui ∂ui∂xj ∂xjdz 0 .

(2.3.9)z030/85

Dissipation integral 2.
We now integrate the heat flux relation across the layer,


Z d
∂T
dz 0 = g αFd
gαh(uz T )iz 0 − κ
∂z
0
z0
or

Z d *
ν
0

∂ui
∂xj

2 +

dz 0 = g αFd − κg α∆T .

(2.3.10)(2.3.11)z0

Now
Fd = Nuκ∆T ,g α∆T = Raκν/d 3

So the dissipation integral can be written +Z d *
∂ui 2
κ2 Ra(Nu − 1)
dz 0 =
.
∂xj
d3
0
0

(2.3.12)z(2.3) Nonlinear Rayleigh-Bénard convection

31/85

Grossmann-Lohse theory 1.

This method is based on estimating where the majority of thedissipation comes from. It may be in the boundary layers, or it may
be in the bulk.
In their picture there is a large scale flow with typical velocity U
which develops a viscous boundary layer.

(2.3) Nonlinear Rayleigh-Bénard convection

32/85

Grossmann-Lohse theory 2.

Depending on the Prandtl number, the viscous boundary layer will
be inside or outside the thermal boundary layer. (2.3.12) gives
Z dνh∇u : ∇ui dz =0

ν3
RaNuPr −2 ,d3

Nu >> 1.

(2.3.13)

If the dissipation is mainly in the boundary layers, we analyse theseboundary layers to estimate ν h∇u : ∇ui. Different estimates also
possible if dissipation is in the bulk.

(2.3) Nonlinear Rayleigh-Bénard convection

33/85

Pr > 1 case, dissipation in boundary layers 1.
Ui
ν

δi
δth
i
δνi UithδiPr > 1 means the thermal boundary layer
is nested inside the viscous boundary layer.
U is large scale flow speed in bulk, δ ν , δ th
viscous and thermal b.l. thicknesses.

Velocity at edge of thermal b.l. = δ th U/δ ν .
Advection balances diffusion in the boundary layers, (2.1.3) is∂u1+ (u · ∇)u = − ∇p − g (1 − αT )ẑ + ν∇2 u.
∂t
ρ0
In the boundary layers dominant horizontal terms are (u · ∇)u andνd 2 u/dz 2 .
This means that

(2.3) Nonlinear Rayleigh-Bénard convection

U2
νU∼ ν 2.
d(δ )(2.3.14)34/85

Pr > 1 case, dissipation in boundary layers 2.
Advection balances diffusion in the boundary layers in thetemperature equation too, (2.1.4),∂T+ u · ∇T = κ∇2 T∂tIn the boundary layers dominant terms here are u · ∇T andκd 2 T/dz 2 . This gives
UT
κT∼ th 2d(δ )givingδ th UT
κT∼ th 2 ,
δν d
(δ )

U2
νU
∼ ν 2,d(δ )(2.3.15)

which givesδ ν /δ th = Pr 1/3 (2.3.16),d/δ ν = (Ud/ν)1/2 = Re 1/2 , (2.3.17)

where Re is the Reynolds number of the flow.
(2.3) Nonlinear Rayleigh-Bénard convection

35/85

Pr > 1 case, dissipation in boundary layers 3.
∆T∆T
Since from (2.3.1) F = K 2δth , Fbasic = K d ,Nu = F /Fbasic ∼dd∼ ν Pm1/3 ∼ Re 1/2 Pm1/3 .
2δ
2δ th

(2.3.18)

Now we use the dissipation integral (2.3.13)
Z dh∇u : ∇ui dz =0

ν2
RaNuPr −2 ,d3

Nu >> 1.

Estimating the dissipation in the b.l.s,Z dh∇u : ∇ui dz ∼=0

2
U2
Re 2 ν 2 Re 1/2
5/2 ν
=
=
Re
δν
d2
d
d3

⇒ RaNuPr −2 ∼ Re 5/2 ⇒ Nu ∼ Ra1/4 Pr −1/12 .
(2.3) Nonlinear Rayleigh-Bénard convection(2.3.19)(2.3.20)36/85

2.4 Rotating flows(2.4) Rotating flows

37/85

Rotating flows

We now consider the onset of thermal convection in a rotating
frame. Experimentally this is done by putting the convecting layer
on a rotating turntable.
Also interested in convection in planets and stars, which usuallyrotate. Here the gravity is directed towards the centre of theplanet or star, which rotates about its polar axis.
We can consider the rotational acceleration as providing a Coriolis
force 2ρΩ × u, u being the fluid velocity, and a centrifugal forceρΩ × (Ω × r). Centrifugal force can be written as a gradient of apotential, and so can be combined with the gravity force. Coriolis
force is the crucial new ingredient.

(2.4) Rotating flows

38/85

The centrifugal terms

The effects of the Coriolis force are often more significant than thecentrifugal force. If the rotation vector is in the z direction, and wetake cylindrical polars (R, φ, z) then Ω × (Ω × r) = Ω2 R R̂. But
this is just the gradient of the potential ∇ 12 Ω2 R 2 .
We already have a term like this: gravity. So the centrifugal term
just changes the gravity force. On Earth, the centrifugal term isquite small compared to the g , so its more convenient to think of
the centrifugal force just changing gravity a bit rather than being afundamentally new term.
In some experiments, the centrifugal term can be larger thangravity, which is a way of mimicking a radial gravity as found in aplanet or star.

(2.4) Rotating flows

39/85

When are Coriolis terms important?
Compare u · ∇u and 2Ω × u. If the flow is varying on a lengthscale L, then u · ∇u ∼ U 2 /L very roughly. 2Ω × u ∼ ΩU again
roughly, so the ratio of these terms is the Rossby number
Ro =

U
.
ΩL

(2.4.1)

If the Rossby number is large, Coriolis force won’t be importantbut if it is small it will be an important part of the acceleration.
For the Earth Ω ∼ 7 × 10−5 , a weather system can easily be 1000km (= 106 m) across, and wind speeds are typically of order 10
m/s, giving Ro ∼ 1/7 quite small, so rotation is important.
Same is true for ocean currents. Flows in rapidly rotatingmachinery can also have Ro small.

(2.4) Rotating flows

40/85

Taylor-Proudman theorem: geostrophic flow
Suppose we have slow, (low Ro), steady (no ∂u/∂t), inviscid flow.
Then 2ρΩ × u = −∇p + ρg, including the centrifugal term ingravity.
The curl of this equation gives∇ × 2ρ(Ω × u) = 2ρΩ(∇ · u) − 2ρ(Ω · ∇)u = 0⇒ (Ω · ∇)u = 0(2.4.2)since Ω(∇ · u) = 0 if the fluid is incompressible. The velocity in arotating fluid is independent of the rotation axis direction, z. The
fluid has to flow in columns, independent of z.
This is the Taylor-Proudman theorem, and is the most importantand remarkable result in rotating fluids. It is why rotating fluids
behave so differently from non-rotating fluids.

(2.4) Rotating flows

41/85

Taylor’s penny experiment

Taylor built an experiment inwhich he towed a penny along
the bottom of a rotating tankof water. The column of water
is made visible by adding dye.
All the water in the column of fluid above the penny moved alongwith the penny.
The fluid near the penny had to move with the penny as it was
towed along, and by the z-independence of the flow, all the fluid in
the column above has to move too.

(2.4) Rotating flows

42/85

Dye in non-rotating fluid movie

The fluid is at rest and the tank is stirred by hand, then dye isadded. Fluid disperses in three dimensions, very different from the
2D rotating case.

(2.4) Rotating flows

43/85

Dye in rotating fluid movie

The fluid is on a rotating turntable, and the camera is fixed in theframe of the turntable. The tank is stirred by hand, then dye isadded. Main view is from above, but side view shown at the start.
Fluid moves in two-dimensional sheets.

(2.4) Rotating flows

44/85

2.5 Plane layer rotating convection(2.5) Plane layer rotating convection

45/85

Rotating Convection problems
The rotating plane layer

Rotation parallel to gravity, centrifugal acceleration negligible. Thisis the most common configuration for doing laboratory
experiments, and the one where most is known about the nonlinear
behaviour.
A standard Rayleigh-Bénard convection cell is put on a turntablerotating fast enough for the Rossby number to be small, but not sofast that centrifugal terms are important.
(2.5) Plane layer rotating convection

46/85

Plane layer experiments

Left: Bénard convection with no rotation. Right: Bénard
convection with mild rotation. H.T.Rossby, JFM 36, 309 (1969).
Note the smaller length scale cells in the rotating case.
(2.5) Plane layer rotating convection

47/85

Plane layer model
Chandrasekhar, 1961. The linearised equation of motion is (2.2.2)
with the Coriolis term,∂u1+ 2Ωẑ × u = − ∇p 0 + g αθẑ + ν∇2 u.
∂t
ρ

(2.5.1)

The temperature equation is∂θ= βuz + κ∇2 θ,∂t(2.5.2)β = ∆T /d, and the Boussinesq continuity equation is∇ · u = 0,

(2.5.3)since the fluid is incompressible. Rayleigh-Bénard + Coriolis term.

(2.5) Plane layer rotating convection

48/85

Vorticity equation

Take the curl to eliminate pressure, forming the vorticity equation.
Recall∇ × (a × b) = (b · ∇)a − (a · ∇)b + a∇ · b − b∇ · a,so ∇ × (ẑ × u) = −∂u/∂z.
So
∂u
∂ζ
− 2Ω
= g α∇ × θẑ + ν∇2 ζ.
(2.5.4)∂t∂z
Note that z-vorticity can be created by fluid stretching in thez-direction. Vertical vorticity is now coupled into the linear systemin the rotating case.

(2.5) Plane layer rotating convection

49/85

Inertial waves 1.
We first retain the ∂/∂t term but ignore the others (assume
temporarily the fluid is homogeneous and inviscid) to get
∂ζ
∂u
− 2Ω
= 0.
∂t
∂z

(2.5.5)∂ζ∂ 2∇ u + 2Ω=0∂t∂z

(2.5.6)The curl of this gives

and eliminating ζ gives the wave equation
2
∂2 2
2∂ u
∇
u
=
−4Ω
.
∂t 2
∂z 2

(2.5.7)

These are inertial waves. Vortex lines behave like stretched strings:
if you bend them they vibrate about the straight position.
(2.5) Plane layer rotating convection

50/85

Inertial waves 2.

The solution isu = u0 exp i(kx x + ky y + kz z − ωt),(2.5.8)and then4Ω2 kz2, k 2 = kx2 + ky2 + kz2 ,(2.5.9)k2
is the inertial wave dispersion relation. Note that columnar motion
has kz small, so ω is small. Also we must have |ω| < 2Ω.
ω2 =If viscosity is restored, waves are damped.
Note that if we let ω/2Ω = cos θ, kz /|k| = cos θ, so waves of given
frequency travel with a definite angle θ to the z axis.

(2.5) Plane layer rotating convection

51/85

Inertial waves in the laboratory

Inertial waves excited by a vibrating
cylinder in a rotating water tank. IfΩ · k = |Ω||k| cos θ, cos θ is uniquely set byω, so waves travel at specific angle.
These waves have not been seen in the
Earth’s core, because we have no means ofdetecting them, but they are seen inrotating MHD experiments.
If Ω nearly perpendicular to k, the waves are much slower. If
k · Ω = 0, no z-dependence, the motion is called geostrophic. If
kz << k⊥ it is called quasi-geostrophic (QG). QG waves are called
Rossby waves.

(2.5) Plane layer rotating convection

52/85

Convection in a plane parallel layer
The z-component of the vorticity equation is
∂ζz
∂uz
− 2Ω
= ν∇2 ζz ,∂t∂z

(2.5.10)and the z-component of the double curl gives∂ 2∂ζz∇ uz + 2Ω= g α∇2H θ + ν∇4 uz ,∂t∂z

(2.5.11)and then (2.5.2), (2.5.10) and (2.5.11) give(∂∂
∂ 2 uz∂− ν∇2 )2 ( − κ∇2 )∇2 uz + 4Ω2 ( − κ∇2 ) 2
∂t
∂t
∂t
∂z
= g αβ∇2H (∂− ν∇2 )uz .
∂t(2.5.12)

where ∇2H = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 .
(2.5) Plane layer rotating convection

53/85

Plane layer model: dispersion relation
For stress-free, constant temperature boundaries, the solution is
uz = A sin πz exp i(kx x + ky y ) exp σt,(2.5.13)which satisfies all the stress-free boundary conditions, and we get(σ + π 2 + a2 )2 (Pr σ + π 2 + a2 )(π 2 + a2 ) + (Pr σ + π 2 + a2 )E −2 4π 2
= a2 Ra(σ + π 2 + a2 )(2.5.14)

where
E=

ν
,
Ωd 2

Ra =

g αβ
,
κν

a2 = kx2 + ky2 ,

Pr =

ν
,
κ

(2.5.15)and E is the Ekman number, Ra the Rayleigh number, and Pr the
Prandtl number.
(2.5) Plane layer rotating convection

54/85

Steady modes 1.
First look for neutral modes for the onset of steady convection,
σ = 0. Then (2.5.14) simplifies to
Rasteady =(π 2 + a2 )3
4π 2
+
.
a2
E 2 a2(2.5.16)

Note that as E → ∞ we get the classic Rayleigh-Bénard dispersion
relation. If E is finite, the onset of convection is delayed byrotation, so rotation is normally a stabilising influence.
Letting x = a2 /π 2 ,dRasteady= 0 whendx

2x 3 + 3x 2 = 1 +

4E −2
.
π4

(2.5.17)

As E → 0 , x 3 ∼ 2E −2 /π 4 ora ∼ (2π 2 )1/6 E −1/3 ,(2.5) Plane layer rotating convection(2.5.18)
55/85

Steady modes 2.

So a is large for rapid rotation which means onset occurs as tallthin columns, width ∼ E 1/3 .
This is because viscosity is essential for convection, and at small E
viscosity can only act over short length scales. So the columns
have to be tall and thin. Note this columnar structure is the best
the fluid can do to obey the Taylor-Proudman theorem.
In the small E limit,
Ra ∼ 3a4 ∼ 3(2π 2 )2/3 E −4/3

(2.5.19)so the critical Rayleigh number gets very large at small E .

(2.5) Plane layer rotating convection

56/85

Oscillatory modes 1.

To see if we can get oscillatory modes we let σ = iω in (2.5.14)
and look at real and imaginary parts. This gives−ω 2 Pr (π 2 +a2 )+(π 2 +a2 )3 (2+Pr )+4PrE −2 π 2 = a2 Ra, (2.5.20)−ω 2 (1 + 2Pr )(π 2 + a2 ) + (π 2 + a2 )3 + 4E −2 π 2 = a2 Ra. (2.5.21)
Could eliminate ω 2 between these to get the critical Ra for
oscillatory modes, Raosc . However, first note that (2.5.20) - Pr ×
(2.5.21) gives2ω 2 Pr 2 (π 2 + a2 ) + 2(π 2 + a2 )3 = a2 Ra(1 − Pr ),(2.5.22)

which shows that oscillatory modes can only occur if Pr < 1.

(2.5) Plane layer rotating convection

57/85

Oscillatory modes 2.

However, if Pr < 1 it is possible to find oscillatory modes.
Indeed, if Pr is sufficiently small the minimum critical Rayleighnumber for oscillatory convection is lower than the minimum
critical for steady convection, so if the Rayleigh number is slowly
increased, oscillatory modes onset first at low Prandtl number.
In the limit E → 0 the critical value of Pr for oscillatory modes toonset first is Pr = 0.6766.
As with non-rotating case, other boundary conditions, e.g. no-slip
boumdaries can be dealt with numerically. Nothing radicallydifferent happens.

(2.5) Plane layer rotating convection

58/85

2.6 Nonlinear rotating convection(2.6) Nonlinear rotating convection

59/85

Plane layer model: nonlinear simulations

Numerical simulations(Stellmach et al. 2014) at
E = 10−7 . These are forno-slip boundaries, and thetemperature perturbationis shown. (a) just above
critical Ra, (b-d) at
increasing Ra. Pr = 7.
In (a) the local Rossby
number Ro = U/`Ω issmall, but inertia becomes
important at higher Ra.

(2.6) Nonlinear rotating convection

60/85

Near the onset of convection
For Ra just above Rac :Horizontal and vertical slices of the axial vorticity (E = 5 × 10−6 ):Convective structures near the onset:• Multitude of small vortices of either sign, driven directly bybuoyancy• Elongated structures with horizontal size decreasing with Ekman
number• Mid-plane antisymmetry of the axial vorticity(2.6) Nonlinear rotating convection

61/85

Formation of large-scale vortices (LSV)For Ra ≈ 3Rac :During the growing phase:

During the saturated phase:

Clustering of small cyclonic vortices into a fast cyclonic circulation.
Convective vortices advected by a slower anticyclonic circulation.
Regions of intense shear: horizontal stretching of the convective
vortices.
(2.6) Nonlinear rotating convection

62/85

Structure of the large-scale vortices (LSV)
Horizontal and vertical slices of the axial vorticity

Flow dominated by a cyclone: never an anticyclone. Always grows
to the box size.
LSVs are mostly z-invariant, but outside the LSV small structuresare z-dependent.
Periodic horizontal boundary conditions: horizontal average of ζz is
zero.
LSVs essentially horizontal motions (not a convective structure)(2.6) Nonlinear rotating convection

63/85

Nusselt number versus Rayleigh number
King et al. 2012.
100

Ekman

Numerical Simulations:, ... Pr = 1, ... Pr = 7, ... Pr = 100∞-3

Nusselt Number

10

40

10

-4

4x10
10

3x10

10

1 3
10

10

-5

-5
-6

Laboratory Experiments:, ... Present Work
Rossby (1969)
,

-6

4

10

5

10

6

7

10
10
Rayleigh Number

8

10

9

10

10

10

Rotation delays the onset of convection, but the Nusselt number
then grows rapidly.
When Ro ∼ O(1) the Nusselt number curve goes back to thenonrotating value. LSVs only found on the rapidly rotating part ofthe curve.
(2.6) Nonlinear rotating convection

64/85

2.7 Rotating spherical conection: Busse annulus(2.7) Rotating spherical conection: Busse annulus

65/85Rotating convection in spheres

Left. Sketch of the nature of convection between concentricrotating spheres. Right. Experimental apparatus. Centrifugal
acceleration much larger than gravity. Inner sphere cooled. F.H.
Busse and C.R. Carrigan, Science, 191, 81-83 (1976).
(2.7) Rotating spherical conection: Busse annulus

66/85Rotating convection experiment in a sphere

Snapshots of their convection experiments in rapidly rotating
spheres. Left: small inner sphere. Right: large inner sphere. Notethe convection columns.
(2.7) Rotating spherical conection: Busse annulus

67/85

The Busse annulus model
Busse annulus captures theessential feature, the sloping
boundaries. We assume theflow is almost 2D, so the horizontal motion is described by astreamfunction ψ, but there isa small slope χ to the boundaries. Gravity radially inward,
temperature T = T0 at y = 0,
T = T0 + ∆T at y = d, so
T = T0 + (y /d)∆T + θ.
Gap-width is d, annulus height is L.
We let ux = −∂ψ/∂y and uy = ∂ψ/∂x. The temperature gradientdriving the convection is in the y direction.
(2.7) Rotating spherical conection: Busse annulus

68/85

Busse annulus 1.
The temperature gradient can now drive z-vorticity, and thez-component of the linearised vorticity equation is now
∂ζz
∂uz
∂θ
− 2Ω
= −g α+ ν∇2 ζz ,
∂t
∂z
∂x

(2.7.1)since gravity is in the y -direction, perpendicular to the rotation
vector. So buoyancy drives z-vorticity directly.
Average over z, column length H,
∂ 1
∂t H

Z H
0

2Ω
gα ∂
ζz dz −[uz (H/2) − uz (−H/2)] = −
H
H ∂x
ν
+ ∇2
H

(2.7) Rotating spherical conection: Busse annulus

Z H
θ dz
0

Z H
ζz dz

(2.7.2)069/85

Busse annulus 2.

Now since the slope at the boundaries is a small angle χ,uz = ±χuy on the endwalls z = ±H/2, and we replace thevorticity and temperature with its z average,
∂ζz
4Ωχ
∂θ
−
uy = −g α+ ν∇2 ζz .
∂tH∂x

(2.7.3)Since we have averaged out the z-dependence, and thez-component of velocity is small compared to the x and ycomponents, we can introduce a streamfunctionu = −∇ × ψẑ

(2.7.4)so then ζz = ∇2H ψ, uy = ∂ψ/∂x.

(2.7) Rotating spherical conection: Busse annulus

70/85

Busse annulus 3.
So we get after non-dimensionalising the linearised equations
∂ζz
∂ψ
∂θ
−β
= −Ra
+ ∇ 2 ζz ,
∂t
∂x
∂x

(2.7.5)∂ψ∂θ=−+ ∇2 θ,∂t∂x(2.7.6)ζz = ∇2H ψ,(2.7.7)Pr

where

g α∆Td 3
4Ωχd 3, β=.
(2.7.8)κν
νH
Note that in the annulus model β is no longer the temperature
gradient! β plays the role of the inverse Ekman number here. It isalso related to the beta-effect in geophysical fluid dynamics.
Ra =(2.7) Rotating spherical conection: Busse annulus

71/85

Busse annulus dispersion relation 1.
We look for solutions
ψ = A sin πy exp i(kx − ωt)(2.7.9)Substituting in to (2.7.5)-(2.7.7) gives−iω(π 2 + k 2 )(π 2 + k 2 − iωPr ) + ikβ(π 2 + k 2 − iωPr )= k 2 Ra − (π 2 + k 2 )2 (π 2 + k 2 − iωPr ).

(2.7.10)The imaginary part of this givesω=βk(π 2 + k 2 )(1 + Pr )(2.7.11)and the real part givesRa =(π 2 + k 2 )3β 2 k 2 Pr 2+.
k2
(π 2 + k 2 )(1 + Pr )2(2.7) Rotating spherical conection: Busse annulus(2.7.12)72/85

Busse annulus dispersion relation 2.

In the annulus model convection is always oscillatory. Waves
described by (2.7.11) are thermal Rossby waves.
If we drop the buoyancy term and the viscosity, (2.7.5) and (2.7.7)
give (2.7.11) with Pr = 0 which is the standard Rossby wave
dispersion relation. The restoring force comes from the sloping
boundaries. If a column of fluid moves towards the axis, it
stretches and the vortex tension brings it back.
Since ω/k > 0 these Rossby waves travel eastward. In geophysicalfluid dynamics, the beta effect comes from a variation in the
Coriolis parameter at different latitudes, and has the opposite sign,
so in the atmosphere Rossby waves travel westward.

(2.7) Rotating spherical conection: Busse annulus

73/85

Busse Annulus: large β asymptotics

If we take the large β limit (low viscosity), and minimise Ra overk, the optimising k is large, so we have tall thin columns again.
In fact the lowest critical Ra occurs whenk=

β 1/3 Pr 1/3,
21/6 (1 + Pr )1/3
Ra =

ω=

β 2/3 21/6,
Pr 1/3 (1 + Pr )2/3

3β 4/3 Pr 4/3.
22/3 (1 + Pr )4/3

(2.7.13)

Note that the dependencies are very similar to those in plane layer
convection with β playing the role of E −1 .

(2.7) Rotating spherical conection: Busse annulus

74/85

Nonlinear equations in the annulus
If we add in the advection terms for ζz and θ we get the nonlinear
Busse annulus equations
∂ζz
∂ψ
∂θ
+ (u · ∇)ζz − β= −Ra+ ∇2 ζz ,
∂t
∂x
∂x


∂ψ
∂θ
+ (u · ∇)θ = −+ ∇2 θ,Pr∂t∂xζz = ∇2H ψ.

(2.7.14)(2.7.15)(2.7.16)

These can be integrated in time fairly easily, because they are only
two-dimensional in x and y .
It is also possible to add in a term proportional to −ζz in (2.7.14)(Rotvig and Jones, 2006) to represent friction acting on the top
and bottom boundaries.
(2.7) Rotating spherical conection: Busse annulus

75/85Zonal jet formation in the annulus

The x direction is eastward, andthe y direction is towards therotation axis.
These simulations have
β = 7 × 105 and Ra = 2.5Racrit .
Multiple jet solutions occur when
there is some bottom friction,
stress-free boundaries give fewer
jets.

(2.7) Rotating spherical conection: Busse annulus

76/85

Quasi-geostrophic approximation

The Busse annulus model assumes the angle of the boundary χ issmall. The vertical velocity is then small, and the 2D equations are
asymptotically valid as χ → 0.
In spherical geometry, the angle χ is not small, but neverthelessflow behaves similarly to that in full spherical simulations.
In the QGA, the 2D annulus equations in s and φ are integrated
forward in time, but the varying slope of the boundary is taken into
account.
The QGA model is popular, because 2D simulations are much lessdemanding than 3D simulations.

(2.7) Rotating spherical conection: Busse annulus

77/85Onset of convection in a sphere

Contours of axial vorticity.
Sequence with
E = 3 × 10−5 ,10−5 ,E = 3 × 10−6 , 10−6 .
Picture by Emmanuel
Dormy.
Note convection onsets
first near tangent cylinder(where field is strongest in
the dynamo). Convection
columns get thinner at
smaller E .

B(vi). Onset of convection in a sphere78/85

Rapidly rotating convection: equatorial sections

Contours of axial vorticity
at E = 10−7 . Top,numerical calculations.
Bottom asymptotic theory.
Left internal heating, rightdifferential heating.
For radius ratio 0.35,
internal heating sets infirst in the interior. Notegood agreement between
asymptotics and numerics.
Note also spiralling natureof solution.

B(vi). Onset of convection in a sphere79/85

2.8 Scaling laws in convection(2.8) Scaling laws in convection

80/85

Scaling laws

No simulations can reach the very small Ekman number found inplanetary cores.
The aim of scaling law work is to understand how key quantities
such as the convective velocity and heat flux scale with Ekman
number. The hope is that we can then use simulations at relatively
large E to estimate what will happens at small E .
Physical arguments are based on the force balance:
Inertia

Coriolis

B(vii). Scaling laws in Rotating Convection

Buoyancy

Lorentz

Viscous

81/85

Heat transport and typical velocity
In nonlinear theories, the temperature fluctuation has to be
estimated from heat flux,
Z
Fconv =
ρcp Ur θ dS/4πr 2 ∼ ρcp U∗ T∗ ,(2.8.1)S

where U∗ and T∗ are root mean square velocity and temperaturefluctuations, and cp is the specific heat.
This assumes that there is a strong correlation between hot fluid
and rising fluid. In rotating convection this is not so clear!
In strongly nonlinear convection, balance is between inertia andbuoyancy, the mixing length theory, soU∗2 /d ∼ g αT∗ ∼ g αF /ρcp U∗ .

(2.8.2)

In compressible convection d is usually taken as the density scaleheight, in Boussinseq convection as the distance between theboundaries.
B(vii). Scaling laws in Rotating Convection

82/85

Mixing length theory

Mixing length theory: U∗2 /d ∼ g αT∗ ∼ g αF /ρcp U∗ .
Gives the Deardorff velocity
U∗ ∼



g αFd
ρcp

1/3
(2.8.3)

which works well in laboratory experiments.
In the core, for 1TW of convective heat flux this gives U∗ about 10
times too big. Suggests that rotation/magnetic field is slowing
down the convection.

B(vii). Scaling laws in Rotating Convection

83/85

Inertial theory of rotating convection
Vorticity equationu · ∇ζ − 2(Ω · ∇)u = ∇ × g αθr̂,(2.8.4)

ingoring viscosity. This gives
ΩU∗
g αθ
U∗2
∼
∼
.
2
d
L⊥
L⊥

(2.8.5)Here U∗ is typical convective velocity. |ω| ∼ U∗ /L⊥ ,d = rcmb − ricb ,L⊥ is length scale perpendicular to z, the roll axis.

L⊥ ∼



U∗ d
Ω

1/2

∼



5 × 10−4 × 2 × 106
7 × 10−5

1/2

∼ 4km

(2.8.6)L⊥ is Rhines length, balance of inertia and Coriolis. On longer
length scales, inertia << Coriolis.
B(vii). Scaling laws in Rotating Convection

84/85

Heat flux estimate: Inertial theory of rotating convection
Convective heat flux per square metre Fconv ∼ ρcp U∗ θ.
Eliminate θ to getU∗= Ro ∼Ωd

g αFconvρcp Ω3 d 2

2/5= (RaQ )2/5 .

(2.8.7)For compositional convection, g αFconv is replaced by the buoyancy
flux.
Fitting data from dynamo simulations, CA2006 obtained
0.41
Ro = 0.85RaQ(2.8.8)very close to inertial scaling. Predicted velocity can be comparedwith westward drift velocity in the Earth’s core.

B(vii). Scaling laws in Rotating Convection

85/85