Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://num-meth.srcc.msu.ru/english/zhurnal/tom_2006/ps/v7r110.ps Äàòà èçìåíåíèÿ: Thu Mar 23 13:57:15 2006 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 22:59:59 2012 Êîäèðîâêà: IBM-866 Ïîèñêîâûå ñëîâà: m 31

Numerical Methods and Programming, 2006, Vol. 7 (http://numímeth.srcc.msu.su) 85
UDC 532.517:537.584
SHELL MODELS IN RAPIDLY ROTATING DYNAMO SYSTEMS
M. Reshetnyak 1 and B. Steffen 2
A typical feature of convection of rapidly rotating planets is a geostrophic state. This state is charí
acterized by the following two different scales: the large one along the axis of rotation and the small
transverse one. For the liquid cores of planets, their ratio can be some orders of magnitude already
at the onset of convection. This phenomenon complicates simulations of convection and dynamo
processes and requires a special analysis, essentially in the nonlinear regime. We consider the main
features of the spectra for various intensities of heat sources in the existing spherical models of coní
vection and geodynamo and propose a simple shell model which can mimic some properties of the
original dynamo system.
1. Introduction. The last decade saw some fascinating progress in the modeling of the planetary dynamo,
the Earth's dynamo in particular. While some years ago the main question was whether the heating from below
in the simple Boussinesq system of a rotating sphere can sustain the selfgenerating magnetic field, now there are
a variety of models which, at least at largeíscales, provide a reasonable behavior of a magnetic field comparable
with observations [1]. As it usually happens in astrophysical bodies, the convection and the magnetic field
generation take place at a very wide range of scales. So, using commonly accepted estimates of the Reynolds
number (Re ¦ 10 9 ), the Prandtl number (Pr ¦ 10 8 ) for the Earth's liquid core, for the 3D geometry one
immediately obtains an estimate for the number of degrees of freedom of order ¦ Re 9=4 ¦ 10 20 , which for sure
is out of reach of any modern computer for foreseeable time. The situation with direct numerical simulations
(DNS) of a magnetic field is not so dramatic, since Rm ¦ 10 3 , but even this simulation is at the level of the
modern computer's capacity. So, the DNS cannot handle the desired range of parameters and one should switch
to different methods.
Usually, in many physical and technical problems the large eddy simulations (LES) is a good choice. This
approach provides a description of the effective influence of the fields at subgrid scales on the largeíscale fields.
There is a number of papers on this subject, starting from the rigorous renormalization group analysis [2] to the
simple but robust Prandtl and Smagorinsky models [3]. At least in principle, this approach allows one to provide
a dynamic (in time) inhomogeneous (in space) distribution of turbulent coefficients for use in the original partial
differential equations (PDE). By now, however, practically all the papers on the geodynamo are based on the
assumption that the turbulent coefficients are constant parameters in time and space. Their particular values
are defined by the numerical resolution available in the model [4]. This strange feature of geodynamo models
has a deep background. As we know, without any exception all the planets of the solar system which have their
own magnetic field due to the dynamo process belong to the soícalled class of rapidly rotating bodies, where the
ratio of the solar day œ sd (œ sd
¦\Omega \Gamma1 ,
where\Omega is the daily angular velocity of the body) to the typical convective
time œ c ¦ l
v
(v and l are the typical velocity and scale) is extremely small (the Rossby number Ro = œ sd
œ c
œ 1).
From the early theoretical works of Chandrasekhar summarized in [5] it follows that in the presence of
rapid rotation the thermal convection degenerates to a quasiítwoídimensional form with small gradients of the
fields along the axis of rotation and with large gradients in the perpendicular directions. As a result, the kinetic
energy maximum occurs at very small scales. For the full dynamo problem with a magnetic field, at least for
the moderate modified Rayleigh numbers Ra (in units [6, 7] of its critical value for a given Ekman number E ),
the situation is similar: the position of the maximum of the magnetic and kinetic spectra is very close to the
estimate for the leading mode, so that the critical azimuthal wavenumber m cr ¦ E \Gamma1=3 [8, 9]. On the other
hand, such a coincidence cannot be expected for the regimes where the magnetic diffusion is much stronger than
the viscous dissipation and the magnetic spectrum is shorter than the kinetic energy spectrum, as it is believed
to be in the planets.
1 Institute of the Physics of the Earth, Russian Academy of Sciences, 123995 Moscow, Russian Federation;
eímail: reshetnyak@ifz.ru; Research Computing Center of Moscow State University, 119992, Moscow, Russian
Federation; eímail: rm@uipe.srcc.msu.su
2 Central Institute for Applied Mathematics (ZAM) of Forshungszentrum Jèulich, Dí52425, Jèulich, Germany;
eímail: b.steffen@fzíjuelich.de

86 Numerical Methods and Programming, 2006, Vol. 7 (http://numímeth.srcc.msu.su)
Calculations with higher Ra, close to the ones in the Earth's core, show that the spectrum is nearly flat
up to m cr and then breaks. The wellíknown fact that the spectrum of the observed geomagnetic field decays
as ¦ e \Gamma0:1m [10] gives information on the poloidal part of the field at the surface of the core only. At the
same time, the contribution of small scale fields in the bulk of the core can be already quite considerable, at
least in the directions perpendicular to the rotation axis. The simple estimate of m cr for the Earth's core with
E = 10 \Gamma15 gives m cr ¦ 10 5 . This means that already at the onset of thermal convection the DNS cannot
be used. At the same time, the application of an isotropic and homogeneous semiempirical turbulence model
(like Kolmogorov's one) is not suitable for the description of the geostrophic quasiítwoídimensional turbulence
in the limit of E ! 0. Below we describe the dynamo equations in terms of the shell models. This approach
has been well developed for the analysis of the isotropic turbulence and provided spectral properties of the
stochastic fields for Re AE 1 [11]. We propose its extension to quasiígeostrophic flows. The model developed is a
toy model and can mimic only a few properties of the MHD systems: the energy exchange between the different
components of the fields, the transition of heat energy to the kinetic energy and then to the magnetic energy,
an increase of the critical modified Rayleigh number due to the rapid rotation (¦ E \Gamma1=3 ), and the flat behavior
of the spectrum for m ! m cr , comparable with that in geodynamo DNS.
2. Basic equations and linear analysis.The geodynamo equations for an incompressible fluid (r \Delta V = 0)
in a volume of scale L rotating with an angular
velocity\Omega in the Cartesian coordinate system (x; y; z) in its
traditional dimensionless form can be written as follows:
@B
@t
= r \Theta (V \Theta B) + q \Gamma1 \DeltaB;
E Pr \Gamma1
ß @V
@t + (V \Delta r)V
Ö
= \GammarP \Gamma 1 z \Theta V + RaT z1 z + (r \Theta B) \Theta B + E \DeltaV ;
@T
@t + (V \Delta r)(T + T 0 ) = \DeltaT :
(1)
Here the velocity V , the magnetic field B, the pressure P , and the typical diffusion time t are measured in
the units of ß
L
,
p
2\Omega ßïae,
aeß 2
L 2 , and L 2
ß
, respectively, where ß is the thermal diffusivity, ae is the density, ï is
permeability, Pr = Ú
ß is the Prandtl number, E = Ú
2\Omega L 2 is the Ekman number, Ú is the kinematic viscosity,
j is the magnetic diffusivity, q = ß
j
is the Roberts number, Ra = ffg 0 ffi TL
2\Omega ß
is the modified Rayleigh number,
ff is the coefficient of volume expansion, ffi T is the unit of temperature (for more details see [4]), and g 0 is
the gravitational acceleration. In the case of magnetoconvection, the modified Elsasser number \Lambda appears (it
describes the amplitude of the imposed magnetic field B 0 ): \Lambda = B 2
0
ïaej\Omega
Pr, where the first multiplier is the
Elsasser number in its traditional form. Note that T 0 describes external heating.
Now we are looking for the linear solution \Theta, V = (U x ; U y ; U z ), B = (B x ; B y ; B z ) with the imposed
magnetic field B 0 = (0; 0; B 0 ), neglecting all small quadratic nonlinear terms. The application of the curl and
double curl operators to the Navier--Stokes equation and the curl operator to the induction equation yields:
@B z
@t
= @U z
@z
+ q \Gamma1 r 2 B z ;
@J z
@t = @W z
@z + q \Gamma1 r 2 J z ;
E Pr \Gamma1 @r 2 U z
@t = Er 4 U z + Ra \Delta 1 T \Gamma @W z
@z + \Lambdar 2 @B z
@z ;
E Pr \Gamma1 @W z
@t
= Er 2 W z + @U z
@z
+ \Lambda @J z
@z
;
@T
@t
+ U z = \DeltaT ;
(2)
where W z = rot z V , \Delta 1 = @
@x 2
+ @
@y 2
, and @T 0
@z
= 1. By definition, W is the vorticity of the velocity
field V and J = r \Theta B is an electric current, rT 0 j @T 0
@z
= 1. The substitution of (B z ; J z ; U z ; W z ; T ) =

Numerical Methods and Programming, 2006, Vol. 7 (http://numímeth.srcc.msu.su) 87
(b z ; j z ; u z ; w z ; `)e flt + i(k x x + k y y + k z z) into (2) leads to
flb z = ik z u z \Gamma q \Gamma1 k 2 b z ; flj z = ik z w z \Gamma q \Gamma1 k 2 j z ; fl` = u z \Gamma k 2 `;
\GammaE Pr \Gamma1 flk 2 u z = E k 4 u z \Gamma Ra `k 2 \Gamma ik z w z \Gamma \Lambdaik z k 2 b z ; E Pr \Gamma1 flw z = \GammaE k 2 w z + ik z u z \Gamma i\Lambdak z j z ;
(3)
where 2 \Gamma1=2 k = k x = k y AE k z ¦ 1. From the continuity condition v i k i = 0 (b i k i = 0) it follows that u x ' \Gammau y
(b x ' \Gammab y ). The relations (3) can be represented in matrix form as follows:
A \Delta (b z ; j z ; u z ; w z ; `) T = 0; (4)
where
A =
0
B B B B B B B @
fl + q \Gamma1 k 2 0 \Gammaik z 0 0
0 fl + q \Gamma1 k 2 0 \Gammaik z 0
i\Lambdak z 0 \GammaE (Pr \Gamma1 fl + k 2 ) ik z k \Gamma2 Ra
0 i\Lambdak z \Gammaik z E (Pr \Gamma1 fl + k 2 ) 0
0 0 \Gamma1 0 fl + k 2
1
C C C C C C C A
: (5)
The main results of a dispersion analysis of (3 -- 5) can be found in [12]. In the purely thermal regime without
rotation, we have Ra cr ¦ O(1) and k cr ¦ O(1). The effect of rotation leads to a decrease of the transverse scales,
to an increase of dissipation, and, as a result, to an increase of the critical Rayleigh number: Ra cr ¦ E \Gamma1=3 and
k cr ¦ E \Gamma1=3 . Depending on a value of the Prandtl number, the first marginal mode can be either noníoscillating
(Pr ? 1), or oscillating (Pr ! 1) with ! cr ¦ E \Gamma2=3 . The magnetic field in the zídirection imposed alone without
rotation leads to the same suppression of convection: Ra cr ¦ \Lambda and k cr ¦ Ha 1=3 , where the Hartman number is
Ha = (\LambdaE \Gamma1 ) 1=2 .
For the combined effects (rotation with an imposed magnetic field) and for the weak magnetic field state,
one has Ra cr ¦ E \Gamma1=3 , when \Lambda ! E 1=3 and k cr ¦ Ra cr ¦ O(1) for \Lambda ¦ O(1), the soícalled strong field regime. A
later decrease of Ra cr is a consequence of MHD instabilities which can occur when transition from the weak to
strong regimes takes place. This corresponds to the case where both effects neutralize each other. The reverse
scenario, when the arising smallíscale turbulence damps the largeíscale magnetic field generation, is known as
the dynamo catastrophe (see a review of the problem in [4]). From the very beginning, it is clear that such
an idealized model can only predict the possibility of instabilities, but cannot be a serious foundation for the
faríreaching conclusions on the strong nonlinear regimes of MHD systems without checking it in DNS.
3. Experience from DNS. We classify results on the Boussinesq dynamo into two classes: the weak
convection regime with Ra ¦ (1; : : : ; 10) Ra cr and the strongly convective regime with Ra ¦ (10 2 ; : : : ; 10 3 ) Ra cr .
The detailed statistics of various features of these regimes are presented in [4]. Further we will focus our attention
on the spectral properties of the fields averaged over the bulk of the liquid core.
The typical behavior of the kinetic energy spectrum for E = 10 \Gamma5 without a magnetic field is presented in
Figure 1 a. These calculations are based on the controlívolume model [13] in a rotating spherical shell for the
problem (1) with a grid of 250 3 points. For low Rayleigh numbers, the maximum corresponds to the critical
wavenumber m cr ¦ 10 1 predicted in linear analysis. A similar picture was observed in DNS by [6] and [7] in full
dynamo models. The magnetic energy spectrum was very similar to the kinetic energy spectrum. So, as in the
previous case, the maxima of these spectra agree with the prediction of linear analysis for the purely thermal
convection. Convection at largeíscales is too weak to start an efficient magnetic field generation at these scales
and both maxima are localized at m cr . This is the reason why the columnar structure of flow is not strongly
disturbed by the magnetic field when the dynamo process is switched on.
An increase of the Rayleigh number destroys the regular columnar structure of flow and leads to a flattening
of the spectrum at m ! m cr (Figure 1 b). This phenomenon is well observed in dynamo simulations discussed
in [14], where the break in both the spectra for E ¦ 10 \Gamma6 corresponds to m cr ¦ 20. The typical form of the
spectrum is presented in Figure 1 b. Below we consider the shell model technique, which helps to proceed to
smaller values of E and to longer spectra. It would not be reasonable to expect that this approach will reproduce
all the details for the weak turbulent flows for moderate Rayleigh numbers; however, it is a helpful tool when
the fields a rescalable at the wider range of scales.
4. The dynamo shell models. The idea of the shell model approach is to mimic the original partial
differential equations by a dynamical system with a system of ordinary differential equations. Usually, such an
idea is implemented for the isotropic homogeneous turbulence, see review in [11]. Following the shell approach,

88 Numerical Methods and Programming, 2006, Vol. 7 (http://numímeth.srcc.msu.su)
1 10 100
m
1x10 í5
1x10 í3
1x10 í1
1x10 1
1x10 3
1x10 5
1x10 7
1
2
3
1 ~m í7
1 10 100
m
1x10 í3
1x10 í1
1x10 1
1x10 3
1x10 5
1x10 7
1
2
3
1
~m í7
a) b)
Fig. 1. The Fourier spectra in the azimuthal direction at the equatorial plane of V r ícomponents (1) ,
V ` ícomponents (2) and V'ícomponents (3) of the velocity field for E = 10 \Gamma5 , Pr = 1 and Ra = 150 (a),
Ra = 600 (b) in the largeíscale spherical model. The V ` ícomponent is antisymmetric relative to the equator
plane and its amplitude is considerably smaller than the other two components
one has to write the equations in the Fourier space, making the two assumptions: only local (in kíspace)
interactions occur (e.g., the mode with wavenumber k is produced by harmonics with wavenumbers k \Gamma 1 and
k + 1) and k n = ffk n\Gamma1 (with k 0 = 1), where ff is a constant, usually ff = 2. This logarithmic distribution of the
wave vectors makes it possible to cover a very large range of scales needed in the simulations. As an example,
we present the most popular GOY shell model, named after the three authors of [15, 16]:
dU n
dt
= k n
i
\Gamma 1
8 U n\Gamma2 U n\Gamma1 \Gamma 1
4 U n\Gamma1 Un+1 + U n+1 Un+2
j
\GammaÚ k 2
n U n (6)
with n = 0; : : : ; nmax and with the cutíoff nmax ¦ Re 3=4 . Since this model was developed under the assumption
of isotropy for the velocity field U , there is no Coriolis force in (6). In spite of the apparently simple form of
(6), this model satisfies the two conservation laws of the original Navier--Stokes equation: the conservation of
the kinetic energy E and of the hydrodynamic helicity H in the limit of vanishing viscosity Ú:
E =
X
n
jUn j 2 ; H =
X
n
(\Gamma1) n kn jUn j 2 : (7)
The integration of (6) in time can be done for the very large Reynolds numbers of astrophysical applications,
say Re ¦ 10 10 -- 10 14 , and gives reasonable scaling laws for turbulent flows [11].
Though the decomposition of the nonlinear terms accepted in the framework of the shell models is not unique
in a strict mathematical sense, this approach was very fruitful for the modeling of the classical Kolmogorov'sílike
turbulence. In this kind of turbulence, energy propagates through the spectrum in a ``continuous'' way: from one
scale to the next closest one. For the Navier--Stokes equation and the heat flux equation in three dimensions,
the energy transfer takes place from the large scales to the small ones (the direct cascade).
The next crucial contribution for our modeling was made by the authors of [17, 18] who solved the Boussií
nesq problem (the Navier--Stokes equation in the form (6) but with an Archimedean force ¦ Ra ` n ) and the
shell equation for the temperature `. The ideas of the truncated Fourier series used in the shell models were sucí
cessfully applied to the dynamo equations (the Navier--Stokes equation and the induction equation combined),
see [19]. The main point of this model is to provide the conservation of the total energy (the sum of the kinetic
and magnetic energies) as well as the conservation of the cross helicity in the system. Using these components,
the combined dynamo shell model for the isotropic homogeneous turbulence was presented in [20]. Note that
the shell models can also be used for the description of subgrid fields and eddy diffusion in the controlívolume
models of the largeíscale thermal convection [21].
However, the isotropic shell models can hardly be used for the modeling of the quasiígeostrophic (or
magnetostrophic) turbulence, which is believed to be in the liquid core [22, 23]. As follows from Section 2,
the rapid rotation (as well as the strong magnetic field) leads to the twoídimensionalization of the flow, so

Numerical Methods and Programming, 2006, Vol. 7 (http://numímeth.srcc.msu.su) 89
that the two quite different scales can exist even at the threshold of magnetic field generation: the large scale
along the axis of rotation (or in the direction of the imposed magnetic field B z ) and the small scale in the
perpendicular plane. There is a possibility of a third scale along the strong azimuthal magnetic field [22] which
we do not consider here. This is our motivation to build a twoídimensional shell model. Such a model was first
given in [24] for the thermal convection equation without a magnetic field. Here we will use it as a part of the
full dynamo model, so we recall some basic points. We consider an extension of the linear system (2) to the
nonlinear regime, using a quasiígeostrophic approximation where it is needed. All the fields are represented by
their twoídimensional Fourier modes, where the second index corresponds to the zídirection along the axis of
rotation and the first one corresponds to the perpendicular xídirection. The corresponding wavenumbers are k z
and k x . We start from the linear set of equations (3), rewrite it for the velocity and magnetic fields (instead
of the vorticity w and the electric current j), and add nonlinear terms. The full system of equations for the
magnetic field b (a = b x ; b = b z ), the velocity field u (u = v x ; w = v z ), and the temperature ` are proposed in
the form
da ij
dt
= i''k xj
`
u \Lambda
i+1;j a \Lambda
i+2;j + u \Lambda
i\Gamma1;j a \Lambda
i+1;j \Gamma
u \Lambda
i\Gamma2;j a \Lambda
i\Gamma1;j
2 + u \Lambda
i+2 j a \Lambda
i+1 j \Gamma
u \Lambda
i+1;j a \Lambda
i\Gamma1;j
2 \Gamma
u \Lambda
i\Gamma1;j a \Lambda
i\Gamma2;j
4
'
+ ik zj
`
w \Lambda
i;j+1 a \Lambda
i;j+2 + w \Lambda
i;j \Gamma1 a \Lambda
i;j+1 \Gamma
w \Lambda
i;j \Gamma2 a \Lambda
i;j \Gamma1
2 + w \Lambda
i;j+2 a \Lambda
i;j+1 \Gamma
w \Lambda
i;j+1 a \Lambda
i;j \Gamma1
2
\Gamma
w \Lambda
i;j \Gamma1 a \Lambda
i;j \Gamma2
4
'
\Gamma ik zj u \Lambda
ij (b \Lambda
ij + b \Lambda
i;j \Gamma1 + b \Lambda
i;j \Gamma2 + b \Lambda
i;j+1 + b \Lambda
i;j+2 ) \Gamma q \Gamma1 (k 2
xi + k 2
zj )a ij ;
db ij
dt = i''k xj
`
u \Lambda
i+1;j b \Lambda
i+2;j + u \Lambda
i\Gamma1;j b \Lambda
i+1;j \Gamma
u \Lambda
i\Gamma2;j b \Lambda
i\Gamma1;j
2 + u \Lambda
i+2;j b \Lambda
i+1;j \Gamma
u \Lambda
i+1;j b \Lambda
i\Gamma1;j
2 \Gamma
u \Lambda
i\Gamma1;j b \Lambda
i\Gamma2;j
4
'
+ ik zj
`
w \Lambda
i;j+1 b \Lambda
i;j+2 + w \Lambda
i;j \Gamma1 b \Lambda
i;j+1 \Gamma
w \Lambda
i;j \Gamma2 b \Lambda
i;j \Gamma1
2 + w \Lambda
i;j+2 b \Lambda
i;j+1 \Gamma
w \Lambda
i;j+1 b \Lambda
i;j \Gamma1
2 \Gamma
w \Lambda
i;j \Gamma1 b \Lambda
i;j \Gamma2
4
'
\Gamma i''k xi w \Lambda
ij (a \Lambda
ij + a \Lambda
i\Gamma1;j + a \Lambda
i\Gamma2;j + a \Lambda
i+1;j + a \Lambda
i+2;j ) \Gamma q \Gamma1 (k 2
xi + k 2
zj )b ij ;
d` ij
dt
= i''k xj
`
u \Lambda
i+1;j ` \Lambda
i+2;j + u \Lambda
i\Gamma1;j ` \Lambda
i+1;j \Gamma
u \Lambda
i\Gamma2;j ` \Lambda
i\Gamma1;j
2 + u \Lambda
i+2;j ` \Lambda
i+1;j \Gamma
u \Lambda
i+1;j ` \Lambda
i\Gamma1;j
2 \Gamma
u \Lambda
i\Gamma1;j ` \Lambda
i\Gamma2;j
4
'
+ ik zj
`
w \Lambda
i;j+1 ` \Lambda
i;j+2 +w \Lambda
i;j \Gamma1 ` \Lambda
i;j+1 \Gamma
w \Lambda
i;j \Gamma2 ` \Lambda
i;j \Gamma1
2 + w \Lambda
i;j+2 ` \Lambda
i;j+1 \Gamma
w \Lambda
i;j+1 ` \Lambda
i;j \Gamma1
2
\Gamma
w \Lambda
i;j \Gamma1 ` \Lambda
i;j \Gamma2
4
'
\Gamma (k 2
xi + k 2
zj )` ij + w ij ;
2 Pr \Gamma1 E du ij
dt
= 2 Pr \Gamma1 E i
''
''k xi
`
u \Lambda
i+1;j u \Lambda
i+2;j \Gamma
u \Lambda
i\Gamma1;j u \Lambda
i+1;j
4 \Gamma
u \Lambda
i\Gamma2;j u \Lambda
i\Gamma1;j
8
'
+ k zj
` u \Lambda
i;j \Gamma1 w \Lambda
ij
2 + 2u \Lambda
i;j+1 w \Lambda
ij
'
\Gamma ''k xi
` w \Lambda
i\Gamma1;j w \Lambda
ij
2 + 2w \Lambda
i+1;j w \Lambda
ij
' #
\Gamma 2E (k 2
xi + k 2
zj )u ij \Gamma w ij
k zj
k xi
\Gamma ik zj a \Lambda
ij (b \Lambda
ij + b \Lambda
i;j \Gamma2 + b \Lambda
i;j \Gamma1 + b \Lambda
i;j+1 + b \Lambda
i;j+2 );
Pr \Gamma1 E dw ij
dt = Pr \Gamma1 E i
''
k zj
`
w \Lambda
i;j+1 w \Lambda
i j+2 \Gamma
w \Lambda
i;j \Gamma1 w \Lambda
i;j+1
4 \Gamma
w \Lambda
i;j \Gamma2 w \Lambda
i;j \Gamma1
8
'
+ 2''k xi
` w \Lambda i\Gamma1;j u \Lambda ij
2 + 2w \Lambda
i+1;j u \Lambda
ij
'
\Gamma 2k zj
` u \Lambda i;j \Gamma1 u \Lambda ij
2 + 2u \Lambda
i;j+1 u \Lambda
ij
' #
\Gamma E (k 2
xi + k 2
zj )w i j + u ij
k zj
k xi
+ Ra ` ij \Gamma i''k xi b \Lambda
ij (a \Lambda
i\Gamma2;j + a \Lambda
i\Gamma1;j + a \Lambda
ij + a \Lambda
i+1;j + a \Lambda
i+2;j ):
(8)
The form of the nonlinear terms in the heat flux equation is adopted from the oneídimensional model [18]
and '' is 1 or E 1=3 , depending on the convection regime. This form of the discrete operator provides the
conservation of the heat energy ¦ ` 2 in the limit of vanishing thermal diffusion. The only energy source in this
system is described by the last term in the rightíhand side of the equation.
As regards to the convective and magnetic parts, we stress the following points. In both the equations

90 Numerical Methods and Programming, 2006, Vol. 7 (http://numímeth.srcc.msu.su)
for the velocity and the magnetic field we take into account only one conservation law: the conservation of
energy. Here we do not consider the helicity conservation and the cross helicity conservation (at least, not all
the chains in (8) provide the conservation of these integrals). This makes our model simpler and we hope to
take into account these effects in the future. Note that at least for the three dimensional case, the spectral
characteristics do not change much when one considers either the GOY model with the conservation of energy
and helicity together, or the simpler model [25] with the only one integral, the kinetic energy. Note also that
the helicity conservation h ¦ v z w z in a rapidly rotating body, where both the multipliers are defined at very
different scales, is in principle impossible in terms of the oneídimensional shell model with local interactions in
the kíspace only [26].
To construct shell equations for the velocity and magnetic fields, we ``uncurl'' equations for w z and j z in (3)
and rewrite them for the xícomponent of the fields, adding the nonlinear terms and assuming that w ¦ ik x u
and j ¦ ik x a. For the isotropic convection without the Coriolis force, we have to get two similar equations which
should reproduce properties of the fields similar to the oneídimensional shell model, see (8) with '' = 1.
0
4
8
1
2
1
6
n
1
x
1
0
í
1
4
1
x
1
0
í
8
1
x
1
0
í
2
q
2
/
2
1
x
1
0
í
1
4
1
x
1
0
í
8
1
x
1
0
í
2
1
x
1
0
4
U
2
/
2
1
x
1
0
í
1
4
1
x
1
0
í
8
1
x
1
0
í
2
1
x
1
0
4
W
2
/
2
1
x
1
0
í
1
6
1
x
1
0
í
1
0
1
x
1
0
í
4
1
x
1
0
2
A
2
/
2
1
x
1
0
í
1
7
1
x
1
0
í
1
1
1
x
1
0
í
5
1
x
1
0
1
B
2
/
2
Fig. 2. Spectra of the fields for the
shell dynamo problem (8) without
Coriolis force, E = 10 \Gamma4 ,
Ra = 1:3 \Delta 10 3 , Pr = q = 0:1. Squares
correspond to the zíspectrum and
circles to the xíspectrum
This set of equations (8) provides the conservation of kinetic and
magnetic energies in the limit of zero diffusion, the exchange of kií
netic and magnetic energies, and the feedback of the magnetic field on
the flow. In Figure 2 we present simulations for the regimes without
rotation, '' = 1. The slight asymmetry of the xí and zícomponents
of the spectrum is due to the fact that the Archimedean force in the
system has only the zícomponent. All the spectra have the Kraichnan--
Iroshnikov slope \Gamma3=2 [27, 28]. We call to attention that this regime
mimics the full dynamo model (1) with a prescribed heat source energy
(T 0 ) at n = 0 and with the equipartition state when the kinetic and
magnetic energies V 2
2 and Pr B 2
2E are of the same order of magnitude.
The next step is to take into account rotation. This means that in
the limit of rapid rotation our model should imitate the geostrophic or
magnetostrophic regime. The point is that even for the fully developed
turbulence, when Re ¦ 10 9 , the convective term for the Navier--Stokes
equation in the Earth's core is much less than the Coriolis term. Basing
on the dominant scale of the core and the west drift velocity, one has
that this difference is of three orders of magnitude. In the zero order,
hence, there is a geostrophic balance (the balance of the Coriolis force
and the pressure) or a magnetostrophic balance when the Lorentz
force is taken into account. Note that because Pe
Rm = q \Gamma1 ¦ 10 5 is
quite large and the spectrum of the magnetic field generated should be
shorter than the convective one, the two forms of the balance can exist
at different scales. Then, for a description of turbulence we need to
exclude the large terms from the equations and retain only the terms
which will be of order of convective ones. (Otherwise, the Coriolis
force blocks the energy transfer in the shell model and the spectrum
collapses.) The main features of the geostrophic balance are that w z ¦
2ik x u x and jk x u x j ¦ \Gammajk y u y j AE jk z u z j (similar for j and a). This
means that the convective transport of a quantity in the transverse
direction will be proportional to k z , but not to k x . The physical nature of this phenomenon is trivial. In spite of
the strong gradients of the fields, which would provide an intensive convective flux in the case without rotation,
rotation leads to a swirling of the flow and damps fluxes in the transverse directions. This mechanism is wellí
observed in the different thermal regimes in the Taylor cylinder (where rotation does not prevent a heat flux
in zídirection) and outside the cylinder, where the convective thermal flux is blocked by the Coriolis force [29].
That is why we put '' = E 1=3 when we come to the geostrophic state: jk x uj AE jk z vj. This trick blocks the direct
cascade to the high wavenumbers in xídirection and leads to the flattening of the spectra. We emphasize that
because we started from (3) -- (5), our shell model is able to reproduce features of the system (2) known from
linear analysis (see Section 2).
In Figure 3 we present the dynamo simulations for large and small Roberts numbers q, which can change the
length of the magnetic spectra and, more important, the location of the magnetic spectra relative to the break
of the convective spectra at ¦ m cr . In both the cases, we observe a flat behavior of the fields for wavenumbers

Numerical Methods and Programming, 2006, Vol. 7 (http://numímeth.srcc.msu.su) 91
0
5
1
0
1
5
2
0
2
5
n
1
x
1
í
3
1
1
0
í
1
1
x
1
0
1
1
x
1
0
3
q
2
/
2
1
x
1
0
9
1
x
1
0
1
1
1
x
1
0
1
3
1
x
1
0
1
5
1
x
1
0
1
7
1
x
1
0
1 9
U
2
/
2
1
x
1
0
7
1
x
1
0
9
1
x
1
0
1
1
1
x
1
0
1
3
1
x
1
0
1
5
1
x
1
0
1
7
1
x
1
0
1 9
W
2
/
2
1
x
1
0
0
1
x
1
0
2
1
x
1
0
4
1
x
1
0
6
2
/
2
1
x
1
0
0
1
x
1
0
2
1
x
1
0
4
1
x
1
0
6
2
/
2
0
5
1
0
1
5
2
0
2
5
n
1
x
1
0
í
4
1
x
1
0
í
1
1
x
1
0
2
q
2
/
2
1
x
1
0
7
1
x
1
0
1
1
1
x
1
0
1 5
2
/
2
1
x
1
0
8
1
x
1
0
1
2
1
x
1
0
1 6
2
/
2
1
x
1
0
í
7
1
x
1
0
í
3
1
x
1
0
1
1
x
1
0
5
A
2
/
2
1
x
1
0
í
5
1
x
1
0
í
1
1
x
1
0
3
1
x
1
0
7
B
2
/
2
a) b)
Fig. 3. Spectra of the fields for the dynamo problem (8) with Coriolis force, E = 10 \Gamma10 , Ra = 1:8 \Delta 10 4 ,
Pr = 0:1, q = 10 \Gamma3 (a) and q = 10 2 (b). Squares indicate the zíspectrum and circles the xíspectrum
up to ¦ m cr with a subsequent decrease. Calculations including the Coriolis force show less stability than the
nonrotating simulations, resulting in an irregular form of the spectra. Moreover, the rotation introduces its
own (short) time scale and the approximation of the Coriolis force is local in the wave space. Experiments with
different approximations of the nonlinear terms show that stability gets better when the number of the products
in the chain increases. From this point of view, the effects of the Coriolis force should be very unstable.
The flat behavior of the spectra means that the diffusion in the system is estimated by ¦ m 2
cr (u 2 ; j 2 ). This
is the reason why the eddy viscosity should be increased by a factor of E \Gamma1=3 and can be used in LES [30].
We also note that the flat part of the kinetic energy spectra has faríreaching consequences for the energy
budget of the whole system, since the main contribution to the dissipation in the system is of order ¦ E 1=3 u 2
0 ,
not ¦ Eu 2
0 !
5. Discussion. We have considered a very important property of the rotating turbulence: the flattening
of the spectra for the scales smaller than the diameter of the columns (¦ E 1=3 ). This flattening is caused by
nonlinear interactions. Now, simple linear analysis reveals that the growth rate oe of the eigenmodes with realistic
Ra for the small k x ¦ 10 1 is negative in the Earth's core [31]. Hence, the existing convection in this range of
scales is due only to the nonlinear terms. The DNS of purely thermal convection and of the full dynamo problem
demonstrates that already for the moderate values of Ra the kinetic energy distribution is flat for m ! m cr
and rapidly decays for larger m. So, since in the existing DNS the molecular diffusion coefficients are large, the
magnetic spectrum obtained has a form similar to the kinetic one. To check what will happen for a small q
and a kinetic energy spectrum long enough, one needs methods other than DNS, such that they help to resolve
spectra with different lengths. With this end in view, we used the shell models. Our modeling shows that in this
case the magnetic field is generated mainly at the large scales where the microscale Reynolds number r m = u n
k n j
is sufficient for the magnetic field generation. Usually, the appearance of a magnetic field does not change the
kinetic energy spectrum. This means that our nonlinear system comes to a stable state without the instabilities
predicted in linear analysis. However, a further thorough analysis of the aboveíproposed model is required.
Acknowledgements. The first author is grateful to the Central Institute for Applied Mathematics (ZAM)
of Forschungszentrum in Jèulich for hospitality. This work was supported by the INTAS Foundation (grant # 03--
51--5807) and by the Russian Foundation for Basic Research (grant # 06--05--64619a).

92 Numerical Methods and Programming, 2006, Vol. 7 (http://numímeth.srcc.msu.su)
REFERENCES
1. Kono M., Roberts P. Recent geodynamo simulations and observations of the geomagnetic field // Reviews of Geoí
physics. 2002. 40, N 10. B1--B41.
2. McComb W.D. The physics of fluid turbulence. Oxford: Clarendon Press, 1992.
3. Monin A.S., Yaglom A.M. Statistical fluid dynamics. NY: The MIT Press, 1971.
4. Jones C.A. Convectionídriven geodynamo models // Phil. Trans. R. Soc. London. 2000. A 358. 873--897.
5. Chandrasekhar S. Hydrodynamics and hydromagnetic stability. NY: Dover Publications, 1981.
6. Kutzner C., Cristensen U.R. From stable dipolar towards reversing numerical dynamos // Phys. Earth Planet. Inter.
2002. 131. 29--45.
7. Simitev R.Convection and magnetic field generation in rotating spherical fluid shells (Ph.D. Thesis). Bayreuth:
University of Bayreuth, 2004.
8. Roberts P.H. On the thermal instability of a highly rotating fluid sphere // Astrophys. J. 1965. 141. 240--250.
9. Busse F.H. Thermal instabilities in rapidly rotating systems // J. Fluid Mech. 1970. 44. 441--460.
10.Lowes F.J. Spatial power spectrum of the main geomagnetic field // Geopys. J. R. Astr. Soc. 1974. 36. 717--725.
11.Bohr T., Jensen M., Paladin G., Vulpiani A. Dynamical systems approach to turbulence. Cambridge : Cambridge
University Press, 1998.
12.Hollerbach R., Rèudiger R. The magnetic Universe. Weinheim: WileyíVCH Verlag GmbH & Co. KGaA, 2004.
13.Reshetnyak M., Steffen B. Dynamo model in the spherical shell // Numerical Methods and Programming. 2005. 6,
N 1. 32--38.
14.Roberts P.H., Glatzmaier G.A. A test on the frozeníflux approximation using a new geodynamo model // Phil. Trans.
R. Soc. Lond. 2000. A358. 1109--1121.
15.Gledzer E.B. System of hydrodynamic type admitting two quadratic integrals of motion // Sov. Phys. Dokl. 1973.
18, N 4. 216--225.
16.Ohkitani K., Yamada M. Temporal intermittency in the energy cascade process and local Lyapunov analysis in fully
developed model turbulence // Prog. Theor. Phys. 1989. 81, N 2. 329--334.
17.Mingshun J., Shida L. Scaling behavior of velocity and temperature in a shell model for thermal convective turbuí
lence // Physycal Review E. 1997. 56, N 1. 441--446.
18.Lozhkin S.A., Frick P.G. Inertial Obukhov--Bolgiano interval in shell models of convective turbulence // Fluid Dyí
namics 1988. 33, N 6. 125--140.
19.Frick P., Sokoloff D. Cascade and dynamo action in a shell model of magnetohydrodynamic turbulence // Phys.
Rev. E. 1998. 57. 4155--4164.
20.Frick P., Reshetnyak M., Sokoloff D. Cascade models of turbulence for the Earth's liquid core // Doklady Earth
Sciences. 2002a. 387, N 8. 988--991.
21.Frick P., Reshetnyak M., Sokoloff D. Combined gridíshell approach for convection in a rotating spherical layer //
Europhys. Lett. 2002b. 59, N 2. 212--217.
22.Braginsky S.I., Meytlis V.P. Local turbulence in the Earth's core // Geophys. Astrophys. Fluid Dynamics. 1990. 55.
71--87.
23.Matsushima M., Nakajima T., Roberts P.H. The anisotropy of local turbulence in the Earth#s core // Earth Planets
Space. 1999. 51. 277--286.
24.Reshetnyak M. The estimate of the turbulent viscosity of the Earth's liquid core // Doklady Earth Sciences. 2005a.
400, N 1. 105--109.
25.Desnianskii V.N., Novikov E.A. Simulation of cascade processes in turbulent flows // J. Appl. Math. Mech. 1974.
38. 507--513.
26.Verma M. Statistical theory of magnetohydrodynamic turbulence: recent results // Phys. Reports. 2004. 401. 229--
380.
27.Iroshnikov P.S. Turbulence of a conducting fluid in a strong magnetic field // Sov. Astron. 1964. 7. 566--571.
28.Kraichnan R.H. Inertialírange spectrum of hydromagnetic turbulence // Phys. Fluids. 1965. 8. 1385--1387.
29.Glatzmaier G.A., Roberts P.H. A threeídimensional convective dynamo solution with rotating and finitely conducting
inner core and mantle // Phys. Earth Planet. Inter. 1995. 91. 63--75.
30.Reshetnyak M., Steffen B. The subgrid problem of the thermal convection in the Earth's liquid core // Numerical
Methods and Programming. 2004. 5, N 1. 45--49.
31.Reshetnyak M. Dynamo catastrophe, or why the geomagnetic field is so longílived // Geomagnetism & Aeronomy.
2005b. 45, N 4. 571--575.
Received 11 February 2006