Surface zonal flows induced by thermal convection in a rapidly
rotating thin spherical shell
Youhei SASAKI1 , Shin-ichi Takehiro2 , Ken-suke Nakajima3 , and Yoshi-Yuki Hayashi4121Department of Mathematics, Kyoto University, Kyoto, Japan
Research Institute for Mathematical Sciences, Kyoto University, Kyoto, Japan
3
Faculty of Sciences, Kyushu University, Fukuoka, Japan
4
Center for Planetary Science, Kobe University,JapanIntroduction

Surface flows of Jupiter and Saturn are characterized
by the broad prograde jets around the equator and thenarrow alternating jets in mid- and high-latitudes. It is
not yet clear whether those surface jets are produced byconvective motions in the “deep” region, or are the resultof fluid motions in the “shallow” weather layer.
“Shallow” models, where the thickness of the atmospheric layer is quite smaller than the radius of the planet,
can produce narrow alternating jets in mid- and highlatitudes, while the equatorial jets are not necessarily prograde [1]. On the other hand, “deep” models, where thethickness of the atmospheric layer is comparable to the radius of the planet, can produce equatorial prograde flows
easily, while it seems to be difficult to generate alternating
jets in mid- and high-latitudes [2]. Recently, Heimpel andAurnou [3] consider a thin spherical shell model and showthat the equatorial prograde jets and alternating jets in
mid- and high-latitudes can produce simultaneously when
the Rayleigh number is sufficiently large and convection
becomes active inside the tangent cylinder. However, inthat study, eight-fold symmetry in the longitudinal direction is assumed and fluid motion is calculated only in theone-eighth sector of the whole spherical shell. The artificial
limitation of the computational domain may influence on
the structure of the global flow field. For example, zonal
flows may not develop efficiently due to the insufficient upward cascade of two-dimensional turbulence, or stability ofmean zonal flows may change with the domain size in thelongitudinal direction.
In the present study, we perform numerical simulationsin the whole thin spherical shell domain while coarse spatial resolution and slow rotation rate compared to [3] areused due to the limit of computational resources.
2Formulation and experiments setup

We consider Boussinesq fluid in a spherical shell rotatingwith angular velocity Ω. The non-dimensionalized governing equations consist of equations of continuity, motion, and temperature [4]. The non-dimensional parameters appearing in the governing equations are the Prandtl
number, Pr = ν/κ, the Ekman number, E = ν/(ΩD2 ),and the Rayleigh number, Ra = αgo ∆T D3 /(κν), whereν, D, κ, α, ro , go , and ∆T are kinematic viscosity, the shell
thickness, thermal diffusivity, the outer radius of the shell,
the thermal expansion coefficient, the acceleration of gravity at the outer boundary, and the temperature contrast

between the boundaries, respectively. The spherical shellgeometry is defined by the radius ratio, χ = ri /ro , whereri is the inner radius of the shell. Both inner and outerboundaries are isothermal, impermeable and stress free.
The initial conditions are zero velocity field with 4th-foldsymmetry temperature perturbation used in the dynamo
benchmark[4]. The control parameters used in this study
are summarized in Table 1. The radius ratio χ is slightly
reduced compared to the previous study [3] due to thelimit of computational resources.

3Numerical method

The numerical method used in the present study is a traditional spectral method, where nonlinear terms are evaluated with the transform method [5].
The toroidal and poloidal potentials of velocity are introduced in order to satisfy the equation of continuity.
The velocity potentials and the temperature field are expanded horizontally by the spherical harmonic functions
and radially by the Chebychev polynomials. The nonlinear terms of the governing equations are evaluated in
the physical space and are converted back into the spectral
space. The time integration is performed using the CrankNicolson scheme for the diffusion terms and the secondorder Adams-Bashforth scheme for the other terms.
Our numerical code is parallelized using both MPI andOpenMP. The MPI implementation is quite simple. EachProcessor Element (P.E.) keeps the whole spectrum data
and the grid data divided in the latitudinal direction. The
inverse transformation are performed in a parallel way oneach P.E. After inverse transformation, the spectrum data
are synchronized with all P.E.
In the present study, The total wave number of sphericalharmonics is truncated at 170, and the Chebychev polynomials are calculated up to the 48th degree. We solve
the governing equations in a full-spherical shell without
assuming any symmetry. The numbers of grid points inthe physical space are 512, 256, and 48 in the longitudinal, latitudinal, and radial directions, respectively. In
order to save computational resources, we use hyperdiffusion with the same functional form as the previous studies[3, 6]: ν = ν0 for l ≤ l0 while ν = ν0 [1 + ε(l − l0 )2 ] forl > l0 , where l is total horizontal wave number. In thisstudy, we set l0 = 85 and ε = 10−2 in all calculations. The
model was run for over 4000 planetary rotations (100000
time step) and total CPU time is about 2 weeks.

