Mean zonal flows
induced by Boussinesq thermal convection
in rotating spherical shells

Shin-ichi Takehiro∗1Youhei Sasaki1 , Ken-suke Nakajima2 , Yoshi-Yuki Hayashi3

1 Kyoto Univ. , 2 Kyushu Univ.,3 Kobe Univ.

Outline
Background
Formulation
Linear theory
Weak nonlinear theory : generation of mean zonal flows
Full nonlinear calculations
Our calculations
Summary

Outline
Background
Formulation
Linear theory
Weak nonlinear theory : generation of mean zonal flows
Full nonlinear calculations
Our calculations
SummarySurface flows of gas giant planets• Surface flows of Jupiter andSaturn are characterized by thebroad prograde jets around theequator and the narrow alternating
jets in mid- and high-latitudes.
• It is not yet clear whether thosesurface jets are produced byconvective motions in the “deep”
region, or are the result of fluid
motions in the “shallow” weather
layer.
(Vasavada and Showman, 2005)

“Deep” models and “Shallow” models• “Shallow” models◦ 2D turbulence on a rotating sphere◦ Primitive model• Result: Narrow alternating jets in mid-and high-latitudes.
• Problem: the equatorial jets are notnecessarily prograde• “Deep” models◦ Convection in rotating sphericalshells• Result: Produce equatorial prograde

flows easily• Problem: difficult to generatealternating jets in mid- andhigh-latitudes

Outline
Background
Formulation
Linear theory
Weak nonlinear theory : generation of mean zonal flows
Full nonlinear calculations
Our calculations
SummaryModel setup• Boussinesq fluid in a rotating spherical shell.
◦ scaling: the shell thickness, viscous diffusion time, temperature difference.
∇·u = 0, ∂u
RaE2 rE+ (u · ∇)u − ∇2 u + 2k × u + ∇p =
T,
∂t
Pr ro
∂T
1
+ (u · ∇)T = ∇2 T.
∂tPr• Parameters:ν◦ Plandtl number: Pr =κ
αg ∆T D3
◦ Rayleigh number: Ra = oκνν◦ Ekman number: E =ΩD2r◦ radius ratio: η = iroOutline
Background
Formulation
Linear theory
Weak nonlinear theory : generation of mean zonal flows
Full nonlinear calculations
Our calculations
SummarySlowly rotating cases
Busse (1970a)• Asymptotic expansion with small parameter, small rotation rate
• 0th order◦ convection structure without rotation rate
◦ degeneration of Ynm modes with the same degree n• 1st order◦ retrograde propagation• 2nd order◦ preference of the sectorial mode Ynn◦ banana-type structureRapidly rotating cases• Due to the Taylor Proudman theorem, convective
motion is along the rotation axis (columnar
convection, Busse 1970b).
• The columns are distorted by the curvature of theouter boundary and form spiral structure.
• The tilting causes cylindrically outward Reynolds
stress.
• Convection cells propagate in the prograde
direction(Zhang 1992)

Banana or column
Takehiro (2010)• Slowly rotating cases◦ Banana type structure → retrograde propagation• Rapidly rotating cases◦ Columnar type structure → prograde propagation• Explanation with conservation of potential vorticity

Slowly rotating case
Rapidly rotating case
(Ek = 10−1 , Pr = 1)
(Ek = 10−2 , Pr = 1)
Contour surfaces of absolute value of vorticity

Propagation mechanism
Takehiro (2010)• Importance of meridional structure◦ Slow (banana) : vortex tubes bend along the shell◦ Rapid (column) : vortex tubes extend along the rotation axis, suppressed by
the outer spherical boundary.

Slowly rotating case

Rapidly rotating caseSpiral structure and Rossby waves
Reason why spiral structure appears?
• Dispersion relation of Rossby wavesβk
ω=− 2 2k +l
• Incorrect explanation◦ Tilting of the boundaries larger in the outer region→ Topographic β large→ Phase speed of the waves faster→ Tilting of the phase• Phase speed must be identical in critical modes• Correct explanation◦ an increase of the radial wavenumber in the outer region due to conservationof the phase speed.

Spiral structure and Rossby waves• Variation of β in y direction → variation of
y-wavenumber l (ω, k are conserved).

rl(y) = ±

−k2 −

