Planetary interiors:
Magnetic fields, Convection and Dynamo Theory
5. How numerical dynamo models are constructed
and what they produce
Chris Jones, Department of Applied Mathematics
University of Leeds UK

FDEPS Lecture 5, Kyoto, 30th November 2017

Section 5.
How numerical dynamo models are constructed
and what they produce
5.1 Spherical geodynamo models(5.1) Spherical geodynamo models

2/60

Spherical geometry Boussinesq dynamo models

Model for the geodynamo.
Spherical shell rotating about
z-axis. Gravity radially inward,g = g0 r . Centrifugal acceleration
small. Length scale d is
gap-width from inner to outerboundary. Convection onsets
outside tangent cylinder.
Tangent cylinder (TC) is imaginary cylinder touching the innercore. Convection columns inside the TC are divided.
Boussinesq fluid. Sometimes consider basal heating through theinner core boundary (ICB), sometimes internal heating.

(5.1) Spherical geodynamo models

3/60

Spherical dynamo equationsDu1+ 2Ω × u = −∇p + j × B + ν∇2 u + g0 αθr,Dt
ρ

(5.1.1)∂B= η∇2 B + ∇ × (u × B),∂t(5.1.2)∂θ+ u · ∇θ = κ∇2 θ + β(r )ur + H,∂t

(5.1.3)∇ · B = ∇ · u = 0,µj = ∇ × B.

(5.1.4, 5, 6)β(r ) - Basic state temperature profile, H is heat source in the core.
No-slip and stress-free boundaries considered, but assume θ = 0
there.

(5.1) Spherical geodynamo models

4/60

Basic state in dynamo models
Basic state temperature satisfies
κ∇2 T̄ + H = 0,

1 d d T̄
κ 2 r2+ H = 0,r dr dr(5.1.7)so

C2 Hr 2−(5.1.8)r
6κ
where C1 and C2 are determined by the boundary conditions on T ,e.g. for constant temperature boundaries T̄ = Ti at r = ri andT̄ = Ti − ∆T at r = ro .
T̄ = C1 +Then the temperature gradient β(r ) = −d T̄ /dr . The fulltemperature T = T̄ + θ, so fixed temperature boundaries have
θ = 0 at both boundaries. Alternatively, fixed flux conditions have−κd T̄ /dr = F a prescribed flux, and then dθ/dr = 0.
(5.1) Spherical geodynamo models

5/60

Dimensionless equations
Non-dimensionalise with d = ro − ri as unit of length, d 2 /η as unitof time, (Ωρηµ)1/2 as the unit of magnetic field, κ∆T /η as unit oftemperature, or Fd/κ in fixed flux case.
E Du
ERaPm
+ 2ẑ × u = −∇p + j × B + E ∇2 u +θr, (5.1.9)
Pm Dt
Pr
∂B
= ∇2 B + ∇ × (u × B),∂t(5.1.10)∂θPm 2=∇ θ − u · ∇θ + β(r )ur ,∂tPr(5.1.11)∇ · B = ∇ · u = 0,j = ∇ × B.

(5.1.12, 13, 14)β(r ) - Basic state temperature profile. No internal heating
assumed here.
(5.1) Spherical geodynamo models

6/60

Dimensionless parameters

Ekman number E = ν/Ωd 2 , Core value ∼ 10−15 .
g0 αβd 3
Rayleigh number Ra =, Core value very large.
νκ
Prandtl number Pr = ν/κ, Core value ≈ 0.1.
Magnetic Prandtl number Pm = ν/η,Core value ∼ 10−5 .

η magnetic diffusivity, ν kinematic viscosity,κ is the thermal diffusivity, d = ro − ri .
Even if the molecular diffusivities are enhanced by turbulenceE ∼ 10−9 .
Typical simulation values are Pm = Pr = 1, E = 5 × 10−5 , small
but not so small as to make numerical simulation inconvenient.

(5.1) Spherical geodynamo models

7/60

5.2 Pseudo-spectral method for Boussinesq equations(5.2) Pseudo-spectral method for Boussinesq equations

8/60

Poloidal-Toroidal expansion
Poloidal-Toroidal decompositionu = ∇ × T r + ∇ × ∇ × Pr.