Table 1. Values of the control parameters, χ, E, Ra, Pr,
and the output parameters Re, Ro in the numerical experiments. The modified Rayleigh number Ra∗ is the ratio ofbuoyancy force to the Coriolis force[2]. The Reynolds andRossby numbers given here are based on the maximum
value of mean zonal flow.

Control parameterχ = ri /roE = ν/ΩD2Pr = ν/κRa = αgo ∆T D3 /κν(Ra∗ = RaE2 /Pr
Output parameter
ReRo

0.75
10−40.1
5 × 105 , 1060.05, 0.1 )

80

60

60

40402020Latitude

(1)

80

0-20-40-40-60-60-80-80-300-200-100

0100

200

0200

400-1500

-750

0750

1500

Fig. 2. Snapshot of azimuthally-averaged zonal velocity
distributions in a meridional cross section. The left andright panels show the cases with Ra∗ = 0.05 and 0.1, respectively. Red and blue represent prograde (eastward)
and retrograde (westward) flows, respectively.

0-20-400-200

Zonal velocity (Ra* = 0.1)

Zonal velocity (Ra*=0.1)

Zonal velocity (Ra*=0.05)

Latitude-400

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

(1)

-500

Zonal velocity (Ra* = 0.05)

300

400

Reynolds number

500
(1)

-15-10-50

51015(×100 1)

Reynolds number

Fig. 1. Snapshot of azimuthally-averaged zonal velocity
profiles at the outer surface. The left and right panels
shows the cases with Ra∗ = 0.05, and 0.1, respectively.
The dashed lines show the latitude at which the tangent
cylinder intersects the outer surface.

of the National Astronomical Observatory of Japan.
The library for spectral transform ’ISPACK’(http://www.gfd-dennou.org/library/ispack/)andits Fortran90 wrapper library ’SPMODEL library’(http://www.gfd-dennou.org/library/spmodel/) wereused for the numerical calculations.
The products of the Dennou-Ruby project(http://www.gfd-dennou.org/arch/ruby/) were used todraw the figures.
The numerical codes used in the present study arealso available from the ’SPMODEL library’ web site:(http://www.gfd-dennou.org/library/spmodel/).
References4Results

Table.1 summarizes the output results of our numerical simulations. Fig.1 shows the azimuthally-averaged
zonal velocity profiles at the outer surface for two different values of the Rayleigh number. Alternating zonal jets
emerge in mid- and high-latitudes when Ra∗ = 0.1 while
they cannot be observed in the case of Ra∗ = 0.05. Fig.2
shows azimuthally-averaged zonal velocity distributions ina meridional cross section. It can be seen that zonal flows
are almost uniform along the rotation axis. Convective
motions can be observed in the tangent cylinder whenRa∗ = 0.1 while they are not active when Ra∗ = 0.05.

5Concluding remark

Even though the spatial resolution and the rotation rate
are reduced and radius ratio is slightly increased in order tosolve the whole domain of the shell, we have confirmed thatthe equatorial prograde flows and alternating jets in midand high-latitudes are produced simultaneously in thermal
convection in a rapidly rotating thin spherical shell. However, the results of forced two-dimensional turbulence on arotating sphere [7] suggest that those zonal flow structure
might not be maintained when time integration is furthercontinued, which should be resolved in the future.
Acknowledgments
Numerical computations were carried out on the Cray
XT4 at the Center for Computational Astrophysics, CfCA,[1] Schneider, T., and Liu, J., Formation of Jets and Equatorial Superrotation on Jupiter. Journal Atmospheric Science, 66, 579–601, 2009.
[2] Christensen, U.R., Zonal flow driven by strongly supercritical convection in rotating spherical shells, Journal ofFluid Mechanics, 470, 115-133, 2002.
[3] Heimpel, M., Aurnou, J., Turbulent convection in rapidly
rotating spherical shells: A model for equatorial and highlatitude jets on Jupiter and Saturn, Icarus, 187, 540–557,
2007[4] Christensen, U., Aubert, J., Cardin, P., Dormy, E., Gibbons, S., Glatzmaier, G., Grote, E., Honkura, Y., Jones,C., Kono, M., Matsushima, M., Sakuraba, A., Takahashi,F., Tilgner, A., Wicht, J., Zhang, K., A numerical dynamo
benchmark. Physics of the Earth and Planetary Interiors
128, 25–34, 2001
[5] Glatzmaier, G.A., Numerical simulations of stellar convective dynamos. I - The model and method. Journal ofComputational Physics 55, 461–484, 1984.
[6] Kuang, W., Bloxham, J., Numerical Modeling of Magnetohydrodynamic Convection in a Rapidly Rotating Spherical
Shell: Weak and Strong Field Dynamo Action, Journal ofComputational Physics, 153, 51–81, 1999[7] Obuse, K., Takehiro, S., Yamada, M., Long-time asymptotic states of forced two-dimensional barotropic incompressible flows on a rotating sphere. Physics of Fluids, 22,056601, 2010