β(y)k
ω

• Large β(y) → Large l(y)→ strong tilting of the phase in theoutward region → spiral structure
emerges• Direction of the tilting◦ Rossby waves propagate outward (the inner
region is convective unstable)◦ Group velocity Cgy = 2βkl/(k2 + l2 )2 must be
positive → positive l(y)(Takehiro and Irie, 2001)

Spiral structure and Rossby waves :dependency of the Prandtl number
Takehiro (2008)• Small Pr → small viscosity → longer propagation of waves → spiralstructure• Large Pr → large viscosity → waves dissipate quickly, shorter
propagation → columnar structure• Calculation with WKB theoryRψ = ψ0 ei[kx+ l(y)dy−ωt] · e−

Pr = 0.1

Pr = 0.3

R

Pr (k2 +l(y)2 )/Cgy (y)dy

Pr = 0.5

Pr = 1.0

Outline
Background
Formulation
Linear theory
Weak nonlinear theory : generation of mean zonal flows
Full nonlinear calculations
Our calculations
SummaryWeak non-linear analysis• Obtaining linear responses to the non-linear terms into which aresubstituted the critical solution(cf. Takehiro and Hayashi 1999)◦ Original equations :∂x+ Lx + N[x] = 0,∂t
◦ Critical solution (linearized equations) :∂x(1)+ Lx(1) = 0,∂t
◦ Weak non-linear fields :Lx(2) = −N[x(1) ].

Generation of mean zonal flow
∼ weak nonlinear calculations
Takehiro and Hayashi (1999)

Generation of mean zonal flow
∼ weak nonlinear calculations• Large Pr : meridional
circulation⇓
equatorial sub-rotation• Small Pr : Reynolds
stress⇓equatorial
super-rotation

Slowly rotating cases with large Pr

AM transport by mean meridional circulations• Convective heat transport⇒ High temp. in the upper layer around theequator⇒ Latitudinal temperature gradients⇒ Induced mean meridional circulation⇒ Transport stellar AM from the equator⇒ Equatorial sub-rotation√
Ek = 1/ 10, Pr = 102

Rapidly rotating cases with large Pr
Thermal wind balance∂u1+ u · ∇u + 2Ω × u = − ∇p − αT g+ν∇2 u∂tρ
• Taking rotation,−(2Ω · ∇)u = −α∇T × g• Longitudinal components:2Ω∂vφ
αg(r) ∂T
=−
∂z
r ∂θ

• Thermal wind balance :◦ Torque by Coriolis force = Torque by
buoyancy

Ek = 10−3 , Pr = 102

Rapidly rotating cases with small Pr• Curvature of the outer boundary⇒ tilting of convection cells⇒ velocity correlation v0φ v0r⇒ Reynolds stress transports AM outward⇒ equatorial super-rotation (cf. Busse 1983)(Busse 2002)

Ek = 10−3 , Pr = 1

Slowly rotating cases with small Pr• Sectorial type horizontal structure⇒ longitudinal velocity is twisted by Coriolis⇒ velocity correlation v0φ v0θ
⇒ Reynolds stress transports AM equatorward⇒ equatorial super-rotation (cf. Busse 1970, 1973)√
Ek = 1/ 10, Pr = 1

Outline
Background
Formulation
Linear theory
Weak nonlinear theory : generation of mean zonal flows
Full nonlinear calculations
Our calculations
SummaryPioneering calculations

Radial velocity
Mean zonal flows
Surface zonal flows• Sun and Schubert (1995)◦ Numerical calculations with large Ra, low Ek, high resolution◦ Alternating band structure emerges?

Higher resolution, longer time integrations

Temperature

Axial vorticity
Mean zonal flows(equatorial, meridional)• Christensen (2002)◦ Systematic parameter study for high-Rayleigh number convection◦ Alternating band structure do not emerge!
◦ Interpret alternating bands of SS1995 as residual of initial conditions

“Thin” spherical shell model• Heimpel and Aurnou(2007)(hereafter, HA2007)◦ “Thin” spherical shell modelwith large Rayleigh number,small Ekman number.
◦ Prograde jets and alternating
jets in mid- and high-latitudes
can produce simultaneously◦ However, eight-fold symmetry in
the longitudinal direction isassumed.

