Anisotropic turbulent thermal diffusion and thermal convection in a rapidly rotating fluid sphere

Physics of the Earth and Planetary Interiors 190–191 (2012) 1–9

Contents lists available at SciVerse ScienceDirect

Physics of the Earth and Planetary Interiors journal homepage: www.elsevier.com/locate/pepi

Anisotropic turbulent thermal diffusion and thermal convection in a rapidly rotating fluid sphere D.J. Ivers a,⇑, C.G. Phillips b,⇑⇑ a b

School of Mathematics and Statistics, University of Sydney, NSW 2006, Australia Mathematics Learning Centre, University of Sydney, NSW 2006, Australia

a r t i c l e

i n f o

Article history: Received 28 July 2011 Accepted 21 October 2011 Available online 28 October 2011 Edited by Keke Zhang Keywords: Magnetohydrodynamics Turbulence Anisotropic Dynamo Earth’s core Convection

a b s t r a c t Estimates of the molecular values of magnetic, viscous and thermal diffusion suggest that the state of the Earth’s core is turbulent and that complete numerical simulation of the geodynamo is not realizable at present. Large eddy simulation of the geodynamo with modelling of the sub-grid scale turbulence must be used. Current geodynamo models effectively model the sub-grid scale turbulence with isotropic diffusivities larger than the molecular values appropriate for the core. In the Braginsky and Meytlis (1990) picture of core turbulence the thermal and viscous diffusivities are enhanced up to the molecular magnetic diffusivity in the directions of the rotation axis and mean magnetic field. We neglect the mean magnetic field herein to isolate the effects of anisotropic thermal diffusion, enhanced or diminished along the rotation axis, and explore the instability of a steady conductive basic state with zero mean flow in the Boussinesq approximation. This state is found to be more stable (less stable) as the thermal diffusion parallel to the rotation axis is increased (decreased), if the transverse thermal diffusion is fixed. To examine the effect of simultaneously varying the diffusion along and transverse to the rotation axis, the Frobenius norm is used to control for the total thermal diffusion. When the Frobenius norm of the thermal diffusion tensor is fixed, it is found that increasing the thermal diffusion parallel to the rotation axis is destabilising. This result suggests that, for a fixed total thermal diffusion, geodynamo codes with anisotropic thermal diffusion may operate at lower modified Rayleigh numbers. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction The small estimates for the magnetic, viscous and thermal molecular diffusion suggest that the state of the Earth’s core is turbulent. In particular, the Ekman number E  1014. Thus direct numerical simulation of the geodynamo with Earth-like parameter values is not possible with the present generation of computers. Scales in the magnetic field, the velocity and the temperature less than the presently achievable numerical grid-scale cannot be resolved although they may contain important physics. Geodynamo models (Glatzmaier and Roberts, 1995; Kageyama and Sato, 1995; Kuang and Bloxham, 1997; Christensen et al., 1999; Sakuraba and Kono, 1999) use some form of large eddy simulation with non Earth-like parameter values, typically isotropic diffusivities larger than the molecular values or, in pseudo-spectral codes, hyperdiffusivities, i.e. viscous and thermal diffusivities which increase with wavenumber.

⇑ Principal corresponding author. ⇑⇑ Corresponding author. E-mail addresses: [email protected] sydney.edu.au (C.G. Phillips).

(D.J.

Ivers),

collin.phillips@

0031-9201/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.pepi.2011.10.006

Various investigations analyse and model the sub-grid scale turbulence (Phillips and Ivers, 2000; Matsushima, 2001; Buffett, 2003; Matsui and Buffett, 2005), which include filtering techniques and parameterised models. We have adopted the Braginsky and Meytlis (1990) picture of turbulence in which the turbulence enhances the molecular viscous and thermal diffusivities in preferred directions to values comparable to the molecular magnetic diffusivity (see Phillips and Ivers, 2003). Localized numerical simulations (St Pierre, 1996; Matsushima et al., 1999), which find that the turbulent diffusion is enhanced parallel to the rotation axis and the toroidal magnetic field, give some support to this picture. The present work concentrates on the effects of anisotropic turbulent thermal diffusivity enhanced or diminished parallel to the rotation axis. To separate the effects of the anisotropic thermal and viscous diffusion the viscous diffusivity is kept isotropic. Anisotropic viscous diffusion will be the subject of future work. Anisotropic thermal and viscous diffusion in a rapidly rotating system have been studied previously. Soltis et al. (2004) consider the linear stability of a plane layer rapidly rotating about an axis perpendicular to the layer with vertical gravity, an imposed azimuthal magnetic field about the rotation axis and anisotropic thermal diffusion of the same form as considered herein. The boundaries are perfectly conducting and stress free. Modal

2

D.J. Ivers, C.G. Phillips / Physics of the Earth and Planetary Interiors 190–191 (2012) 1–9