(5.2.1)This guarantees that ∇ · u = 0.
We use (5.1.11) fror the temperature equation, and expand themagnetic field as
B = ∇ × T r + ∇ × ∇ × Pr.

(5.2.2)so ∇ · B = 0 exactly.
So we now have 5 scalar functions, P, T , P, T and θ to integrate
forward in time.
The three components of u and B have been reduced to two by
this expansion.
(5.2) Pseudo-spectral method for Boussinesq equations

9/60

Expansion in spherical harmonics
Each of these 5 scalars is expanded in spherical harmonics, e.g.
P=

`=L X
X

P`m (r , t)Y`m (θ, φ)(5.2.3)`=0 |m|≤`where Y`m are the Schmidt normalised spherical harmonics. L is
the truncation parameter, typically L = 128 or L = 256.
Nowr̂ · ∇ × ∇ × rP = −r ∇2H P =`=L X
X`(` + 1)`=0 |m|≤`rP`m (r , t)Y`m (θ, φ)(5.2.4)so each individual component is just multiplied by `(` + 1)/r , andthe components are not coupled.
(5.2) Pseudo-spectral method for Boussinesq equations

10/60

Deriving the scalar equations 1.
The radial component of the induction equation (5.1.10) and its
curl give two equations for P and T .
The induction equation is written as(∂t − ∇2 )B = NB = ∇ × (u × B).

(5.2.5)

NB is the nonlinear term, the other terms are linear. When we
insert the expansion (5.2.3) into this and take the radial
component, the toroidal term gives nothing, so for each `, mcomponent we can write(∂t − ∇2 )P =rr̂ · NB .
`(` + 1)(5.2.6)

So if we can evaluate the ` and m components of the nonlinear
term, we can use (5.2.6) to timestep P in a simple way.
(5.2) Pseudo-spectral method for Boussinesq equations

11/60

Deriving the scalar equations 2.
The radial component of the curl of the induction just gives a very
similar equation for T , only the nonlinear term is more
complicated.
The same method works for the toroidal components and poloidal
components of the velocity.
We write (5.1.9) in the form(∂t − Pm∇2 )u = Nu ,(5.2.7)

where the Coriolis and buoyancy terms are included in the Nu term.
We then insert the expansion (5.2.1), and take the radial
component of the curl and double curl of this equation.
As with induction equation, the P and T equations separate out,
making it convenient to time-step them forward.
(5.2) Pseudo-spectral method for Boussinesq equations

12/60

Influence matrix method 1.
There is a slight difficulty with the P equation, as this is fourthorder in r , so we get two equations(∂t − Pm∇2 )P = g ,−∇2 g =

rr̂ · ∇ × ∇ × Nu . (5.2.8)`(` + 1)

However these can be solved simply by using the influence matrix
method described in Peyret’s book.
Each equation is second order in r . There are two boundary
conditions on T , four on P and none on g .
The coupled system for P and g is solved by means of Greens
functions with the no-slip condition on P

X P̄ = ḡ , (ri < r < ro ) ∂r P̄ = 0, (r = ri , r0 )(5.2.9)Q ḡ = f ,ḡ = 0.
We solve (5.2.9) for P̄ and ḡ at each time-step.
(5.2) Pseudo-spectral method for Boussinesq equations

13/60

Influence matrix method 2.
We then solve the time independent equations
X PG = gG ∂rPG = 0
Q gG = 0
gG = 1, 0.


X PG0 = gG0
Q gG0 = 0

∂rPG0 = 0
gG0 = 0, 1.

(5.2.10)(5.2.11)which can be pre-computed. We then add on the linear
combination of PG and PG0 , P = P̄ + aPG + bPG0 , which satisfies
the boundary conditions P = 0 at r = ri , ro ,

 


PG (ri ) PG0 (ri )aP̄(ri )=−
(5.2.12)P̄(ro )
PG (ro ) PG0 (ro )b
P now satisfies the desired boundary conditions P = ∂rP = 0 at
both boundaries.
(5.2) Pseudo-spectral method for Boussinesq equations

14/60

Radial dependence
We need to represent the radial dependence numerically. Some
codes use a high order finite difference non-uniform mesh over r
typically up to 160 points. Other codes expand P`m in Chebyshev
polynomials.
All differentiation done in spectral space, all nonlinear
multiplications in physical space.
At each time-step, the quantities needed to evaluate the nonlinear
terms are evaluated on a grid in r , θ and φ space. The mesh isused is bigger than the number of spectral coefficients used. If
0 ≤ ` ≤ L, the number of θ points is 3L/2, a choice known as thede-aliasing rule.
Once the nonlinear terms are evaluated in physical space they canbe converted back to spectral space, i.e. expansions in sphericalharmonics.
(5.2) Pseudo-spectral method for Boussinesq equations

15/60

Time-step methods

We can now advance the integration one time-step. The linear
terms are advanced implicitly using the Crank-Nicolson scheme,
the nonlinear terms explicitly using Adams-Bashforth. A
predictor-corrector method is used to choose the time-step.
In practice, it is more convenient to treat some of the linear terms
explicitly, i.e. treat them the same way as the nonlinear terms. In
particular, the Coriolis terms are usually treated explicitly.
Slowest part of the code is evaluating the Legendre transform togo from the spectral representation to the grid and back. Modern
codes have optimised this and so are faster.
A fast Fourier transform can be used in φ.

(5.2) Pseudo-spectral method for Boussinesq equations

16/60

Implementing the magnetic boundary conditions
No-slip boundaries imply T = 0, P = 0 and dP/dr = 0 at theboundaries. Stress-free conditions are also possible.
Insulating magnetic boundary conditions are the most common,T`m = 0 on r = ri , ro ,(5.2.13)

∂P`m
P`m
−`
= 0 on r = ri
(5.2.14)
∂r
r
P`m
∂P`m
+ (` + 1)= 0 on r = ro .
(5.2.15)∂rr
these are derived by matching the field to a potential field inside
the inner core and outside the CMB. More complicated options are
possible.

(5.2) Pseudo-spectral method for Boussinesq equations

17/60

5.3 Results from Boussinesq dynamo codes(5.3) Results from Boussinesq dynamo codes

18/60

Results from dynamo codes

Pr = Pm = 1, Ra = 750, E = 10−4 .
Radial magnetic field snapshot at the CMB
No internal heating. No-slip, fixed temperature, insulatingboundaries.
Very dipolar, doesn’t reverse. Field slightly weaker at the poles.
Intense flux patches at high latitudes.
(5.3) Results from Boussinesq dynamo codes

19/60

Velocity field

Pr = Pm = 1, Ra = 750, E = 10−4 .
Radial velocity snapshot at
r = 0.8rCMB .
Note the columnar nature of the convection rolls, local Rossby
number small.
Intense flux patches at the top of these columnar rolls.
Pattern propagates westward.
(5.3) Results from Boussinesq dynamo codes

20/60

Higher local Rossby number

Pr = Pm = 0.2, Ra = 750, E = 10−4 .
Radial magnetic field snapshot at the CMB
Much less dipolar. Field strength is weaker.
This type of dynamo can reverse.
Rossby number is larger, and inertia is playing a significant role.

(5.3) Results from Boussinesq dynamo codes

21/60

Flow at higher local Rossby number

Pr = Pm = 0.2, Ra = 750, E = 10−4 .
Radial velocity snapshot at
r = 0.8rCMB .
More activity near the poles,
less columnar convection rolls.
Between these patterns lies an Earth-like regime.
(5.3) Results from Boussinesq dynamo codes

22/60

Pm

Pm

Variation with Ekman number

dipolar non−dip
dipolar
10
10
5
5
2
2
1
1
0.5
0.5
0.2
0.2
0.1 E=10−3
0.1 E=10−4
.05
.05
1 2
5 10 20 50
1 2
5 10 20

dipolar
10
5
2
1
0.5
0.2
0.1 E=10−5
.05
50
1 2
5 10 20 50

dipolarnon−dip 10
dipolar
10
5
5
2
2
1
1
0.5
0.5
0.2
0.2
0.1 E=3×10−4
0.1 E=3×10−5
.05
.05
1 2
5 10 20 50
1 2
5 10 20

dipolar
10
5
2
1
0.5
0.2
0.1 E=3×10−6
.05
50
1 2
5 10 20 50

Ra/Racrit

Ra/Racrit

Ra/Racrit

Christensen and Aubert 2006 summarise the dependence on
Ekman number E and Pm.
At low E , low Pm dynamos are possible, provided Ra is large
enough. Important, as liquid metals have low Pm.

(5.3) Results from Boussinesq dynamo codes

23/60

Small E , low Pm dynamo, Flow pattern

Pr = 1, Pm = 0.1, Ra = 50Racrit , E = 3 × 10−6 .
Left: radial velocity at r = 0.5r0 . Right: radial velocity at
r = 0.8r0 .
At lower E the convective columns are much thinner, particularly
further out from the ICB.
(5.3) Results from Boussinesq dynamo codes

24/60

Small E , low Pm dynamo, Scale separation

Left: radial magnetic field at z = 0.2. Right: vorticity at z = 0.2.
Notice that the magnetic field is on a much larger scale in these
low Pm calculations. Temperature fluctuations on same scale asmagnetic field.

(5.3) Results from Boussinesq dynamo codes

25/60

Earth-like geodynamo models
Christensen et al. 2010.

(a) Geomagnetic field. (b) Earth-like geodynamo model with
E = 3 × 10−5 , Ra = 3 × 108 , Pm = 2.5, Pr = 1. Zero flux on theouter boundary, compositional driving.
Case (c) is not sufficiently dipolar, (d) is too dipolar, (e) is toostrong near the poles.
(5.3) Results from Boussinesq dynamo codes

26/60

Wedge for Earth-like geodynamo models

Earth-like geodynamo models have Rm and Eη = E /Pm that lie ina particular wedge in the Rm − Eη plane. Earth may lie in thiswedge too.
(5.3) Results from Boussinesq dynamo codes

27/60

Evolution of magnetic energy with time
Sreenivasan and Jones 2011.
(ii)

106
10

4

Em

102
10

(i)

0

(iii)
(iv)

(vi)

10-2
10-4

(v)

10-6
10-8
0.0

0.5

1.0
t

1.5

2.0

E = 10−4 , Pr = Pm = 1 and stress-free boundary conditions. (i)Ra = 400, dipole dominated solution. (ii) Ra = 600 strong field
dipole dominated solution. (iii) Ra = 600, initial small field
solution grows into a relatively weak quadrupolar solution. (iv) as
case (iii) but with a different small initial field. (v) Ra = 500, an
initial small field decays away. (vi) Ra = 550, an small initial field
grows, eventually resulting in a quadrupolar field(5.3) Results from Boussinesq dynamo codes

28/60

Helicity of the strong and weak field solutions

6.0e+05

6.0e+05

4.0e+05

4.0e+05

2.0e+05

2.0e+05

0.0e+00

0.0e+00-2.0e+05-2.0e+05-4.0e+05-4.0e+05-6.0e+05-6.0e+05

E = 5 × 10−5 , Pr = Pm = 1, with stress-free boundaries. Helicity
is antisymmetric about equator.
Left: is the helicity u · ζ for the strong field solution on the planez = 0.5, Ra = 400. Right: is the helicity without magnetic field at
Ra = 500. At Ra = 500 an initial seed field decays.
Note that in the pure convection case, there is helicity, but little
net helicity. With magnetic field, there is far more net helicity, as
predicted by linear theory.
(5.3) Results from Boussinesq dynamo codes

29/60

Scaling laws: magnetic field strength
In the Earth’s core, magnetic energy is much greater than kinetic
energy, so a simple balance as used in astrophysics won’t work here.
We start with
Ohmic dissipation + Viscous dissipation = rate of working of
buoyancy forces,
The rate of working work done by the buoyancy forces isZ
ρg αT 0 ur dv
which can be written in terms of the heat flux.
Ignore the viscous dissipation and we obtainZZ
g αFconv
dv ∼ ηµj2 dv
cp
Now need to relate magnetic energy to magnetic dissipation.
(5.3) Results from Boussinesq dynamo codes

30/60

Christensen-Tilgner law
The dissipation time is the time taken for the magnetic energy to
be dissipated through ohmic loss,
ZZτdiss ηµj2 dv = B2 /2µ dv
Equivalently, the magnetic dissipation length

δB =

τdiss
η

1/2
.

Christensen and Tilgner (Nature, 2004) proposed that
δB ∼ dRm−1/2
mainly on the basis of simulations and laboratory experiments. It
also has some theoretical support, because at high Rm flux ropesof this thickness are formed.
(5.3) Results from Boussinesq dynamo codes

31/60

Predicted field strength
We now have
ηµj2 ∼ η(∇ × B)2
ηB2
g αF
∼ 2 ∼
,
µ
cp
µδB

givingB∗ ∼g αFconv µd
U∗ cp

1/2
.

With the inertial theory scaling for U∗ this givesB∗ ∼ µ1/2 2/5 1/5

d

ρ

Ω

1/10



g αF
cp

3/10
.

Remarkable feature is weak dependence of B on Ω. We areassuming though that the planet is in the rapidly rotating low Ro
regime.
(5.3) Results from Boussinesq dynamo codes

32/60

5.4 Compressible convection equations(5.4) Compressible convection equations

33/60

Compressible convection equations 1.
In compressible flow, mass conservation∂ρ+ ∇ · (ρu) = 0∂t(5.4.1)and this replaces ∇ · u = 0, the incompressible flow equation.
The momentum equation is unchanged,ρDu+ 2ρΩ × u = −∇p + gρ + j × B + Fν
Dt

(5.4.2)

except that now the density is no longer constant, but a variable tobe determined. Fν is the viscous force. B is given by (5.1.2), theusual induction equation.
We also need an equation of state giving the pressure in terms ofthe temperature and density p = p(ρ, T ) and an equation for theentropy, S = S(p, ρ, T ). We need to use entropy as it appears in
the compressible heat transfer.
(5.4) Compressible convection equations

34/60

Compressible convection equations 2.
For a perfect gas,p = RρT ,S = cv lnpT= cv ln (γ−1) + const,γρρ

(5.4.3, 4)

where R is the gas constant, S is the entropy and cv and cp are
the specific heats at constant volume and constant pressure,respectively. γ = cp /cv .
Giant planets are not perfect gas, because the hydrogen can be
metallic. But the perfect gas case helps for understanding tbe
effects of compressibility.
The entropy equation is∂S+ (u · ∇)S = ∇ · ρcp κm ∇T + Qj + Qν .
(5.4.5)ρT∂t
To get to the Boussinesq equation, entropy S is replaced bytemperature, ρ and T take their constant reference values and therate of heating by viscous dissipation Qν and ohmic dissipation Qj
are neglected.
(5.4) Compressible convection equations

35/60

Compressible convection equations 3.
This is now a complete set of equations. If we start with an initialcondition for all the variables, we use (5.4.1), (5.4.2) and (5.4.5)
to evolve ρ, u and S. Then (5.4.4) determines the temperature T ,and (5.4.3) determines the pressure p.
Some numerical simulations of the dynamo equations do work like
this, but there are difficulties.
The term ∂ρ/∂t in (5.4.1) means that sound waves are present in
the system. To see this look at small perturbations about theconstant entropy, uniform state case.
Sound waves are typically very fast, so we need a very small
timestep to integrate the equations, which is not efficient. We use
the anelastic equations to avoid this problem.

(5.4) Compressible convection equations

36/60

5.5 Anelastic convection equations(5.5) Anelastic convection equations

37/60

Anelastic models
Compressible systems that are convecting are often close to an
adiabatic state.
The anelastic approximation is based on assuming thermodynamic
quantities are close to their adiabatic values. Sometimes called
thermodynamic linearisation.
p = p̄(r ) + p 0 ,ρ = ρ̄(r ) + ρ0 ,T = T̄ (r ) + T 0 ,(5.5.1)

where p̄, ρ̄ and T̄ are (in spherical geometry) functions of radius
only, and the primed quantities are small compared to thereference state values.
The reference state is chosen to satisfy the hydrostatic equationd p̄/dr = −g ρ̄. For a centrally condensed object whereg = GM0 /r 2 , a polytropic reference state is a popular choice.

(5.5) Anelastic convection equations

38/60

Polytropic Reference StateThe reference state is thenp̄ = p0 ζ n+1 ,ρ̄ = ρ0 ζ n ,T̄ = T0 ζ,

ζ=

c1
+ c0
r

(5.5.2)where c1 and c0 are constants, which satisfies the hydrostaticbalance equation.
A polytropic index n ≈ 1/(γ − 1) ensures that the reference stateis close to adiabatic.
Other more complicated, but more realistic models are possible.
If ri and ro are the radius of the inner boundary (ICB) and outerboundary (CMB) the ri /ro is the radius ratio.
 ρiNρ = ln(5.5.3)ρomeasures the density ratio across the spherical shell.
(5.5) Anelastic convection equations

39/60

Anelastic approximation 1.
The basic assumption is that the entropy difference across the layeris small, i.e.
S=  << 1.
(5.5.4)cp
Since for perfect gas
p0ρ0 T 0
= +
,
p̄
ρ̄
T̄

S
p0
ρ0
=
− ,
cp
γp̄
ρ̄

p0
ρ0
T0
∼
∼
∼ O().
p̄ρ̄T̄
(5.5.5)Equation of motion contains terms of O(),Du+ 2ρ̄Ω × u = −∇p 0 + gρ0 + j × B + Fν ,(5.5.6)Dt√√
so for consistency, |u| ∼ O(1/2 ) gd, where gd is the free fallvelocity across the shell. The Alfvén speed must also be small
compared to the free fall velocity.
ρ̄(5.5) Anelastic convection equations

40/60

Anelastic approximation 2.

Now consider mass conservation (5.1.1),∂ρ0+ ∇ · (ρ̄u) = 0,∂t(5.5.7)

since ρ̄ is independentp√ of time.
Since |u| ∼ O(1/2 ) gd, and ∂/∂t ∼ O(1/2 ) g /dand ρ0 /ρ̄ ∼ O(), ∂ρ0 /∂t is very small compared to ∇ · (ρ̄u) = 0.
So mass conservation equation is simply∇ · (ρ̄u) = 0,(5.5.8)in the anelastic approximation. The sound wave term is removed.

(5.5) Anelastic convection equations

41/60

Anelastic approximation 3.

The anelastic equations are now∇ · (ρ̄u) = 0,Duρ̄+ 2ρ̄Ω × u = −∇p 0 + gρ0 + j × B + Fν ,
Dt


∂S
+ (u · ∇)S = ∇ · ρ̄cp κm ∇T 0 + Qj + Qν ,
ρ̄T̄∂t

(5.5.9)(5.5.10)(5.5.11)and these, together with the induction equation (5.1.2), can besolved more easily than the full compressible equations.
There are two further modifications to get the system into its mostfrequently used form.

(5.5) Anelastic convection equations

42/60

Entropy diffusion

What is the effect of small scale turbulence on the entropy
equation?
In Boussinesq convection, turbulence enhances the molecular
diffusion. The form of the equation remains the same though. We
just have κT ∇2 T , with a turbulent value of κ, instead of κ∇2 Twith a laminar value of κ.
In anelastic convection, it is different. We use the gradientdiffusion model, where a small scale ut induces a small-scale S 0 .

(5.5) Anelastic convection equations

43/60

Gradient diffusion model

We assume a small-scale turbulence, ut which gives rise to a
turbulent entropy fluctuation S t . S t is forced by the termρ(ut · ∇)S so its proportional to ∇S. The turbulent entropy flux isthen∂Sρut S t = It = −ρκij.
(5.5.12)∂xj
We now make the usual assumption that κij is an isotropic tensor,
κt δij . So
It = −ρκt ∇S.
(5.5.13)(5.5) Anelastic convection equations

44/60

Heat transport with turbulent diffusionl
Enhanced diffusion of entropy also gives a source term of entropy,
which turns out to beσt = −

1 t(I · ∇)T
T(5.5.14)the only form consistent with energy conservation. So, we get


∂S
ρ̄T̄
+ (u · ∇)S = ∇ · ρ̄T̄ κt ∇S + ∇ · ρ̄cp κm ∇T 0 + Qj + Qν ,∂t
(5.5.15)
that is turbulent diffusion of entropy and molecular diffusion oftemperature.
In compressible convection, the turbulence does not just enhance
the molecular diffusion of temperature, it diffuses entropy instead.

(5.5) Anelastic convection equations

45/60

Second law of thermodynamics?
Landau and Lifshitz’s equation is compatible with the second lawof thermodynamics, entropy always increases. Is this true of the
turbulence model?
Source is1
σ t = − (It · ∇T )T
so1
σ t = κt (∇S · ∇T )(5.5.16)
T
Entropy gradient and temperature gradient both point outward, so
generally OK.
For a stably stratified layer, the entropy gradient and temperature
gradient can be opposite, and then entropy diffusion becomes
invalid.
(5.5) Anelastic convection equations

46/60

Lantz-Braginsky-Roberts approximation

The terms1
ρ0− ∇p 0 + gρ̄ρ̄

can be rewritten as
 0


 0
ρ0
p
S p 0 1 d p̄ 1 d ρ̄
p
S
1
0
−g +
−
≈∇
−g
− ∇p +g = ∇ρ̄ρ̄ρ̄
cp p̄ γp̄ dr
ρ̄ dr
ρ̄
cp
(5.5.17)
the last term vanishing because the reference state is close toadiabatic. The momentum equation can now be written
 0
Du
p
S
j × B Fν
+ 2Ω × u = −∇+.
(5.5.18)
−g +
Dt
ρ̄
cp
ρ̄
ρ̄

(5.5) Anelastic convection equations

47/60

Viscous and magnetic terms
The magnetic field B and the current j are governed by thestandard induction equation: compressibility makes no difference.
The viscous terms do change,


∂uj
∂ui
2 ∂ ∂uj∂µ+µ,(5.5.19)−Fν,i =∂xj∂xj∂xi
3 ∂xi ∂xj
where µ = ρ̄ν is the dynamic viscosity and ν is the kinematic
viscosity. The models below are constant kinematic viscosity, butsometime constant dynamic viscosity is used. The viscous
dissipation rate is∂uj∂ui ∂ui2Qν = µ+− δij ∇ · u .
(5.5.20)∂xj ∂xj∂xi
3
and the ohmic dissipation rate is
QJ = µηj2 .

(5.5.21)

The anelastic equations are non-dimensionalised in the same way
as the Boussinesq equations.
(5.5) Anelastic convection equations

48/60

5.6 Pseudospectral method for anelastic convection(5.6) Pseudospectral method for anelastic convection

49/60

Pseudospectral method for solving the anelastic equations 1
Poloidal-Toroidal decompositionρ̄u = ∇ × (ρ̄T )r + ∇ × ∇ × (ρ̄P)r.
This guarantees that ∇ · ρ̄u = 0.
Equations for T and P given by r̂ · ∇× momentum equation andr̂ · ∇ × ∇× momentum equation, giving two equations for the two
unknowns T and P.
We use (5.5.15) with no temperature diffusion for the entropy
equation, and expand the magnetic field as
B = ∇ × T r + ∇ × ∇ × Pr.
so ∇ · B = 0 exactly. The radial component of the induction
equation and its curl give two equations for T and P.
(5.6) Pseudospectral method for anelastic convection

50/60

Pseudospectral method for solving the equations 2.
We now have 5 equations for the 5 unknowns T , P, S, T and P,
which are closed because of the entropy diffusion ansatz.
P
Expand all variables as S = S`m (r , t)P`m (cos θ)e imφ andsubstitute these expansions into the equations. Typically need `from 0-192.
Finite difference non-uniform mesh over r typically up to 160
points. All differentiation done in spectral space, all nonlinear
multiplications in physical space.
At each time-step, the quantities needed to evaluate the nonlinear
terms are evaluated on a a grid in r , θ and φ space. The mesh isused is bigger than the number of spectral coefficients used. If
0 ≤ ` ≤ L, the number of theta points is 3L/2, a choice known as
the de-aliasing rule.
Once the nonlinear terms are evaluated in physical spce they canbe converted back to spectral space, i.e. expansions in sphericalharmonics.
(5.6) Pseudospectral method for anelastic convection

51/60

5.7 Anelastic dynamo benchmark(5.7) Anelastic dynamo benchmark

52/60

Anelastic Dynamo Benchmark

To test the various codes, the dynamo community ran an anelastic
benchmark case to see if they gave the same results.
The case E = 10−3 , Ra = 351, 806, Pr = 1, Nρ = 5, ri /ro = 0.35,
polytropic index n = 2 was selected, as this gives simple solutions
and is not very demanding computationally.
The hydrodynamic case with the field switched off was tested, andthen magnetic field put back in.

(5.7) Anelastic dynamo benchmark

53/60

Hydrodynamic benchmark (no magnetic field)
equatorial
plane
ur

spherical
surfaceentropy0.93r0meridional
slice
ur

meridionalslice
zonal uφE = 10−3 , Ra = 351, 806, Pr = 1, Nρ = 5, β = 0.35, n = 2
Solution steady in a drifting frame, has m = 19 symmetry. Critical
Ra for onset of convection Ra = 283, 175, with m = 20.
(5.7) Anelastic dynamo benchmark

54/60

Results of the hydrodynamic benchmark

E = 10−3 , Ra = 351, 806, Pr = 1, Nρ = 5, β = 0.35, n = 2
Code:
Leeds
Glatzmaier
MAGIC
KE
81.858181.8335
81.8385
Luminosity
4.19886
4.19886
4.19876
Zonal KE9.37724
9.374379.37514Meridional KE 0.0220183
0.0220109
0.0220136
Resolution 128 × 192 × 384: 121 × 512 × 1024: 121 × 192 × 384
Timestep
2.5 × 10−6
5.8 × 10−6
5 × 10−6
Agreement to within 0.05%, very satisfactory. Solutions drift
eastwards, period 0.0187. Need small timestep.

(5.7) Anelastic dynamo benchmark

55/60

Steady dynamo benchmark
Steady means steady in a drifting frame. These solutions are very
difficult to find. Flow has to have large Rm for a dynamo, butlowish Re to be hydrodynamically stable. Large Pm only
possibility, but this is numerically difficult. Needed to compromise
on Nρ and E .
E = 2 × 10−3 , Ra = 80, 000, Pr = 1, Pm = 50, Nρ = 3, β = 0.35,
n=2
The dynamo is supercritical, Racrit = 61, 621.682, with azimuthal
wavenumber m = 10.
However, the benchmark dynamo has exact m = 7 symmetry.
Runs did not assume this, but used all wavenumbers, to ensurestability against perturbations of any m.
This compressible dynamo drifts eastwards. The Boussinesq
Christensen benchmark drifts westward, and much slower.

(5.7) Anelastic dynamo benchmark

56/60

Dynamo Benchmark energy
106E

Kinetic energy

Magnetic energy

105
0.0

0.5

1.0

1.5

2.0

2.5

t

E = 2 × 10−3 , Ra = 80, 000, Pr = 1, Pm = 50, Nρ = 3, β = 0.35,
n=2Time here is measured in magnetic diffusion times, so t=2
corresponds to 100 viscous diffusion times. Timestep is 3e-06 (orless) on magnetic diffusion timescale. Fortunately, the solution
does not require very high resolution: 64 radial points and 48 X 48
spherical harmonics gives well-resolved solutions.
(5.7) Anelastic dynamo benchmark

57/60

Steady Dynamo Benchmark Flowequatorial
plane
ur

spherical
surfaceentropyr = 0.86romeridional
slice
ur

meridionalslice
zonal uφE = 2 × 10−3 , Ra = 80, 000, Pr = 1, Pm = 50, Nρ = 3, β = 0.35,
n=2(5.7) Anelastic dynamo benchmark

58/60

Steady Dynamo Benchmark Field

br viewed
from
N. Poleequatorial
slice
bθ

spherical
surfacebrr = 0.86romeridionalslice
zonal bφ

E = 2 × 10−3 , Ra = 80, 000, Pr = 1, Pm = 50, Nρ = 3, η = 0.35,
n=2(5.7) Anelastic dynamo benchmark

59/60

Steady Benchmark Comparisons

Code:
Leeds
Glatzmaier
MAGIC
ASH
KE
4.19407 × 105 4.19390 × 105 4.19438 × 105
4.19106 × 105
ME3.20185 × 105 3.20083 × 105 3.19749 × 105
3.17330 × 105
Luminosity 11.5030
11.5030
11.5030
11.5033
Resolution 128X144X252: 65X128X256: 65X128X256:
97X256X512
Timestep
10−6
2.1 × 10−6
10−6
−6
1.4 × 10
Generally OK, though ASH is a bit of an outlier.

(5.7) Anelastic dynamo benchmark

60/60