Outline
Background
Formulation
Linear theory
Weak nonlinear theory : generation of mean zonal flows
Full nonlinear calculations
Our calculations
SummaryPurpose• HA2007: eight-fold symmetry in the longitudinal direction is assumed.
◦ The artificial limitation of the computational domain may influence on the
structure of the global flow field.
• Zonal flows may not develop efficiently due to the insufficient upward cascade of

two-dimensional turbulence• Stability of mean zonal flows may change with the domain size in the longitudinal

direction.

• In the present study:◦ Numerical simulations in the whole thin spherical shell domain.
◦ Coarse spatial resolution and slow rotation rate are used due to the limit ofcomputational resources.

Experimental setup• Boundary condition: Isothermal, Impermeable and Stress free.
• Input parameters:

parameters
Prandtl number: Pr
Radius ratio: η
Ekman number: E
Modified Rayleigh number: Ra∗present study

HA2007

0.1
0.75
10−4
0.05, 0.1

0.1
0.85, 0.9
10−6
0.05

◦ the definition of modified Rayleigh number: Ra∗ =• the ratio of Coriolis term and buoyancy term

RaE2 αg∆T D
=
Pr
Ω2 D

• Output parameters:◦ (local) Reynolds number, Re, is equivalent to the non-dimensional velocity inthe chosen scaling.

◦ (local) Rossby number: Ro = EReNumerical methods• Traditional spectral method.
◦ Toroidal and Poloidal potentials of velocity are introduced.
◦ The total wave number of spherical harmonics is truncated at 170, and the
Chebychev polynomials are calculated up to the 48th degree.
• The numbers of grid points:512, 256, and 48 in the longitudinal, latitudinal, andradial directions, respectively.

• In order to save computational resources, we use hyperdiffusion with thesame functional form as the previous studies

