An application of a mapping method to asymmetric kinematic dynamos

~,

4<

,.

PHYSICS O F T t | E EARTH ANDPLANETARY INTERIORS

+."

ELSEVIER

Physics of the Earth and Planetary Interiors 91 (1995) 53-61

An application of a mapping method to asymmetric kinematic dynamos Takahiro Nakajima, Paul H. Roberts * School of Mathematics and Statistics, University of Sydney, Sydney, N.S. I,E. 2006, Australia

Received 8 November 1994; revision accepted 28 February 1995

Abstract

The Earth's magnetic field is generated by electromagnetic induction in the core, where the Ekman and Rossby numbers are small. In consequence, the magnetostrophic equation should be an excellent approximation to the equation of motion. In an earlier paper, we argued that the simultaneous integration of the magnetostrophic and induction equations would be simplified if a certain mapping method were employed. In this method, the grid points of one set of coordinates in the meridional plane lie on lines parallel to the rotation axis, and one curve of the other (non-orthogonal) set coincides with the surface of the core. By this arrangement, not only is the non-local condition on the magnetic field at the surface of the core easily implemented but also the magnetostrophic equation is readily solved without interpolation and to second-order accuracy. The mapping method was applied in that paper to axisymmetric dynamos, not only to kinematic a 2- and ato-models but also to a dynamic (intermediate) dynamo of model-Z type. The results, especially for the kinematic models, were in good agreement with earlier integrations that used spectral methods. In this sequel to that paper, we generalize the mapping method to kinematic dynamos in three dimensions. We suppose that the velocity is axisymmetric but that the magnetic field is proportional to exp(im~b + At), where ~b is longitude and t is time. We recover the growth rates A for the decay modes, and find the critical magnetic Reynolds number, Rm, for dynamo models proposed by other workers, the critical magnetic Reynolds number being the smallest value of the Reynolds R (or of - R) for which cr = ~ A is zero. These values of R m agree well with values derived in earlier studies.

1. Introduction

This p a p e r is c o n c e r n e d with asymmetric kinematic dynamos, i.e. self-sustaining solutions o f

* Corresponding author present address: Institute of Geophysics and Planetary Physics, University of California, Los Angeles, CA 90095, USA. Elsevier Science B.V. SSDI 0031-9201 (95)03 048 -4

the induction equation, which in dimensionless form is b B / O t = R V × ( V × B ) + V2B,

(1)

where the fluid velocity, V, is steady and completely axisymmetric, whereas the magnetic field, B, is purely asymmetric, being proportional to exp(im~b + A t ) , where ~b is longitude, m is a n o n z e r o integer and A = tr + il~ is the (complex) growth rate. R, the magnetic Reynolds number, is

T. Nakajima, P.H. Roberts/Physics of the Earth and Planetary Interiors 91 (1995) 53-61

54

