Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://num-meth.srcc.msu.ru/english/zhurnal/tom_2004/ps/v5r104.ps Äàòà èçìåíåíèÿ: Thu Jan 8 14:06:24 2004 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 22:59:40 2012 Êîäèðîâêà:

Numerical Methods and Programming, 2004, Vol. 5 (http://num­meth.srcc.msu.su) 41
UDC 532.517:537.584
THE SUBGRID PROBLEM OF THERMAL CONVECTION
IN THE EARTH'S LIQUID CORE
M. Reshetnyak 1 and B. Steffen 2
The problem of turbulent thermal convection in the Earth's liquid core is considered. Following
assumptions on decreasing of the spatial scales due to rapid rotation, we propose a subgrid model
of eddy diffusivity used in the large­scale model. This approach makes it possible to model realistic
regimes with small Ekman and Rossby numbers (E ¸ 10 \Gamma14 , Ro ¸ 10 \Gamma8 ) and with a sufficiently large
Rayleigh number Ra ¸ 10 12 . An estimate obtained for the averaged kinetic energy is comparable
with observations. The model includes rotation of the solid core due to viscous torque.
1. Introduction. Convection in the Earth's liquid core caused by radioactive heating and compositional
processes [1] is the subject of numerous researches, usually concerned with geomagnetic field generation as
well. The last few decades are demonstrating a fascinating development in this area of research [2]. Based
on the MHD large­scale equations, numerical simulations can reproduce different geomagnetic and geophysical
phenomena: various properties of the geomagnetic field (e.g., its reversals and spectrum), eastward rotation of
the Earth's inner core as well as the realistic ratio of kinetic and magnetic energies [3].
However, a wide range of spatial and time scales makes the direct numerical simulations (DNS) very
cumbersome. The difficulty is caused by small values of transport coefficients: e.g., the kinematic viscosity of
the liquid core is š M = 10 \Gamma6 m 2 s \Gamma1 and the thermal diffusivity Ÿ M = 10 \Gamma5 m 2 s \Gamma1 (here the superscript M
indicates the molecular values). This gives the following estimates of molecular Reynolds and Peclet numbers:
Re M = VwdL
š M ¸ 10 9 and Pe M = VwdL
Ÿ M ¸ 10 8 , where Vwd = 0:2 ffi year \Gamma1 is the west drift velocity and
L = 3 \Delta 10 6 m is the liquid core scale [4], which corresponds to the regime of highly developed turbulence. In
the case of Kolmogorov's turbulence 3D DNS, simulations require ¸ Re 9=4 = 10 20 grid nodes [5]. Attempts
to use the exact values of these parameters on coarse grids lead to numerical instabilities. The first intuitive
models in geodynamo theory that suppress instabilities at small scales (e.g., the model of hyperdiffusivity [3])
gave rise to new questions concerned with interpretation of the results obtained [6]. The more consistent way
is an application of semiempirical turbulence models [7]. Usually, these models are based on assumptions on
the cascade transfer similar to Kolmogorov's, which gives descriptions of the average effect of small­scale field
fluctuations on the large­scale flow in terms of eddy diffusivity. The recent studies on the subgrid [8] and more
complicated models [9, 10] for the thermal convection and dynamo problems in rotating spheres revealed the
principal possibility of describing the small­scale fluctuations in turbulence with desired Reynolds and Peclet
numbers, much like Kolmogorov's model.
These models work up to the regime of moderate rotation speed. The further increase in the Coriolis
force can reduce the total kinetic energy and even suppress convection at all. From linear analysis it follows
that the critical Rayleigh number depends on the Ekman number like Ra cr ¸ E \Gamma1=3 [11]. Even though the
molecular estimate of the Rayleigh number gives huge numbers Ra M ¸ 10 14 [12, 13], this value is only 5 \Delta 10 2
times larger than the critical value Ra cr [2]. Due to rapid rotation of the Earth, the situation in the liquid core
is more complicated and assumptions on similarity of the field spectral characteristics must be checked very
carefully. We show that the direct applications of the traditional turbulence models based on the mixing­length
assumptions lead to results that differ from the observations by orders of magnitude. The cause of such a
disagreement lies in the daily rapid rotation of the Earth, which gives rise to new characteristic spatial scales
in the core [14]. As a result, the energy distribution in the spectrum changes, which leads to difficulties for the
application of Prandtl--Kolmogorov's approach to eddy diffusivity estimation. Convection at these new scales
plays a crucial role in the energy balance of the whole system and changes the estimate of total energy by orders
1 United Institute of the Physics of the Earth, Russian Acad. Sci, 123995, Moscow, Russian Federation;
e­mail: maxim@uipe­ras.scgis.ru; Research Computing Center of Moscow State University, 119992, Moscow,
Russian Federation; e­mail: reshet@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
c
fl Research Computing Center of Moscow State University

42 Numerical Methods and Programming, 2004, Vol. 5 (http://num­meth.srcc.msu.su)
of magnitude. Even a simple account of these effects leads to substantial change in the rate of energy dissipation
and thus to a better agreement of LSS models with observations.
In Section 2 we introduce the large­scale equations of thermal convection and consider Prandtl--Kolmogorov's
assumptions on eddy diffusion. In Section 3 we recall the foundations of convection in a rapidly rotating body
and obtain an estimate for subgrid diffusion. Afterwards, this estimate is used in the large­scale model, Section 4.
Our discussion of results is in Section 5.
2. The large­scale equations. The problem of thermal convection in the Earth's core can be reduced
to that in a spherical shell. Let the surface of a sphere of radius r 0 (in the spherical system of coordinates
(r; `; ')) rotate with an angular
velocity\Omega about the z­axis. This sphere contains a concentric solid inner
sphere of radius r i ; the outer spherical layer (r i ! r ! r 0 ) is filled with an incompressible liquid (V = 0). The
solid inner sphere may rotate freely about the z­axis due to viscous torque. Convection in the outer sphere
is described in the Boussinesq approximation by the Navier--Stokes equation and by the heat­flux equation.
Choosing L = r 0 as unit of length, we can measure the velocity V, the time t, and the pressure p in units of
Ÿ M =L; L 2 =Ÿ M ,
and2\Omega aeŸ M , respectively. Then, the governing equations can be written in the form
Ro M
Ÿ
@V
@t
+ (V \Delta r)V
–
= \Gammarp \Gamma 1 z \Theta V + Ra M T r1 z + E M r\Delta
$
S; (1)
@T
@t
+ (V \Delta r) (T + T 0 ) = r \Delta (rT ); (2)
where 1 z is the unit vector in z­direction,
$
S is the strain tensor, and T is the temperature fluctuations relative
to the imposed profile T 0 = r i =r \Gamma 1
1 \Gamma r i
. The molecular Rossby Ro M , Ekman E M , and Rayleigh Ra M numbers
appear in these equations:
Ro M = Ÿ M
2\Omega L 2 ; E M = š M
2\Omega L 2 ; Ra M = ffg 0 ffi TL
2\Omega Ÿ M :
Here ff is the thermal expansion coefficient, g is the gravity acceleration, and ffi T is a temperature unit (ffiT ¸
10 \Gamma4 K, see [2]). It should be mentioned that the Rayleigh number for nonrotating bodies is usually given in
the form f
Ra M
= ffgffiT L 3 =Ÿ M š M and that Ra M = E M f
Ra M
.
The nondimensionalized momentum equation for the angular velocity ! of the inner sphere (0 ! r ! r i ) is
of the form
Ro M I
@!
@t
= r i E M
I
S
S 'r
fi fi
r=r i
sin # dS; (3)
where I is the moment of inertia of the inner sphere S and S 'r is the strain tensor component in the spherical
system of coordinates [2]. Equations (1), (2), and (3) are accompanied by the nonpenetrating and no­slip
boundary conditions for the velocity V and for zero temperature fluctuations at the shell boundaries.
The system (1), (2), (3) was successfully studied in the regimes of laminar convection with the use of
different numerical approaches [2, 16]. However, these regimes are still very far from the desired estimates
for the Earth's liquid core Ro M = 10 \Gamma8 , E M = 10 \Gamma14 , and Ra M = 10 14 [2]. Attempts to approach to these
parameters using DNS caused numerical instabilities and required the application of turbulence models [8].
However, even the direct usage of the known turbulence models is not trivial.
To support this point, we offer a simple estimate of eddy diffusivity on the basis of the most popular
mixing­length assumption. Following Prandtl--Kolmogorov's hypothesis, the eddy diffusion at scale l can be
estimated as š T = (''l 4 ) 1=3 , where '' = v 3 =l is the energy dissipation rate and v is the velocity at scale l.
Even the upper­bound estimate based on the main scale l = L and the west drift velocity V = 3 \Delta 10 \Gamma3 m s \Gamma1
gives š T = 2 \Delta 10 3 m 2 s \Gamma1 for an Ekman number of order E T = 2 \Delta 10 \Gamma5 . The more realistic estimate with
V l = ffi V ¸ V
`
l
L
' 1=3
and a usual grid scale of l ¸ 3 \Delta 10 \Gamma2 L yields š T = 15 and E T = 10 \Gamma7 . On the other hand,
this estimate of E T would require resolution of about N' ¸ 2 \Delta 10 2 cylindrical columns in the flow [11], which
necessitates the use of the most powerful modern computers. All this means that the above estimate of š T will
not provide the smooth field behavior assumed in Kolmogorov's turbulence when E T ? 1. Thus, the traditional
methods underestimate the eddy diffusion š T .
Such a situation corresponds to the case for which the classical ideas on the direct cascade of energy from
the main scale L to the dissipative scale are violated and additional information on '' at dissipative scale is

Numerical Methods and Programming, 2004, Vol. 5 (http://num­meth.srcc.msu.su) 43
needed. In Section 3 we show that in the case of a rapidly rotating body the energy in the spectrum is shifted
to small scales unresolved by DNS even at the onset of convection. This is the reason why any attempt to
estimate š T in the turbulent regime at scales compared to the grid resolution lead to the not selfconsistent
behavior of the turbulent model.
The way out of such difficulties is to formulate proper assumptions on the spectral properties of the solution
in the range of high wavenumbers.
3. Predictions of asymptotic analysis. The origin of the problem can be seen from the analysis of
linearized system (1), (2) at the onset of convection in the limit of small Rossby and Ekman numbers. As
shown in [11] (see also the recent paper [17]), at the onset of convection the flow structure tends to develop
columns along z­direction such that @=@' ¸ O(E \Gamma1=3 ), @=@s ¸ O(E \Gamma1=6 ), and @=@z ¸ O(1) when E = Ro ! 0.
Linearization of system (1), (2) leads to the balance of Archimedean and viscous terms in the Navier--Stokes
equation: Ra T ¸ E \Gamma1=3 V . The balance of convective and viscous terms in the heat­flux equation gives
V ¸ E \Gamma2=3 T from which follows the estimate of the critical Rayleigh number Ra cr ¸ E \Gamma1=3 . (For convenience
we omitted superscript M.) Thus, at the onset of convection for system (1), (2), the flow is anisotropic with the
smallest scale l E ¸ E 1=3 L defined by the balance of Coriolis and viscous forces. Note that the scale l E ¸ 10 \Gamma5
is beyond the level of DNS. If this asymptotic behavior is valid, the critical Rayleigh number in the Earth's core
is Ra cr ¸ 10 5 [2]. As we show below, the predicted column­like form of flow is very important for estimates of
subgrid dissipation in the liquid core.
The main assumption is that even in the turbulent regime believed to be in the Earth's liquid core, the
flow tends to elongated structures with the smallest scale l E predicted by linear analysis. It is from the scale
l E , the ideas of the direct cascade of energy are applicable. To simplify the problem, we estimate the isotropic
eddy diffusion on the basis of the scale (l E ). In particular, instead of the estimate of velocity gradient at the
subgrid scale l: V 0 ¸ ffi V =l, we use V 0 ¸ E \Gamma1=3 ffi V , where ffi V ¸ 0:3V is the average variation of velocity at
scale l. In this case the estimate of eddy diffusion gives š T ¸ l 2 V 0 ú 5 \Delta 10 4 m 2 s \Gamma1 and E T ¸ 5 \Delta 10 \Gamma4 . This
estimate of turbulent Ekman number corresponds to N' ¸ 10 columns that can be resolved in the large­scale
models with desired accuracy. To demonstrate these arguments, we propose simulations of system (1), (2), (3)
with the given eddy diffusion š T estimated above.
4. A turbulent model. Numerical results. Equations (1), (2), and (3) are solved using the control­
volume method (Simple algorithm) [18] on the staggered grid (n r ; n ` ; n' ) = (45; 45; 64). This method is
based on a finite­difference approximation and demonstrates very high numerical stability for the regimes
with intensive convection 3 . For ease of calculation, we renormalize these equations, using turbulent diffusion
units, so that instead of Ÿ M the value bŸ = 1 m 2 s \Gamma1 is used. Then, the dimensionless parameters become
Ro T = bŸ
2\Omega L 2 = 4 \Delta 10 \Gamma2 and E T = š T
2\Omega L 2 = 10 \Gamma3 . We consider the three regimes with turbulent Rayleigh
numbers Ra T = ffg 0 ffi TL
2\Omega bŸ = 10 6 , 10 7 , and 10 8 (see the time evolution of kinetic energy EK in Figure 1). The
corresponding Reynolds numbers averaged over the shell volume Re M = b Ÿ
š M
p
2EK are 3 \Delta 10 9 , 6 \Delta 10 9 , and 2 \Delta 10 10
(c.f. with the molecular Reynolds number for the Earth's core based on the west drift velocity Re Earth ¸ 10 9 ).
Some characteristic snapshots of large­scale velocity r­, `­, and '­components are presented in Figure 2.
The observed curls in r­ and '­projections correspond to the columns parallel to the z­axis. These columns may
drift in '­direction. In its turn, the nonzero viscous gradient š @
@r
`
V'
r
'
r=r i
causes rotation of the inner core
(see the evolution of the inner core angular velocity ! in Figure 1). Here the positive value of ! corresponds to
the eastward direction known to occur in the Earth [20]. We emphasize that these maps are a result of averaging
of small­scale (l E ¸ O(E M 1=3
) = 10 \Gamma5 ) structures. So far, the micro­scale Reynolds number r e at scale l E is
still larger than unity and the inertial spectrum for the scales smaller than l E exists. An estimate of r e = vl
š M
with l = E M 1=3 and v = 0:1V gives r e ¸ 10 3 . This spectrum has two parts with the transition point defined
by the balance of inertial and Coriolis terms:
l\Omega ¸ Ro v. The turbulence in the range l E 6
l\Omega is influenced by
rotation; the kinetic energy spectrum is E l ¸ l 2 [21]. For the scales smaller than
l\Omega up to the dissipative scale
l d = Re \Gamma3=4 , Kolmogorov's spectrum E l ¸ l 5=3 reappears.
Summarizing the above results, we conclude that based on the realistic values of Rossby and Rayleigh
numbers and on assumptions on the flow spectrum in the liquid core we obtained a value of kinetic energy EK
comparable to observations. Keeping in mind that the velocity field and the eddy diffusion are associated in
3 See [19] for some special considerations of the control­volume method for the full dynamo problem in a sphere.

44 Numerical Methods and Programming, 2004, Vol. 5 (http://num­meth.srcc.msu.su)
Fig. 1. Evolution of the liquid core angular velocity ! and the kinetic energy EK averaged over the volume for
Ro = 4 \Delta 10 \Gamma7 and E = 10 \Gamma14 ; (1): Ra T = 10 6 --- thick line, (2): Ra T = 10 7 --- thin line,
(3): Ra T = 10 8 --- dashed line
Fig. 2. The snapshots of the velocity field components v r ; v# ; and v ' (from left to right) for the equatorial
sections (the left half­plane): (\Gamma700; 1200), (\Gamma4700; 4100), (\Gamma1200; 1800) and meridional sections for the
axisymmetric parts of the fields (the right half­plane): (\Gamma2400; 1800), (\Gamma3800; 3100), (\Gamma400; 1000). Here the
numbers in parentheses indicate ranges
our model, we consider this agreement to be worth to note.
5. Conclusions. We propose a scenario of turbulent thermal convection in a rapidly rotating body when
the Coriolis force shifts the system to the origin of small scales even at the onset of convection. We also
show that the further increase in intensity of heat sources leads to a turbulent regime which is still far from
Kolmogorov's case. It appears that the linear analysis predictions at the onset of convection are applicable
to our eddy diffusion estimate in the regime of fully developed turbulence. Though the original problem is

Numerical Methods and Programming, 2004, Vol. 5 (http://num­meth.srcc.msu.su) 45
highly anisotropic, the ``isotropic'' estimate of eddy diffusion gives a kinetic energy of the system comparable to
observations. Note that the introduction of the magnetic field will not change the problem in essence, because at
the scales l E considered the corresponding micro­scale magnetic Reynolds number is r m Ü 1 and the magnetic
field decays due to the ohmic dissipation process. On the other hand, it is not yet clear how the west drift
velocity relates to the flow at scales l E , and different interpretations of observations can exist. This question
requires the solution of the full dynamo problem.
Acknowledgements. The first author is grateful to the Central Institute for Applied Mathematics (ZAM)
of Forshungszentrum 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 03--05--64074).
References
1. Braginsky S.I., Roberts P.H. Equations governing convection in the Earth's core and the geodynamo // Geophys.
Astrophys. Fluid Dynamics. 1995. 79. 1--95.
2. Jones C.A. Convection­driven geodynamo models // Phil. Trans. R. Soc. London. 2000. A 358. 873--897.
3. 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.
4. Gubbins D., Roberts P.H. Magnetohydrodynamics of the Earth's core. In: Geomagnetism (ed. by Jacobs). 2. 1--184.
New York: Academic Press, 1988.
5. Frisch U. Turbulence: the legacy of A.N. Kolmogorov. Cambridge: Cambridge University Press. 1995.
6. Zhang K., Jones C.A. The effect of hyperviscosity on geodynamo models // Geophys. Res. Lett. 1997. 24,
2869--2872.
7. Kollman W. Prediction methods for turbulent flows. New York: Hemisphere Publishing Corporation, 1980.
8. Buffett B.A. A comparison of subgrid­scale models for large­eddy simulations of convection in the Earth's core //
Geophys. J. Int. 2003. 153. 753--765.
9. Frick P., Reshetnyak M., Sokoloff D. Combined grid­shell approach for convection in a rotating spherical layer //
Europhys. Lett. 2002. 59. 212--217.
10. Frick P.G., Reshetnyak M.Yu., Sokoloff D.D. Cascade models of turbulence for the Earth's liquid core // Doklady
Earth Sciences. 2002. 387. 988--991.
11. Roberts P.H. On the thermal instability of a rotating­fluid sphere containing heat sources // Phil. Trans. R. Soc.
1968. A 263. 93--117.
12. Kono M., Roberts P.H. Definition of the Rayleigh number for geodynamo simulation // Phys. Earth Planet. Int.
2001. 128. 13--24.
13. Gubbins D. The Rayleigh number for convection in the Earth's core // Phys. Earth Planet. Int. 2001. 128. 2--12.
14. Braginsky S.I., Meytlis V.P. Local turbulence in the Earth's core // Geophys. Astrophys. Fluid Dynam. 1991.
A 55. 71--87.
15. Landau L.D., Lifshits E.M. Gidrodinamika (Hydrodynamics). Moscow: Nauka, 1988.
16. Christensen U.R., Aubert J., Cardin P., Dormy E., Gibbons S., Glatzmaier G.A., Grote E., Honkura Y., Jones
C., Kono M., Matsushima M., Sakuraba A., Takahashi F., Tilgner A., Wicht J., Zhang K. A numerical dynamo
benchmark // Phys. Earth Planet. Int. 2001. 128. 25--34.
17. Jones C.A., Soward A.M., Mussa A. The onset of thermal convection in a rapidly rotating sphere // J. Fluid. Mech.
2000. 405. 157--179.
18. Patankar S.V. Numerical heat transfer and fluid flow. New York: Taylor & Francis, 1980.
19. Hejda P., Reshetnyak M. Control­volume method for the dynamo problem in the sphere with the free rotating inner
core // Studia Geoph. et. Geod. 2003. 47. 147--159.
20. Song X., Richards P.G. Observational evidence for differential rotation of the Earth's inner core // Nature. 1996.
382. 221--224.
21. Zhou Y. A phenomenological treatment of rotating turbulence // Phys. Fluids. 1995. 7, N 8. 2092--2094.
Received 27 October 2003