(

ν=

ν0 ,
for l ≤ l0 ,
ν[1 + ε(l − l0 )2 ], for l > l0 .

◦ we choose l0 = 85, ε = 10−2
• The time integration is performed using the Crank-Nicolson scheme forthe diffusion terms and the second-order Adams-Bashforth scheme forthe other terms.

Results: Ra∗ = 0.05
Zonal velocity (Ra* = 0.05)
Zonal velocity (Ra*=0.05)

Axial Vorticity

0

0

-400

-200

0

200

400

Zonal velocity (Ra*=0.05)(1)80• Broad prograde equatorial jet

60• Alternating zonal jets do not
emerge in mid- andhigh-latitudes.

20

Latitude

40

0
-20
-40
-60
-80
-500

-400

-300

-200

-100

0

100

Reynolds number

200

300

400

500
(1)

Results: Ra∗ = 0.1
Zonal velocity (Ra* = 0.1)
Zonal velocity (Ra*=0.1)

Axial Vorticity

0

0

-1500

-750

0

750

1500

Zonal velocity (Ra*=0.1)(1)• Broad prograde equatorial jet

8060• Outer surface: alternating
zonal jets in high latitudes ?

40
20

Latitude• Mean zonal flow: alternating
zonal jets emerge in mid- andhigh-latitudes

0
-20
-40
-60
-80
-15

-10

-5

0

Reynolds number

5

10

15
(×100 1)

Results: output parameters• Input parametersparameters
Prandtl number: Pr
Radius ratio: η
Ekman number: E
Modified Rayleigh number: Ra∗present study

HA2007

0.1
0.75
10−4
0.05, 0.1

0.1
0.85, 0.9
10−6
0.05

• Output parametersparameters
local Reynolds numnber: Re
local Rossby number: Ro◦ small Re, i.e. weak jet◦ Ro, i.e. slow rotation ratepresent study

HA2007

3.59 × 102 , 1.13 × 103
3.59 × 10−2 , 1.13 × 10−1

5 × 104
1.2 × 10−2 , 2.5 × 10−2Comparison with the Rhines scale• Rhines wavenumber and wavelength :ppkβ = β/2U , λβ = 2π 2U/β.
• Inside the tangent cylinder, assuming 2-dimcolumnar motion,

2 1 dh,β=
E h dsh(s) =qro2 − s2 −q

ri2 − s2 .

• Substituting s = ro cos ϕ, we can get jet
latitudinal wavelength ϕβ(Heimpel and Aurnou 2007).

λβ
ϕβ ≡
= 2π
ro sin ϕsUE(χ2 − cos2 ϕ).
ro sin ϕ cos ϕComparison with the Rhines scale
Zonal velocity (Ra*=0.1)
(1)
80
60
40

Latitude

20
0
-20
-40
-60
-80
-15

-10

-5

0

5

10

15
(×100 1)

Reynolds number

Blue line : ϕβ (ϕ), red cross : jet wavelength from the numerical result

Rhines scale can explains the jet wavelength (qualitatively?)Outline
Background
Formulation
Linear theory
Weak nonlinear theory : generation of mean zonal flows
Full nonlinear calculations
Our calculations
SummarySummary and discussion

.

.

.

.

.

.

.
Summary..
• the profile of mean zonal flow:◦ Broad prograde equatorial jet◦ Alternating zonal jets emerge in mid- and high-latitudes• Outer surface: alternating zonal jets in high latitudes ?
◦ thick shell? , small Ra∗ ?, large E ?
. ◦ hyperdiffusivity?
..
.
In the future...
..
• More ’thin’, ’fast rotating’ spherical shell convection◦ η = 0.75 η = 0.8, 0.85, 0.9◦ E = 10−4 → E = 3 × 10−5.• Investigation of generation mechanism..

References
Busse, F. H., 1970a : Differential rotation in stellar convection zones. Astrophys. J., 159, 629–639
Busse, F. H., 1970b : Thermal instabilities in rapidly rotating systems. J. Fluid Mech., 44, 441–460
Busse, F. H., 1973 : Differential rotation in stellar convection zones. II. Astron. Astrophys., 28, 27–37
Busse, F. H., 1983 : A model of mean zonal flows in the major planets. Geophys. Astrophys. Fluid
Dyn., 23, 153–174[5] Busse, F. H., 2002 : Convective flows in rapidly rotating spheres and their dynamo action. Phys.
Fluids, 14, 1301–1314[6] Christensen, U. R., 2002 : Zonal flow driven by strongly supercritical convection in rotating
spherical shells. J. Fluid Mech., 470, 115–133.
[7] Heimpel, M., Aurnou, J., 2007 : Turbulent convection in rapidly rotating spherical shells: A model forequatorial and high latitude jets on Jupiter and Saturn Icarus, 187, 540–557[8] Sun, Z.-P., Schubert, G., 1995 : Numerical simulations of thermal convection in a rotating sphericalshell at high Taylor and Rayleigh numbers., Phys. Fluids, 7, 2686–2699.
[9] Takehiro, S., Hayashi, Y.-Y., 1999 : Mean zonal flows excited by critical thermal convection inrotating spherical shells. Geophys. Astrophys. Fluid Dyn., 90, 43–77[10] Takehiro, S., Irie, M., 2001 : Two-dimensional spiraling convection in a rapidly rotating cylindricalannulus and the dynamics of Rossby waves. NAGARE Multimedia,
[1]
[2]
[3]
[4]

http://www2.nagare.or.jp/mm/2001/takehiro/index.htm[11] Takehiro, S. 2008 : Physical interpretation of spiralling-columnar convection in a rapidly rotating
annulus with radial propagation properties of Rossby waves. J. Fluid Mech., 614, 67–86.
[12] Takehiro, S., 2010 : On the retrograde propagation of critical thermal convection in a slowly
rotating spherical shell. J. Fluid Mech., 659, 505–515[13] Zhang, K., 1992 : Spiralling columnar convection in rapidly rotating spherical fluid shells. J. Fluid
Mech., 236, 535–556

Taylor Proudman theorem• Coriolis term dominates in the momentum eq. ⇒ geostrophic balance∂u1+ u · ∇u + 2Ω × u = − ∇p−αT g + ν∇2 u∂tρ

.

• Taking rotation operation,.
Taylor Proudman theorem.
..
(2Ω · ∇)u = 0.
.
..
.
• Fluid motion is uniform along the rotation axis ⇒ two-dimensional motion(e.g. Taylor columns)

Conservation of potential vorticity• Conservation of potential vorticity ⇔ Local conservation of angularmomentum.

.

..

where ζz is the axial vorticity component, h is the hight of the fluid
column.

.

.

Conservation of potential vorticity
2Ω + ζz= const.
h.

..