a non-dimensional measure of IV [. The existence of self-sustaining (o" >_0) solutions of this type was first established by G.O. Roberts, and reported by one of us (Roberts, 1971). Later, Gubbins (1973) worked on a second model, and more recently still Dudley (1988) and Dudley and James (1989) presented detailed results from their investigations of three further axisymmetric models of a rather simple nature. The existence of a preferred axis, Oz, is suggested by the important role played by Coriolis forces in the dynamical equations, when Oz is parallel to the angular velocity of the Earth. The full magnetohydrodynamic (MHD) dynamo problem, in which the induction equation (1) is solved simultaneously together with V.B

= 0

(2)

and the equation of fluid motion (including the Lorentz force), is not examined in the present paper. It is, however, pertinent to remark that the Coriolis force is presumed to be so potent in the core that the full equation of fluid motion is often replaced by the magnetostrophic approximation, in which the inertial and viscous forces are ignored, the dominant balance being between Corioils, Lorentz, buoyancy and pressure forces. The magnetostrophic equation is most easily integrated in cylindrical coordinates (s, ~b, z), but the geometry of the system makes (1) and (2) most easily soluble in spherical coordinates (r, 0, ~b), where r = I r l (where r, is the radius vector from the geocenter) and ~9 is colatitude, it being supposed that the electrically conducting core, V, is spherical, and of unit radius, in our non-dimensional units. For this reason, several studies (starting with Braginsky, (1978)) have made use of both systems of coordinates simultaneously. The equations are stepped forward in time. After the magnetic field has been advanced by integrating (1) and (2) in spherical coordinates, it is interpolated to a cylindrical grid. The magnetostrophic equation is integrated in that coordinate system to obtain the fluid velocity, II, which is then interpolated back to the spherical grid in preparation for the next time step. In a recent paper (Nakajima and Roberts (1995), hereafter referred to as 'Paper 1') we

have explored an alternative procedure in which the governing equations are mapped onto a nonorthogonal grid that combines the virtues of the cylindrical and spherical systems: one set of coordinate lines are parallel to Oz so that the magnetostrophic equations are easily integrated; the other set contains one member that coincides with the surface, S, of the conductor, i.e. r = 1. More precisely, we introduce coordinates (s', ~b, z'), where

s =s(s'),

z =z(s', z')

(3)

and rewrite the governing equations and boundary conditions in terms of s' and z'. In Paper 1, we presented results obtained by integrating the resulting equations for a number of axisymmetric models. These included (i) axisymmetric decay modes, (ii) a number of mean field a 2- and aw-kinematic models, i.e. solutions of (1) and (2) with V and a given, and (iii) an intermediate dynamical model-Z. The agreement with known or previously calculated results was excellent (except in the case of the dynamical model, where the agreement was only fair). This has encouraged us to look ahead towards the solution of non-axisymmetric models, and as a first step we have attempted to extend the mapping method to asymmetric B, and to examine again dynamos maintained by axisymmetric flows, such as those of G.O. Roberts (Roberts, 1971), of Gubbins (1973) and of Dudley and James (1989) to see if we can recover their resuits. As long as we restrict ourselves to the solution of kinematic problems, it is little less than bizarre to make use of the mapping (3). It can only be seen as a rational enterprise when regarded as a component of a larger program of research, that of solving the MHD equations in the fully three-dimensional case. The precise form of the mapping used in this paper is given in the Appendix. In Paper 1 we explained in detail how such mappings may be implemented. In Section 2 we describe how (1) and (2) were attacked for asymmetric B and steady, axisymmetric and divergenceless V. In Section 3 we present our resuits and make comparisons with earlier work. This is followed by a short concluding section (Section 4).

55

T. Nakajima, P.H. Roberts~Physicsof the Earth and PlanetaryInteriors 91 (1995) 53-61

2. Scalar equations

nents of the source-free potential field in V may be written in the form

We split (1) into cylindrical components but, in place of the cylindrical components of our asymmetric B, we shall employ the fields b~, b, and b,, given by

(4)

+smbz( s, z, t ) I z ] e im¢

Here lq is the unit vector in the direction of increasing coordinate q. The factors of s have been extracted to make b,, b~ and b, even in s and O(1) on the z-axis. According to (2) b,=--

1[

1

0

sm---10S

m

n=m

n=rn

Bn(t)

dP~

.

( n + l ) r "+2 dO-Cm~"

(5)

B~ = Br,

Bo = Bo,

aZ ]

( ~ + imRt )b~ 0bs

7~s(srabs) + v -Lz ovs 1

0bz

(6)

+ A"b" + 2-~-z

[ V~ O ,~ Obz ( -~ + im R ~ ) bz = - R -fi ~s ( S bz ) + Vz-~z

(9)

where Pro(O) is the Legendre function. All components of B must be continuous on S and, by (2), this requires that the solutions of (6) and (7) must be such that Br, B o and OBr/Or on r = 1 are, for some B n, continuous with B,, B o and 01~r/Or calculated from (9):

on

+ aZ

E

OBr (smbs) + s 20b~ ]

This allows us to eliminate b6 from the s and z components of (1), so obtaining

=-R

eim4~,

r

"

/30=-

t)l¢,

B = [sm-lbs(s, z, t)ls+ism-lb4~(s , z,

pm

.

J~r = E n n ( t ) ~

--

Or

on r = --

ar

r= 1.

(10)

The way that (6)-(10) were implemented is described in the Appendix. In brief, a mapping of the form (3) transformed the problem from (s, z)-space to (s', z')-space, where it was rewritten in finite difference form using a second-order accurate central difference scheme, as described in Paper 1, there being N uniform intervals in the s' and z' directions. The diffusion operator (8) appearing in (6) and (7) was split into s' and z' parts, which were integrated in each direction implicitly, except for the cross-derivative terms which, like the motional terms, were added explicitly at each time step; the time interval between steps is denoted by h.

1

- sO---;bs + sos bz ] + Ambz (7) where ~ = V,/s and the diffusion operator, Am, is defined by 02

%.

=

0s---

2m + 1 0 +

- -

s

Os

02 +

-0z - 2

(8)

Solutions to (1) and (2) are required to match, at the boundary, S, of the electrically conducting region V, to a source-free potential field,/~, that exists in the insulating exterior, V. If (r, 0, ~b) are spherical coordinates, then the r and O compo-

3. The decay modes and five dynamo models As V is independent of t, the problem formulated above is an eigenvalue problem for the growth rate, A, which is in general complex: A = ~r + iO, where ~ v~0. We call such modes 'AC'; direct non-oscillatory modes (lq = 0) we call 'DC'. We identify 'the preferred mode(s)', defined as the one(s) for which or--trm~ is largest, by integrating (1) forward in time until the solution approaches the form S ( r , t ) = S ( r ) e ~'

(11)

56

T. Nakajima, P.H. Roberts/Physics of the Earth and Planetary Interiors 91 (1995) 53-61

to any desired accuracy. The marginal state is the smallest value, Rm, of IR I such that trm~x is zero. We first tested our program by using it to determine the decay rates, - ( r , of some of the decay modes of the field, for which V ~ 0 in (1). The decay modes are of two types: poloidal (in

which field on S matches a nonzero source-free potential field in ~') and toroidal (in which B r is zero on S and B - 0). The decay rates are in fact known exactly, as the squares of smallest positive zeros of appropriate spherical Bessel functions. The results of our integrations are shown in Table

(a) 20

,

I

,

i

,

,

I

,

t

,

•

-20

-40

-120

,

,

-80

-40

0

,

,

40

80

,

,

20

R

(b)

,

2, I

,

~

Re(X~,

Fig. 1. Dispersion diagram for Model 2 of Dudley and James, as given by (a) the present calculation, (b) the calculation of Dudley and James (1989). This is an AC dynamo, and the real (or) and imaginary (fl) parts of the complex growth rate ()t) are shown as functions of the Reynolds number R. In (a): e, computed tr; zx, computed I~. Permission to reproduce (b) has been granted by the Royal Society, London and by Dr. R.W. James.

T. Nakajima, P.H. Roberts/Physics of the Earth and Planetary Interiors 91 (1995) 53-61 Table 1 Decay rates, - t r , of the magnetic field for several values of m Poloidal modes

Toroidal modes

Mode

Our value

Exact

Mode

Our value

Exact

S~ S~ S~ S~

9.7942 19.993 32.830 48.149

9.8696 20.191 33.218 48.831

T] 7"22 T~ 7"44

20.086 33.001 48.444 66.329

20.191 33.218 48.831 66.954

1. The truncation level was N = 32 and the time step was in each case h = 0.5 × 10 -4, or smaller. The agreement is as good as may be expected for these values of N and h. We now describe the results we obtained from five different choices of V each of which can, for large enough [R I, self-excite a magnetic field. To make comparison with the work of Dudley and James (1989) easier, we use the toroidal-poloidal representation of V for the first four models:

V=tnl+•s~=V×

[tnl(r)Pnl(O)r ] +

• V X V X [Sn2 (r)Pn2 (O)r]

(12) In general, one would expect to see a sum of such terms each involving a different Legendre polynomial, but in fact each of our four examples has only a single term of each type. In all four cases, A(•, R ) = A ( - • , R) so that we may, without loss of generality, restrict attention to cases where • _> 0. For model 3, we have A(•, R) = A ( • , - R), so that we need consider only R > 0 cases. We studied only the m = 1 asymmetry of fields. We will consider these four models first. (1) Dudley and James (1989), Model 1:

nl=n2=2,

tz(r )=sz(r )=r

sin 7rr (13)

The eigenmodes B are not symmetric with respect to the equatorial plane, but by the correct choice of phase we can make the real part of B r symmetric and the imaginary part antisymmetric. The optimal value of e, the one for which R m is least, appears to be near 0.14, and we therefore set • = 0.14 throughout. For nearly all the range of R shown, the preferred mode is DC (l-I = 0); according to Dudley and James (1989), it is AC

57

(1~ :~ 0) if R < - 105 approximately, but we did not attempt to confirm this. The marginal R m we obtained was approximately 53.1, which is close to the value of approximately 54 quoted by Dudley and James (1989). (2) Dudley and James (1989), model 2:

n I = 1,

n 2 = 2,

t1(r ) = sin zrr,

s2(r ) = r sin ~rr

(14)

The modes are symmetric with respect to the equatorial plane; we examined only those having dipole symmetry. The best value of e is near 0.13, the value we adopted. The preferred mode is AC in the range of R we examined. Fig. 1 shows and 1~ as functions of R and the corresponding results of Dudley and James (1989) for comparison. We found that the marginal R m is approximately 96.8 and for this ~ = - 17.5; Dudley and James (1989) obtained Rm = 95. (3) Dudley and James (1989), model 3:

n I = n 2 = 1,

t t ( r ) = t2(r ) -- sin ~rr

(15)

The modes B are not symmetric with respect to the equatorial plane and are oscillatory. Following Dudley and James (1989), we took ~ = 0.17. We found that the marginal Rm is approximately 151.2 and for this 1~ = - 3 2 . 4 ; Dudley and James (1989) obtained R m -- 155. (4) Gubbins (1973): n l m I'/2 ~-- 2,

t z ( r ) = s z ( r ) = - r sin(27rr)tanh[2~-(1 - r)] (16) As for Model 1, the modes B are not symmetric with respect to the equator, and are steady in the range of R studied. We took • = 0.1 so that we could compare our results with those of Dudley and James (1989). We found that R m = 53.7, which agrees fairly well with the Dudley and James value of 54.7. (5) G.O. Roberts model (Roberts, 1971): the last model (and historically the first) would involve an infinite number of velocity harmonics if represented as a sum of expressions of the form (12). It is better written as

10X

10X

V . . . s. Oz 1 s + s ~ l ~ + -s -~s lz

(17)

T. Nakafima, P.tl. Roberts /Physics of the Earth and Planetary Interiors 91 (1995) 53-61

58

Table 2 The G.O. Roberts dynamo at various truncations (m = 1) N

h

R,,

fl

32 48 64

0.5 x 10 -4 0.22 X 10 4 0.1 X 10 -4 0

96.30 100.48 102.12 104.2

4.80 4.91 4.96 5.02

- -

with 2s X = ~ sin ks sin kz tanh k (1 - r ) , sin kz ~"= tanh kz sin ks tanh k ( 1 - r)

(18)

G.O. Roberts principally studied solutions for k = 8; Dudley and James (1989) also concentrated on this case, as we shall, but they also reported dynamo action for k = 2 and 4. The spatial scale of the velocity (18) is roughly 1 / k and, if k is as large as eight, it is necessary to increase the number, N, of grid points in s' and z' to obtain reasonably well-converged results for B and A. (For the other four models we took N = 32.) We obtained A for three values of N; the time step, h, was chosen so that the Courant-Friedrichs-Levy (CFL) condition was satisfied in each case and so that the parameter hN 2, which enters the discretization of the diffusion operators, was approximately the same for each. The results we obtained for m = 1 are shown in Table 2, together with, in the final row, the Richardson extrapolated values of R m and fl: R,no~ = (642Rm,64 - 482Rm,48)/(642 - 482) and similarly for l~=; these assume second-order accuracy and that accuracy is confirmed when we compare the values of Rm,~ and flo~ shown with the similarly extrapolated values using the N = 48 and N = 32 entries, which a r e R m , ~ = 103.8 and II= = 5.00. Using 15 terms of the expansion of V in a series of terms the form of (12), Dudley and James (1989) reported that Rm = 114; the value obtained by G.O. Roberts and reported by Roberts (1971) is 105.

4. Conclusions This paper has described the present status of our project to integrate the full M H D dynamo

equations by a mapping method designed to make it easy to satisfy both the magnetostrophic equation of motion and the non-local boundary condition imposed on the magnetic field at the surface of the core. In Paper 1, attention was confined to axisymmetric mean field models, both kinematic and dynamic. In this paper we report our initial successes in introducing asymmetry, albeit for a single Fourier mode in longitude and at the kinematic level only. We have shown that the mapping method is as capable as spectral methods in recovering the decay modes of the sphere and also in reliably determining the critical magnetic Reynolds numbers for marginal dynamo action created by five steady axisymmetric flows, namely those of G.O. Roberts, Gubbins and the three simple models of Dudley and James.

Acknowledgments T.N. thanks the Japan Society for the Promotion of Science (JSPS) for a post-doctoral Fellowship under the Society's Japanese Junior Scientists Program. P.R. thanks the National Science Foundation of the USA for support initially from Grant E A R 86-12937 and more recently from Grant E A R 91-04924. Much of the work presented here was carried out during a visit to the University of Sydney, and was supported by the Australian Research Committee under Grant A69330372 and by the University of Sydney. We are grateful to Dr. Ron James and to the Royal Society of Great Britain for permission to reproduce Fig. 8(b) from the paper by Dudley and James (1989). We also thank Dr. James and Maureen Roberts for helpful computing advice.

Appendix A

The map and the surface conditions In the case of Model 2 of Dudley and James (1989), the solutions are symmetric with respect to the equator, so that we only need consider how the triangle (0 ~
T. Nakajima, P.H. Roberts~Physics of the Earth and Planetary Interiors 91 (1995) 53-61

59

Z"

,..

,,+

J J

;i I I I llllllllllllllhm

=S

illl IIIll i i i i i i i i i i i i i I I I l l l l l I I I I I l l l l IIIIIIIIII IIIIIIIIII1 IIIIIIIIIIII IIIIIIIIIIIII IIIIIIIIII1111 I l J l l l l l l l l l l l l I I I I I I I I I I 1 1 1 1 1 1

A ~

;

.V }

Fig. A1. The mapping used in this paper. On the left, the location of the grid points is shown in physical space; the modification of the grid points on the z-axis is evident (see text). On the right, the corresponding grid in mapped space is displayed.

(0 ~ r ~< 1, 0 ~< 0 ~< rr/2). T h e m a p p i n g we used in this p a p e r is shown in Fig. A1. W e may describe it in two steps. First, we take M a p 1 of P a p e r 1:

(11 ,

s = sin ~ r r s

z = ------~_, c o-ss/ 1

,~2-rrs'/] 1,

(A1)

Second, to implement easily the b o u n d a r y conditions on the axis, namely, Obs/Os = abz/~s = 0,

on s = 0,

(A2)

we found it convenient to modify the constant-z' contours defined by (A1) so that they are p e r p e n dicular to the axis, s = 0. T h e grid points on i = 0 given by (A1) were therefore c h a n g e d to have the same z-positions as those on i = 1. This illustrates one of the attractive features of the mapping m e t h o d : provided one makes sure that the m a p chosen is reasonably continuous, one can adapt it at will to simplify the technical tasks or, by increasing the density of grid points in certain areas, to represent better the solutions in regions where the solution changes rapidly. (By changing the location of the i = 0 grid points, we have effectively introduced a discontinuity in (the finite difference representation of) az/as'. W e assessed the effect of this discontinuity by constructing a second p r o g r a m in which the i = 0 points were not relocated and the b o u n d a r y condition on s = 0 was i m p l e m e n t e d in a necessarily m o r e complicated way. W e f o u n d that the n u m e r -

ical results were altered by less than 0.5% in all cases at the N = 32 truncation level.) T h e interval 0 < s' ~< 1 is divided into N equal parts, the grid points being labeled 0 ~< i ~< N. O n e further grid point, i = N + 1, is also introduced. T h e grid points for the interval 0 ~i 0, j >/0 and i + j +.
Fig. A2. The computational triangle, showing the grid in s'z'-space and the three critical diagonals.

T. Nakajima, P.H. Roberts/Physics of the Earth and Planetary Interiors 91 (1995) 53-61

60

outer diagonal, i + j = N + 1. The solutions for the Dudley and James (1989) Models 1 and 3 are not symmetric with respect to the equator, but we can reflect the "triangle" shown in Fig. A2 in the equatorial plane z ' = 0 to provide a mapping for the whole core r = 1, (see Appendix 2 of Paper 1). We shall restrict the discussion below to the case of symmetry. The way that the mapping is implemented has been fully explained in Paper 1, and (except in one respect) we supply here only a brief sketch. The first stage of each step of the time integration is to use the finite difference forms of (6) and (7) to obtain b, and b~, at every grid point of the triangle 1 ~ 1; they are not the vacuum fields, /~s and /~z, at those points. Since the time Paper 1 was written, we have improved the accuracy with which the surface conditions on the magnetic field are applied. It is the objective of this Appendix to explain briefly how this was achieved. In what follows, we represent V in the form (17). It follows from (1) and (4) that O ot(rBr) = Ve(rBr) + R f , on r = 1 (A3) where

f=

1 02X 1 S ~ Br +

OX OB r

s Or O0

After representing ics,

rB r

] im~B r] r~

1

(A4)

and f as sums of harmon-

oe

rB~ = Y'~ Sin(r, t)Pnm(O)e ira6,

(A5)

we obtain from (A3)

osm

~2S2 + 20sm - n(n + 1)S 2 + R f f f

0t

0r 2

Or

on r = 1

(A7)

and the first and third of (10) give 0s m

--+(n+l)S~=0,

on r = l

Or

(A8)

If S is symmetric or antisymmetric with respect to the equatorial plane, every other term in the sums (A5) and (A6) is absent. In practice the sum is stopped after L terms where, to avoid aliasing, we generally chose L to be approximately N/2. We may expand S~ about r = 1 to second order in ( r - 1) 2 and use (A7) and (A8) to provide the first and second derivatives of Sm at r = 1 in terms of Sin(l). In this way, we obtain Sin(l)

Sm(r) -

rn+l ( r - 1) 2

+ - -

S~(1) h- s'm(1) -- R f ~ ] (A9)

where the tilde denotes field values obtained at the previous time step; forward differencing is used for the time derivative, h being the time step. We must now use the known values of b s and b~ on the inner diagonal to evaluate Sn" on that diagonal, and then we must use (A9) to obtain S~(1). It may be noted that r is not constant on the inner diagonal, and we therefore cannot derive S m on the inner diagonal by performing an integration that makes use of the orthogonality of Legendre functions. Instead we proceed as follows. We evaluate at each grid point, (ri, 0i) on the inner diagonal

bi =

b~( r i, Oi) + z i b z ( r i ,

Oi)

riBr( ri, Oi) m Si

n=m c~

f=

E fm(t)Pn~(O)eim6, n=m

(A6)

Sm(ri)

enm(oi)

r7 +m+l

sin m 0i

(AlO)

T. Nakajima, P.H. Roberts~Physicsof the Earth and PlanetaryInteriors 91 (1995) 53-61 (We have suppressed dependences on t and ~b.) Next we compute a modified b i.

(r i --

biM = b i + - 2h

1) 2

× Y ' , [ S y ( 1 ) + R h f n m] e'"(O0 n r/'n sinm

Oi

(All)

We may now write (A9) as

Fi - ~_~A.iSm(1)

- b~ = 0

61

we may use the known sm(1) and (A8) to determine B o on the principal diagonal. Through (A7) and (A8) we may similarly find B o on the outer diagonal. Although this procedure is intricate, most of the necessary matrices and their inverses can be computed before the time iteration is commenced. It is therefore not expensive to apply in terms of computer time.

(A12)

n

where

References

[

Ani

l (ri-1)Z]pm(oi) rin+m+--------"~ 2hrF, sin m 0i

We determine sm(1) by minimizing C = 3" F.2

(A14)

i

with respect to sm(1). This gives

~p [ ~i AniAvi]S'~(1)= ~i AnibM

(A15)

which is a set of simultaneous equations that determine Sin(l). Having solved these, we use (A9) again to determine S m on the outer diagonal. As

10(rS~) OPm e im4~ 0-----~ O0-

Bo= E r n

(m16)

Braginsky, S.I., 1978. Nearly axially symmetric model of the hydron,agnetic dynamo of the Earth. Geomagn. Aeron., 18: 225-281. Dudley, M.L., 1988. A numerical study of kinematic dynamos with stationary flows. Ph.D. Thesis, University of Sydney. Dudley, M.L. and James, R.W., 1989. Time-dependent kinematic dynamos with stationary flows. Proc. R. Soc. London, Ser. A, 425: 407-429. Gubbins, D., 1973. Numerical solutions of the kinematic dynamo problem. Philos. Trans. R. Soc. London, Ser. A, 274: 493-521. Nakajima, T. and Roberts, P.H., 1995. A mapping method for solving dynamo equations. Proc. R. Soc. London, Ser. A, 448, 1-28. Roberts, P.H., 1971. Dynamo theory of geomagnetism. In: A.J. Zmuda (Editor), World Magnetic Survey 1957-1969. IAGA Bull. 28, IUGG Publ., Paris, pp. 123-131.