expansions are possible for the solutions, which strongly decouple the magnetic induction, momentum and heat equations to yield a set of ordinary differential equations in the vertical coordinate. Particularly noteworthy is the result of Soltis et al. (2004), that decrease of the thermal diffusion transverse to the rotation axis, with the diffusion parallel to the rotation axis fixed, is destabilizing. Soltis et al. (2006) showed that, for a plane layer rotating about a vertical axis in a uniform vertical magnetic field, the frequencies of axisymmetric instabilities are greater for anisotropic diffusivities than for isotropic diffusivities, if the Elsasser number K [ 30. Moreover, torsional oscillations of decadal periods occur only if the horizontal diffusivities are greater than the vertical diffusivities. Soltis et al. (2006) and Soltis and Brestensky (2010) include the effects of anisotropic viscous diffusion. The present work investigates a spherical model. Geometrically, this resembles the Earth more closely than a plane layer but is mathematically less tractable. In the spherical fluid core V of radius a, the velocity v, magnetic induction B, temperature H in a frame rotating with angular velocity X are governed in the Boussinesq approximation by the equation for momentum, magnetic induction, heat, mass conservation and Gauss’s Law,

  @v þ x  v þ 2X  v ¼ rP þ J  B  qaHg þ qmr2 v ð1:1Þ @t @B ¼ gr2 B þ r  ðv  BÞ ð1:2Þ @t @H Q ; ð1:3Þ þ v  rH ¼ jr2 H þ @t qc p

 ; HÞ ature are modelled by f(r). More general nonlinear models f ðv are possible but are beyond the scope of the present study. We take f(r) = 1 in the present work for simplicity. If there is a stably stratified region at the top of the core (Braginsky, 1993), the model f ðrÞ ¼ r 20  r 2 ; 0 < r 0 < a, is more appropriate and the radial direction is one of the preferred directions. To investigate the effect of anisotropic thermal diffusion on the convection we consider the instability of a steady conductive basic state with zero mean flow in a sphere V of radius a. No-slip, isothermal boundary conditions are imposed on the mean flow and mean temperature,

v ¼ 0;

H ¼ 0;

at r ¼ a:

ð1:9Þ

An alternative boundary condition for the mean temperature is a uniform temperature gradient; see Section 5. In Section 2 the basic state is derived and the problem formulated. The numerical solution method is described in Section 3. In Section 4 results are given for the effect of anisotropic thermal diffusion on spherical thermal convection as the parameter f is varied from 0.95 to 16. The question of whether anisotropy is stabilising or destabilising is considered in Section 5. Concluding remarks are made in Section 6.

q

r  v ¼ 0;

r  B ¼ 0:

ð1:4Þ

where P ¼ p þ 12 qv 2 ; p is a non-hydrostatic reference pressure, x :¼ r  v ; J :¼ r  B=l0 , g is the gravitational acceleration and Q is the heating rate per unit volume. The effect of the centripetal acceleration is neglected. The reference density q, the thermal expansivity a, the viscous diffusivity m, the magnetic diffusivity g, the magnetic permeability l0, the thermal diffusivity j and the specific heat at constant temperature cp are uniform. To further isolate the effects of anisotropic thermal diffusion on the convection, magnetic effects are ignored (B = 0). The velocity, temperature, pressure are decomposed into mean or super-grid ; H; P, which are in scale components indicated by an overline, v principle resolved by numerical computation, and unresolved turbulent or sub-grid-scale components. The equations for the mean  ; H are thus fields v

   @v  v   þ 2X  v  ¼ rP  qaHg þ qmr2 v q þx @t @H Q þ   rH ¼  r  q þv @t qc p r  v ¼ 0:

ð1:5Þ ð1:6Þ ð1:7Þ

The mean anisotropic turbulent thermal diffusion model (Phillips and Ivers, 2000) is

 ¼ jT  rH; q

jT ¼ j½I þ ff ðrÞ1z 1z :

ð1:8Þ

where j is the rank-2 turbulent thermal diffusion tensor and I = 1s1s + 1/1/ + 1z1z is the identity tensor. (In previous work the  and Dj.) Throughout (s,/,z) and (r,h,/) denote authors used q cylindrical and spherical polar coordinates, respectively, with the z-axis along the rotation axis. The only non-zero components of jT are jss = j// = j and jzz = [1 + ff(r)]j. Thermal diffusivity is enhanced (diminished) in the direction of the rotation axis if the anisotropy parameter f > 0 (f < 0), and isotropic in planes transverse to the rotation axis. In the present picture of turbulence f > 0 throughout the core. However, it is possible that f < 0 in stably stratified regions. The turbulent effects of the velocity and temper-

2. Formulation of the problem We consider the instability of a uniformly-rotating conductive  0 ¼ 0; P 0 ; H0 Þ, where P0 is the mean pressure mean basic state ðv and H0 is the steady mean temperature. The associated volume heat production Q 0 is uniform, H0 vanishes on the boundary r = a and g = gar/a, where ga is the gravitational acceleration at the core surface. Eqs. (1.5) and (1.6) for the basic state reduce to

rP0 þ qa H0 g ¼ 0;

r2 H 0 þ f

@ 2 H0 Q0 ¼ : @z2 qjcp

ð2:1Þ

The temperature due to the heat production Q 0 is

1 2

H0 ¼ bða2  r2 Þ;

ð2:2Þ

where the constant

b¼

Q0 : 3qjcp ð1 þ f=3Þ

ð2:3Þ

The temperature gradient is rH0 ¼ br and the radial temperature gradient at the boundary is ba. The pressure is

P0 ¼ 

qag a 2ba

2

H0 ;

discarding an additive constant of integration. Eqs. (1.5) and (1.6), including (1.7) and (1.8), are linearised about the mean basic state (v0 = 0,P0,H0) with Q ¼ Q 0 . To nondimensionalise the problem, length, time, the velocity and the temperature are scaled in terms of the radius a of the sphere V, the viscous diffusion time a2/m,m/a and ba2. Thus the non-dimensional linearised perturbation equations for the mean fields are

 0   0 @v 0  0 þ 1z  v  0 ¼ rP þ R H r E  r2 v @t ! 0 0 @H 1 @2H 2 0 r H þ f 2 ¼ v 0  r  Pr @t @z

ð2:4Þ ð2:5Þ

r  v 0 ¼ 0;

ð2:6Þ 2

3

where E:¼m/2Xa is the Ekman number, R:¼abgaa /2Xm is a modified Rayleigh number, which differs from the usual Rayleigh number due to the rotation, and Pr:¼m/j is the Prandtl number. The non-dimensional turbulent thermal diffusion tensor jT is

3

D.J. Ivers, C.G. Phillips / Physics of the Earth and Planetary Interiors 190–191 (2012) 1–9

jT ¼ Pr1 ½I þ f1z 1z :

ð2:7Þ

From (1.9),

v 0 ¼ 0;

0

H ¼ 0;

at r ¼ 1:

ð2:8Þ

The primes are suppressed below unless otherwise indicated. 3. Numerical solution of the problem The time dependence can be separated out since the basic state fields and boundary conditions are steady,

v ¼ v ðr; h; /Þect ;

P ¼ Pðr; h; /Þect ;

H ¼ Hðr; h; /Þect :

ð3:1Þ

 and H as time-separated The interpretation of the notation v spatial fields will rely on context. The separation constant c is the growth rate of the fields and the problem reduces to a generalised eigenproblem for the eigenvalue c and the associated eigenvector, which are generally complex. The assumption implied by the exponential time-dependence, that the eigenvalues are nondegenerate, holds in the cases considered herein; otherwise polynomial factors in time must be included. The problem can be decoupled into two simpler subproblems  using the symmetry of the mean flow and mean temperature v and H under reflection in the equatorial plane. (The velocity is a polar vector, unlike the magnetic field which is axial.) Since, in this context, an axial dipole [quadrupole] is anti-symmetric [symmetric] in the equatorial plane, anti-symmetry [symmetry] in the equatorial plane will be termed d-symmetry [q-symmetry]. In  is q-symmetric in the terms of spherical polar components, v  r ðr; p  h; /Þ ¼ v  r ðr; h; /Þ, v  h ðr; p  h; /Þ ¼ equatorial plane, if v  h ðr; h; /Þ, v  / ðr; p  h; /Þ ¼ v  / ðr; h; /Þ. The flow is d-symmetric if v v r ðr; p  h; /Þ ¼ v r ðr; h; /Þ, v h ðr; p  h; /Þ ¼ v h ðr; h; /Þ, v / ðr; p  / ðr; h; /Þ. The temperature H is even or e-symmetric h; /Þ ¼ v [odd or o-symmetric], if Hðr; p  h; /Þ ¼ Hðr; h; /Þ½Hðr; p  h; /Þ ¼  0 ; H0 Þ is qe-symmetric the Hðr; h; /Þ. Since the basic state ðv  ; HÞ problem decouples into two subproblems: qe-symmetric ðv  is q-symmetric [d-symmetric], the pressure and do-symmetric. If v P must be even or e-symmetric [odd or o-symmetric]. The time-reduced governing equations are discretised using a Galerkin method in angle and finite-differences in radius to produce an algebraic generalised eigenproblem. The details are given in Appendix A. Table 1 shows the convergence to the most unstable qe-modes for E = 104 at selected values of f. Results for E = 104 are discussed herein. The last two columns of Table 1 show results for the corresponding do-modes, which indicate that do-modes are

more stable than qe-modes, at least for the parameter values shown. This is consistent with the isotropic case (Busse, 1970). An exhaustive search of the dominant do-modes was not conducted. There are no previous studies in the literature of anisotropic diffusion in spherical geometries to provide benchmarks for this study. However, the numerical method can be validated by comparing the results of two sub-problems: the isotropic case f = 0, which corresponds to rotating spherical thermal convection and for which there are various published studies (Roberts, 1968; Zhang, 1992; Ivers and Phillips, 2003; Jones et al., 2000); and the  ¼ 0; f–0, which can be compared stationary anisotropic case v with the analytic solutions for isotropic thermal conduction in spheroids. Isotropic spheroidal thermal conduction problems can be transformed into equivalent spherical thermal conduction problems with an anisotropic diffusion of the form (1.8); see Appendix B.

4. Thermal convection with anisotropic thermal diffusion Fig. 1 shows the critical modified Rayleigh numbers Rc(m,f) and critical angular frequencies xc(m,f) for Pr = 0.1,1,10, E = 104, N = 40, J = 400, qe-mode, and critical or near critical values of the azimuthal wavenumber m. Let Rc ðfÞ be the minimum of Rc(m,f) with respect to m at fixed f and let m⁄(f) be the minimizing value of m. For each Prandtl number, as f increases, Rc ðfÞ increases strictly monotonically and m⁄(f) is monotonic increasing with (unit) jumps. In fact, for all m-modes computed the critical modified Rayleigh numbers Rc(m,f) increase with f. For Pr = 0.1, m⁄ changes from m = 4 to m = 5 near f = 4.6; for Pr = 1, m⁄ changes from m = 5 to m = 6 near f = 0.4; and for Pr = 10, m⁄ changes from m = 7 to m = 8 near f = 2.7 and from m = 8 to m = 9 near f = 11. Further, for all modes, xc(m,f) < 0, i.e. they are prograde or drift eastward. For Pr = 1,10, jxc(m,f)j increases with f for fixed m, but for Pr = 0.1, jxc(m,f)j has a maximum with respect to f for the m = 4, 5, 6 modes and decreases strictly monotonically for m = 7, 8, 9, 10. The jumps in the magnitude of the angular frequency jxc(m⁄,f)j of the m⁄-modes, as f increases, are positive (negative) for Pr = 0.1,1 (Pr = 10). The actual drift rate xc(m⁄,f)/m⁄ > 0 increases with f for Pr = 1,10 and fixed m⁄, but decreases (increases) at jumps in m⁄ for Pr = 0.1,1 (Pr = 10), i.e. the eastward rate of drift of the most unstable m-mode slows down (speeds up) at jumps in m⁄ as f increases. For Pr = 0.1, the drift rate xc(m,f)/m has a maximum for m = 4,5. In the critical modes about 80% of the kinetic

Table 1 Convergence to selected qe-modes for E = 104. The last two columns show the results for the corresponding do-modes. Pr

0.1 0.1 0.1 0.1 0.1 1 1 1 1 1 10 10 10 10 10 10 10

f

0.95 0 4 5 16 0.95 0 0.3 0.4 16 0.95 0 2 3 8 9 16

m

4 4 4 5 5 5 5 5 6 6 7 7 7 8 8 9 9

qe-mode: N = 20, J = 200

qe-mode: N = 40, J = 400

do-mode: N = 20, J = 200

Rc

xc

Rc

xc

Rc

xc

441.886 465.493 568.200 593.133 834.739 105.412 107.844 108.656 108.921 146.684 13.573 13.881 14.587 14.933 16.549 16.851 18.869

466.707 470.118 479.834 513.768 500.036 116.079 118.810 119.684 126.127 155.704 13.549 13.858 14.487 14.383 15.457 14.551 15.357

441.714 465.286 567.787 592.859 834.019 105.494 107.926 108.737 108.997 146.731 13.591 13.899 14.605 14.956 16.578 16.873 18.890

466.884 470.288 479.976 513.897 500.268 116.274 119.004 119.878 126.334 155.941 13.578 13.886 14.516 14.402 15.477 14.619 15.453

1097.506 1161.932 1478.633 1558.647 2483.153 144.240 160.623 164.085 347.024 483.384 33.182 36.919 41.561 44.213 51.720 23.217 27.277

1155.031 1172.549 1207.712 1239.204 1106.288 13.519 15.962 16.350 199.574 318.976 2.499 3.835 5.170 22.540 26.563 7.609 8.857

4

D.J. Ivers, C.G. Phillips / Physics of the Earth and Planetary Interiors 190–191 (2012) 1–9

(a)

(b)

900

−450

(c)

150

−110

20

−12 6

5 4

7 5

800

−470

140

8

−120

4

9 18

6

700

−490

130

−130

5 ωc

R

−510

120

R

6

c

ω

9

c

7 Rc

−14

c

ωc 16

5

600

6

−16

−140

8

7 7 6

500

−530

110

−150

6 14

400

0

5

10

15

−550

100

0

5

ζ

10

15

−160

−18

0

5

10

ζ

15

ζ

Fig. 1. Critical modified Rayleigh numbers (solid curves) and angular frequencies (broken curves) for the qe-mode with (a) Pr = 0.1, (b) Pr = 1, (c) Pr = 10, and E = 104, N = 40, J = 400. The curves are labelled with the azimuthal wavenumber m.

(c)

(b)

(a) 5500

1000 120

5000

900

4500

800

100

4000 9

700 80

3500 600 5

3000 R′ c

500

60

2500 6

400 2000

40

300 1500 200

1000 4

100

500

0

8 20 7

5

0

5

10 ζ

15

0

0

5

10 ζ

15

0

0

5

10

15

ζ

Fig. 2. The adjusted critical modified Rayleigh number R0c for Pr = 0.1, 1, 10, E = 104, N = 40, J = 400, qe-mode. The curves are labelled with the azimuthal wavenumber m.

5

D.J. Ivers, C.G. Phillips / Physics of the Earth and Planetary Interiors 190–191 (2012) 1–9

(a)

(b)

550

−100

500

−150

(c)

120

0

14

0

−20

12

−2

−40

10

−60

8

ωc

110

100 450

−200 90

400

−4

−250 80 −6

350

−300

Rc

ω

70 Rc

−350

60

c

300

50 250

ω

c

R

c

−8

−80

6

−100

4

−400 −10

40 200

−450 30

−12 150

100

−500

0

0.2

0.4 ζ

0.6

−550

0.8

20

10

−120

0

0.2

0.4 ζ

T

0.6

2

0.8

0

0.2

0.4 ζ

T

0.6

0.8

−14

T

Fig. 3. Critical Rayleigh numbers (solid curves) and angular frequencies (broken curves) for the qe-mode with (a) PrT = 0.1, (b) PrT = 1, (c) PrT = 10, and E = 104, N = 40, J = 400. The solid curves correspond to azimuthal wavenumbers: (a) m = 4:7; (b) m = 5:9; (c) m = 7:10; as fT increase near zero.

(a)

(c)

(b)

500

−150

120 −20 110

−200

450

14

−2

12

−4

10

−6

R

c

ω

8

−8

6

−10

4

−12

−30 100 −40

−250

400

90 −50 −300

80 −60

350 ωc

R

c

−350

70 R

c

300

c

−70 ω

c

60

−80

−400 50

−90

250 −450

40

−100

30

−110

20

−120

200 −500

150 0

5

10 ζ

15

−550

0

5

10 ζ

15

2

0

5

10

15

−14

ζ

Fig. 4. Critical Rayleigh numbers (solid curves) and angular frequencies (broken curves) for the qe-mode with (a) PrF = 0.1, (b) PrF = 1, (c) PrF = 10, and E = 104, N = 40, J = 400. The solid curves correspond to azimuthal wavenumbers: (a) m = 4:6; (b) m = 5:9; (c) m = 7:10; as f decreases near 15.

6

D.J. Ivers, C.G. Phillips / Physics of the Earth and Planetary Interiors 190–191 (2012) 1–9

A second hypothesis is that the average thermal diffusivity,

energy is contained in the toroidal flow for 0.95 6 f 6 16. The mean square temperature decreases monotonically with f.

j ¼ 13 tr jT ¼ jð1 þ f=3Þ, together with the other parameters in the first hypothesis except j, are independent of f. The average j , being based on the trace of the thermal diffusion tensor, is inde-

5. Is anisotropic thermal diffusion stabilising or destabilising?

pendent of the orientation of the coordinate system. In this case  mX is explicitly independent of f and is the R :¼ aQ 0 g a a3 =6qcp j Rayleigh number to compare for different values of f. Fig. 1 corresponds to this case. Anisotropy is still stabilizing (de-stabilizing) for f > 0 (f < 0) but the effect is not as strong. For each of the methods above starting with isotropic diffusion and increasing the component of diffusivity in the direction of the rotation axis has the effect of stabilising the system. Moreover, the critical azimuthal wavenumber increases, i.e. the azimuthal lengthscale decreases. A third hypothesis is to start from a system with isotropic turbulent thermal diffusion jT = jTI (in dimensional form), commonly used in geodynamo simulations, and reduce the components of the diffusion tensor in the non-preferred transverse s, / directions toward molecular values. Thus the anisotropic turbulent thermal diffusion tensor may be written as jT = jT[I  fT(1s1s + 1/1/)]. This is equivalent to a re-scaling of (1.8) with jT = j(1 + f) and fT = f(1 + f)1. In dimensionless form

Fig. 1 suggests that increasing f has a stabilizing effect on the fluid, since Rc ðfÞ increases with f. However, such a conclusion, based on the comparison of Rc ðfÞ for different values of f, is deficient since Rc ðfÞ varies explicitly with f. In fact, there is no unique way to compare the stability of systems with different values of f and some hypothesis must be made about the systems. We offer several alternative hypotheses. Our prior hypothesis is to assume that the problem parameters a; Q 0 ; g a ; a; q; cp ; m; X and j are independent of f. In particular, the thermal diffusivity j is regarded as a property of the fluid, not the turbulence, and hence independent of f; the thermal diffusivity tensor jT models the turbulence and depends on f. The relative stability of two systems can be compared using the adjusted critical modified Rayleigh number, R0c :¼ Rc ðfÞð1 þ f=3Þ, at different values of f. Observe that R0c ¼ aQ 0 g a a3 =6qcp jmX is explicitly independent of f and has the same dependence on the problem parameters (apart from f) as the modified Rayleigh number in the isotropic case. Thus anisotropy is de-stabilizing (compared to isotropy) if R0c < Rc ð0Þ and stabilizing (compared to isotropy) if R0c > Rc ð0Þ. Fig. 2 shows plots of R0c versus f for Pr = 0.1,1,10. Under the R0 c-criterion anisotropy is stabilizing (de-stabilizing) for f > 0 (f < 0); the stabilizing effect increases with f. From Fig. 2 the same conclusion holds for a particular m-mode, using a criterion based on the adjusted modified Rayleigh number, Rc(m,f)(1 + f/3), which is also explicitly independent of f.

1

jT ¼ Pr1 T ½I  fT ð1s 1s þ 1/ 1/ Þ;

ð5:1Þ

where PrT:¼m/jT = Pr/(1 + f) is a turbulent Prandtl number. At f = 0, PrT = Pr. In this hypothesis PrT is fixed, f P 0 is varied and the Prandtl number Pr varies according to Pr = (1 + f)PrT. The problem parameters a; Q 0 ; g a ; a; q; cp ; m; X and jT are independent of fT (and f). Results are shown in Fig. 3 for the critical Rayleigh numbers and angular frequencies of the qe-mode with PrT = 0.1, 1, 10,

(a)

(b)

0.5

z

0

−0.5

−1

0

0.2

0.4

0.6

s

0.8

1

0

0.2

0.4

0.6

0.8

1

s

 i for Pr = 01, E = 104, qe-mode and: (a) m⁄ = 5, f = 0.95; (b) m⁄ = 6, f = 16. The meridional and azimuthal components are Fig. 5. The time-averaged convective heat flux hHv represented by the arrows and level contours, respectively.

7

D.J. Ivers, C.G. Phillips / Physics of the Earth and Planetary Interiors 190–191 (2012) 1–9 Table 2 Growth rates of the slowest decaying even and odd modes for m = 1 in the prolate case (1 < f < 0) with a = 1 and in the oblate case (0 < f < 1) with b = 1. f

b(f < 0) or a(f > 0)

0.95 0.6 0.1 0.1 0.5 1 5 11 16

4.472135955 1.581138830 1.054092553 1.048808848 1.224744871 1.414213562 2.449489743 3.464101615 4.123105626

Even

Odd

c (computed)

c (theoretical)

c (computed)

c (theoretical)

15.585430 17.634135 19.784626 18.720482 14.779177 12.046217 6.438837 4.913695 4.419570

15.585443 17.634142 19.784634 18.720489 14.779183 12.046221 6.438839 4.913695 4.419570

17.472041 24.279530 31.786632 31.485970 26.798452 23.489576 16.309292 14.100463 13.330120

17.471992 24.279539 31.786645 31.485983 26.798462 23.489585 16.309296 14.100463 13.330115

E = 104, N = 40, J = 400, as functions of fT. As fT increases Rc decreases quite steeply and the critical azimuthal wavenumber also decreases. Fixing the thermal diffusion along the z-axis and reducing the diffusion in the tranverse s, / directions is destabilising. The scaling of Soltis et al. (2004) and their results are consistent with this case. The above results all suggest the general hypothesis that increasing (decreasing) total thermal diffusion is stabilising (destabilising). To test this suggestion requires a measure of total thermal diffusion, even if anisotropic, such as a norm of jT. Norms usually depend on the q coordinate ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi system which is used. The Frobenius norm, kjT kF :¼ trðjTT jT Þ, where the superscript T denotes the transpose, has the advantage that it is invariant under rotation of the coordinate axes. Using this norm we ffi factor jT = jFj1, where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kj1kF = 1 and jF :¼ kjT kF ¼ j

2 þ ð1 þ fÞ2 . Equating the Frobenius

norms of (2.7) and (5.1) gives, in dimensionless form,

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PrF ¼ Pr= 2 þ ð1 þ fÞ2 ¼ PrT = 1 þ 2ð1  fT Þ2 ;

ð5:2Þ

where PrF:¼jF/m is a new Prandtl number based on jF. Systems with identical jF or PrF, but different f and hence Pr (equivalently fT and PrT) have the same total diffusion as measured by the Frobenius norm. The fourth hypothesis is that the problem parameters a; Q 0 ; g a ; a; q; cp ; m; X and jF are independent of f. Results are shown in Fig. 4 for the critical Rayleigh numbers and angular frequencies of the qe-mode with PrF = 0.1, 1, 10, E = 104, N = 40, J = 400, as functions of f. As f increases Rc decreases quite steeply and the critical azimuthal wavenumber also decreases. Thus fixing the total in the Frobenius norm but reducing the diffusion in the transverse s, / directions is in fact destabilising. It is also possible to compare systems with the alternative boundary condition of a uniform mean-temperature gradient,

@ r H ¼ b;

at r ¼ a;

ð5:3Þ

where b is a prescribed constant, which together with the other problem parameters a,ga,a,q,cp,m,X and j is independent of f. From (2.3) the volume heat production Q 0 now depends on f. The basic state temperature (2.2) solves (2.1)(b) but with the boundary condition (5.3). The perturbation temperature now satisfies the homogeneous boundary condition, 0

@ r H ¼ 0;

at r ¼ 1;

ð5:4Þ

so the computational results given above do not apply. However, for the four hypotheses considered above we expect the same conclusions to hold. In general, (5.3) and (5.4) are not equivalent to a uniform radial heat flux at the boundary: e.g. if the boundary temperature is unir ¼ j½1 þ ff ðaÞ cos2 h@ H=@r, which varform and f(a) – 0, then q ies with co-latitude h. Latitude-dependence of the basic-state mean temperature or temperature gradient on the boundary would produce secondary anisotropic effects on the perturbation

fields through the convective term in the heat Eq. (1.6) and mix the effects of lateral temperature or temperature gradient variation with the anisotropic thermal diffusion. 6. Concluding remarks It has been shown in Section 5 that enhancement of the thermal diffusion in the direction of the rotation axis, with the transverse thermal diffusion fixed, is stabilizing. If the thermal diffusion parallel to the rotation axis is fixed and the transverse thermal diffusion diminished, the effect is to destabilize the fluid. This conclusion might be expected on physical grounds: according to the Proudman–Taylor theorem the rotation constrains the motion to planes transverse to the rotation axis. The rotational constraint on convective heat transfer parallel to the rotation axis is (partially) offset by increased thermal diffusion parallel to the rotation axis. However, if the total thermal diffusion is fixed, as measured by the Frobenius norm, then increasing the thermal diffusion in the direction of the rotation axis, while simultaneously decreasing the transverse thermal diffusion, is destabilizing. This result suggests that, for a fixed total thermal diffusion as measured by the Frobenius norm, geodynamo codes with anisotropic thermal diffusion may operate at lower modified Rayleigh numbers. The effect of the anisotropic thermal diffusion on wavenumbers in our results is small: the azimuthal wavenumber increases (decreases) as f increases (decreases) with the transverse thermal diffusion fixed, and there is little apparent change in the radial lengthscale as f increases. Nonlinear studies are required for weakly and strongly supercritical modified Rayleigh numbers. The time average of the (linearised) heat flux over a period 2p/x vanishes for the critical mode  is nonlinear, and (Rec = 0). However, the convective heat flux Hv its time-average (denoted by angle brackets) in the critical case,  þ ImHImv  Þ, is steady and also independent of  i ¼ 12 ðReHRev hHv  i for the qe-mode with E = 104, f = 0.95, 16, /. Fig. 5 shows hHv Pr = 1. In both cases the meridional flux is strongly radial outward from the z-axis. The major effect of increasing the axial anisotropy of the thermal diffusion, i.e. increasing f, is that the azimuthal com / i becomes essentially eastward. ponent hHv Asymptotic approximation in the rapidly-rotating limit, in particular E ? 0 [see Jones et al. (2000) for the isotropic case], although interesting, is outside the turbulence model considered here. In fact, anisotropic viscous diffusion must be considered. The analogy between the spherical anisotropic turbulence considered here and the spheroidal thermal conduction problem in Appendix B does not carry over to thermal convection in a spheroid (Ivers, 2002) and hence does not allow comparison with the results there, since the equations differ in the pressure gradient. The effect of the magnetic field will be considered in future work.

8

D.J. Ivers, C.G. Phillips / Physics of the Earth and Planetary Interiors 190–191 (2012) 1–9

Appendix A. Details of the numerical method The solenoidal condition (2.6) is satisfied by a toroidal–poloidal  ¼ r  tr þ r  ðr  srÞ. This rerepresentation for the velocity v duces the problem to the four independent scalar fields: the velocity potentials t,s; the modified pressure P; and the temperature H. The no-slip isothermal boundary conditions (2.8) imply

s ¼ 0;

@s=@r ¼ 0;

t ¼ 0;

H ¼ 0;

at r ¼ 1:

ðA1Þ

The pressure is eliminated by using the vorticity equation. The remaining three fields are expanded in mean normalised spherical harmonics,

f ¼

1 X n X

fnm ðrÞY m n ðh; /Þ;

f ¼ s; t; H:

ðA2Þ

n¼0 m¼n

Symmetry in the equator is simply expressed using the parity of the spherical harmonics: the scalar field f ¼ s; t; H is e-symmetric (o is qsymmetric) if fnm ¼ 0 for n  m odd (even); and the velocity v symmetric (d-symmetric) if t is odd (even) and s is even (odd). Two further symmetries can be exploited to reduce the problem size: symmetry under rotation about the z-axis; and the real prop ; H. The qe- and do-symmetric subproblems thus erty of the fields v decouple into the following chains of coefficients for each m P 0,

qe  symmetry :

sm m;

m m m m ; t mþ1 ; smþ2 ;

H

m m mþ2 ; t mþ3 ; . . . ;

H

m m m m m do  symmetry : t m m ; smþ1 ; Hmþ1 ; t mþ2 ; smþ3 ; Hmþ3 ; . . . :

ðA3Þ ðA4Þ

If m ¼ 0; s00 or t 00 are deleted, but not H00 , if f – 0, since it is coupled to H02 in the anisotropic thermal diffusion term. The problem is discretised in angle using a Galerkin method which produces compact spectral equations. All vector fields in Eqs. (2.4) and (2.5) are expanded in vector spherical harmonics Ym n;n1 (see Ivers and Phillips, 2008),

F¼

1 X

nþ1 X

n X

m Fm n;n1 ðrÞY n;n1 ðh; /Þ;

: ; q F ¼ rH0 ; v

ðA5Þ

n¼0 n1 ¼n1 m¼n

Equations are then extracted for the spherical harmonic coefficients of t; s; H. The heat Eq. (2.5) is actually discretised with the more general anisotropic thermal diffusion term, r  ðFG  rHÞ (Phillips and Ivers, 2000). Eq. (1.8) (with f = 1) is the special case F = G = 1z, which greatly limits the coupling since 1z ¼ Y01;0 . The expansions (A2), (A5) are truncated with a triangular truncation level N, i.e. n 6 N, jmj 6 N; expansions (A5) are truncated in n1 with n1 6 N + 1. Values of N are given below. The radial discretisation uses second-order finite-differences on a uniform grid with J + 1 nodes. A non-uniform grid was unnecessary as the convergence tests showed. See Ivers and Phillips (2003, 2008) and Phillips and Ivers (2005) for the details. The spatial discretisation reduces the problem to a generalised algebraic eigenproblem [(A + fD) + RC] x = cBx, for the eigenpair (c,x), where the matrix A is banded block pentadiagonal, except for the last block-row affected by the boundary condition (A1). The eigenproblem was solved using generalised inverse iteration and the implicitly restarted Arnoldi method (Lehoucq et al., 1998), with continuation in the anisotropic parameter f to speed up the calculations. The associated critical value problem [(A + fD) + RcC]x = ixcBx, for the critical modified Rayleigh number Rc and angular frequency xc = Imc were solved using Newton’s method, noting Rec = 0 at criticality. Table 1 shows the convergence to the most unstable qe-modes for E = 104 at selected values of f and indicates that N = 20, J = 200 suffices for the present purposes. This is not inconsistent with the thickness of the Ekman boundary layer, which would suggest about J  E1/2 = 102 points. In the isotropic case f = 0 the critical modified Rayleigh numbers agree with Roberts (1968) and Ivers

and Phillips (2003); the critical angular frequencies agree with Ivers and Phillips (2003) (noting the use there of the semi-diurnal time scale for the angular frequencies). Of course, N must be increased with m to provide sufficient harmonics. Selected calculations suggest that N = 40, J = 400 is sufficient for E = 106 with Pr = 10. Appendix B. Isotropic thermal conduction in spheroids The anisotropic part of the numerical method, i.e. the spectral equations, the radial finite-differencing and the code itself, can be checked with analytic solutions for isotropic thermal conduction in an oblate or prolate spheroid V0 with the boundary S0 in the primed cartesian coordinates (x0 , y0 , z0 ),

x0 2 þ y 0 2 z 0 2 þ 2 ¼ 1; a2 b

ðB1Þ

where the temperature vanishes on the boundary. If the thermal conduction problem is non-dimensionalised using as a typical length scale L the semi-axes b⁄ or a⁄ in the oblate and prolate cases respectively and the thermal diffusion time scale L2 =j, the dimensionless form of the equations for the temperature H are 2

@ t H ¼ r0 H in V 0 ;

H ¼ 0 on S0 :

ðB2Þ

The problem is separable in confocal spheroidal variables (g, n, /) with prolate and oblate solutions (Niven, 1880), im/þct Hp ¼ Smn ðc; gÞRð1Þ ; mn ðc; nÞe

im/þct Ho ¼ Smn ðic; gÞRð1Þ ; mn ðic; inÞe

ð1Þ where Smn(c,g) and Rmn ðc; nÞ are prolate spheroidal angular and radial functions of the first kind (Flammer, 1957), and c = c2/ jb2  a2j. The degree n and order m are integers, 0 6 m 6 n. Confocal prolate (1 6 n < 1) and oblate (0 6 n < 1) spheroidal coordinates are related to cartesian coordinates by

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 d ð1  g2 Þðn2  1Þ cos /; 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 0 y ¼ d ð1  g2 Þðn2  1Þ sin /; 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x0 ¼

z0 ¼

1 dgn; 2

2

where d :¼ 2 jb  a2 j and 1 6 g 6 1. In either case the boundary S0 (B1) is given by n = n1:¼2b/d. The parameter c is determined from ð1Þ the boundary condition (B2) (b), Rmn ðc; n1 Þ ¼ 0 or Rð1Þ mn ðic; in1 Þ ¼ 0, and hence the growth rate c. The isotropic thermal conduction problem in the spheroid V0 can be transformed into an equivalent anisotropic thermal conduction problem in a sphere V by stretching the cartesian coordinates x = x0 /a, y = y0 /a and z = z0 /b. The boundary S0 is transformed to the sphere S, which has unit dimensionless radius. Eq. (B2) becomes

@H ¼ r  ðD  rHÞ; @t

ðB3Þ

where the anisotropic diffusion matrix D in the oblate and prolate cases, respectively, is

D¼

  1 e2 ; I þ 1 1 z z a2 1  e2

D¼

 1 I  e2 1z 1z a2

ðB4Þ

with e2:¼1  (b/a)2 or e2:¼1  (a/b)2, respectively. Expanding D yields the oblate and prolate equations in V (Ivers, 2005),

! @H 1 e2 @2 2 H; ¼ 2 r þ a @t 1  e2 @z2

! 2 @H 1 2 2 @ H: ¼ 2 r e a @t @z2

ðB5Þ

The factors of a2 can be removed by a new time scaling. Thermal conductivity is enhanced (diminished) in the z-direction by the anisotropy in the oblate (prolate) case. The boundary condition (B2)(b) becomes

D.J. Ivers, C.G. Phillips / Physics of the Earth and Planetary Interiors 190–191 (2012) 1–9

H ¼ 0 on r ¼ 1:

ðB6Þ

 ¼ 0; H ¼ H, Pr = a2 Eqs. (B5), (B6) are identical to (2.5), (2.8)(b) if v and f = e2/(1  e2) (oblate) or f = e2 (prolate). Since 0 6 e2 < 1,e2 = f/(1 + f) with 0 6 f < 1 (oblate) or e2 = f with 1 < f 6 0 pffiffiffiffiffiffiffiffiffiffiffi (prolate). In either case, f = (a/b)2  1 and a=b ¼ 1 þ f. Growth rates for the first z-even mode and first z-odd mode calculated for representative values of f are compared to the spheroidal wave function solutions in Table 2. In the oblate case 0 < f < 1, b = 1 and the azimuthal wave-number m = 1. The radial truncation is J = 1600 and the angular truncation N = 10 for the even mode, and for the odd mode if f 6 4. For the odd mode N = 20 if f > 4. In the prolate case 1 < f < 0,a = 1 and m = 1 for radial truncation J = 1600. The angular truncation is N = 10 for f > 0.95 and N = 20 for f = 0.95. References Braginsky, S.I., Meytlis, V.P., 1990. Local turbulence in the Earth’s core. Geophys. Astrophys. Fluid Dynam. 55, 71–87. Braginsky, S.I., 1993. MAC-oscillations of the hidden ocean of the core. J. Geomag. Geoelect. 45, 1517–1538. Buffett, B.A., 2003. A comparison of sub-grid scale models for large eddy simulations of convection in the Earth’s core. Geophys. J. Int. 153, 753–765. Busse, F.H., 1970. Thermal instabilities in rapidly rotating systems. J. Fluid Mech. 44, 441–460. Christensen, U.R., Olson, P., Glatzmaier, G.A., 1999. Numerical modelling of the geodynamo: a systematic parameter study. Geophys. J. Int., 138 393–409. Flammer, C., 1957. Spheroidal Wave Functions. Stanford University Press, Stanford, CA. Glatzmaier, G.A., Roberts, P.H., 1995. A three-dimensional convective dynamo with rotating and finitely conducting inner core and mantle. Phys. Earth Planet. Int. 91, 63–75. Ivers, D.J., 2002. Thermal instability of an oblate spheroid. In: Eighth SEDI (Studies of the Earth’s Deep Interior) Symposium, Granlibakken, USA (abstract). Ivers, D.J., Phillips, G.C., 2003. A vector spherical harmonic spectral code for linearised magnetohydrodynamics. ANZIAM J. 44 (E), C423–C442. Ivers, D.J., 2005. An angular spectral method for solution of the heat equation in spheroidal geometries. ANZIAM J. 46 (E), C854–C870.

9

Ivers, D.J., Phillips, G.C., 2008. Scalar and vector spherical harmonic spectral equations of rotating magnetohydrodynamics. Geophys. J. Int. 175, 955–974. Jones, C.A., Soward, A.M., Mussa, A.I., 2000. The onset of thermal convection in a rapidly rotating sphere. J. Fluid Mech. 405, 157–179. Kageyama, A., Sato, T., 1995. Computer simulation of a magnetohydrodynamic dynamo. Phys. Plasmas 2, 1421–1431. Kuang, W., Bloxham, J., 1997. An Earth-like numerical dynamo model. Nature 389, 371–374. Lehoucq, R.B., Sorensen, D.C., Yang, C., 1998. ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. SIAM, p. 160. Matsui, H., Buffett, B.A., 2005. Sub-grid scale model for convection-driven dynamos in a rotating plane layer. Phys. Earth Planet. Int. 153, 108–123. Matsushima, M., 2001. Expression of turbulent heat flux in the Earth’s core in terms of a second-order closure model. Phys. Earth Planet. Int. 104, 137–148. Matsushima, M., Nakajima, T., Roberts, P.H., 1999. The anisotropy of local turbulence in the Earth’s core. Earth, Planets, Space 51, 277–286. Niven, C., 1880. On the conduction of heat in ellipsoids of revolution. Phil. Trans. R. Soc. Lond. 171, 117–151. Phillips, G.C., Ivers, D.J., 2000. Spherical anisotropic diffusion models for the Earth’s core. Phys. Earth Planet. Int. 117, 209–223. Phillips, C.G., Ivers, D.J., 2003. Strong field anisotropic diffusion models for the Earth’s core. Phys. Earth Planet. Int. 140, 13–28. Phillips, C.G., Ivers, D.J., 2005. Spherical magnetic instabilities in the Earth’s core and equatorial symmetries. Phys. Earth Planet. Int. 153, 83–100. Roberts, P.H., 1968. On the thermal instability of a rotating-fluid sphere containing heat sources. Phil. Trans. R. Soc. Lond. A 263, 93–117. Sakuraba, A., Kono, M., 1999. Effects of the inner core on the numerical simulation of the magnetohydrodynamic dynamo. Phys. Earth Planet. Int. 111, 105–121. Soltis, T., Brestensky, J., 2010. Rotating magnetoconvection with anisotropic diffusivities in the Earth’s core. Phys. Earth Planet. Int. 178, 27–38. Soltis, T., Brestensky, J., Sevcik, S., 2004. Influence of thermal conductivity anisotropy on rotating magnetoconvection. In: Tenth SEDI (Studies of the Earth’s Deep Interior) Symposium, Garmisch-Partenkirchen, Germany (abstract). Soltis, T., Brestensky, J., Sevcik, S., 2006. Axisymmetric hydromagnetic instabilities and their possible relation to torsional oscillations in the Earth’s core. In: 10th SEDI (Studies of the Earth’s Deep Interior) Symposium, Prague, Czech Republic (abstract). St Pierre, M.G., 1996. On the local nature of turbulence in Earth’s outer core. Geophys. Astrophys. Fluid Dynam. 83, 293–306. Zhang, K., 1992. Spiralling columnar convection in rapidly rotating spherical fluid shells. J. Fluid Mech. 236, 535–556.