Self-consistent dynamo models driven by hydromagnetic instabilities

Physics ofthe Earth and Planetary Interiors, 36 (1984) 78—84 Elsevier Science Publishers By.. Amsterdam — Printed in The Netherlands

78

Self-consistent dynamo models driven by hydromagnetic instabilities D.R. Fearn and M.R.E. Proctor Department ofApplied Mathematics and Theoretical Physics, University ofCambridge, Silver St.. Cambridge. CB3 9E W (Gt. Britain) (Received 20 October, 1983; accepted 12 May, 1984)

Fearn, D.R. and Proctor, M.R.E., 1984. Self-consistent dynamo models driven by hydromagnetic instabilities. Phys. Earth Planet. Inter., 36: 78—84. The dynamics of the Earth’s core are dominated by a balance between Lorentz and Coriolis forces. Previous studies of possible (magnetostrophic) hydromagnetic instabilities in this regime have been confined to geophysically unrealistic flows and fields. In recent papers we have treated rather general fields and flows in a spherical geometry and in a computationally simple plane-layer model. These studies have highlighted the importance of differential rotation in determining the spatial structure of the instability. Here we have proceeded to use these results to construct a self-consistent dynamo model of the geomagnetic field. An iterative procedure is employed in which an a-effect is calculated from the form of the instability and is then used in a mean field dynamo model. The mean zonal field calculated there is then input back into the hydromagnetic stability problem and a new a-effect calculated. The whole procedure is repeated until the input and output zonal fields are the same to some tolerance.

1. Introduction It is now generally agreed that the Earth’s magnetic field is maintained against ohmic decay by the dynamo process (see for example Moffatt, 1978). Although the kinematic aspects of the mechanism are now reasonably well understood, rather less progress has been made on the fully magnetohydrodynamic dynamo. There, the velocities that maintain the magnetic field, instead of being prescribed (as in the kinematic problem), are produced by convective instability, and the problem becomes highly nonlinear. Busse (1975, 1976) has made some analytic progress by supposing that differential rotation in the Earth’s core is weak, so that the poloidal component B ~ and the toroidal component BT of the large scale magnetic field are of the same order of magnitude. The primary large scale force balance is then between the Coriolis force and the pressure gradient (the geostrophic balance). However, estimates of the 0031-9201/84/803.00

© 1984 Elsevjer Science Publishers B.V.

westward drift velocity suggest that the toroidal field in the core may be much larger than the poloidal field (see Moffatt, 1978), and that the former plays a significant role in determining the balance of forces (the so called magnetostrophic balance). The idea that a magnetostrophic balance obtains its given support from consideration of the instabilities which occur in a rotating conducting fluid in the presence of a magnetic field. Soward (1979, 1983) has discussed how in weak-field hydromagnetic dynamos, the relaxation of the geostrophic constraint by the magnetic field leads to a greatly reduced instability threshold, and ultimately to runaway field growth. This process continues until Lorentz and Coriolis forces are cornparable in magnitude; i.e., until the Elsasser number (1 A = B2 /2~ /

is of order unity. (Here B is a-measure of the field

79

strength, fi is the rotation frequency, ~sthe magnetic permeability, p the mean density, and ~ the magnetic diffusivity.) For the Earth, A 0(1) corresponds to a field strength of some 20 Gs so a strong-field model of the geodynamo (BT>> B~) is consistent with a magnetostropic balance (A 0(1)). At larger field strengths (A ~ 0(1)), the magnetic field acts to inhibit convection. Consequently dynamo action must be inhibited and this process provides a mechanism to limit the field strength. Moreover, if A is too large, a magnetically driven instability may set in, which can quickly destroy the field (see Fearn, 1983). Thus it seems most likely that the geodynamo operates in the strong-field regime (BT>> B~)with A 0(1) (where A is based on the toroidal field strength). We therefore seek to develop a model of the geodynamo that has a relatively large toroidal field. It turns out that the success of the scheme elaborated below depends on the magnetic Reynolds number =

=

=

R

m

=

Ur

/i~

(2)

of the differential rotation (where r0 is the core radius and U a measure of the zonal velocity) being large. This is a necessary condition for the existence of a strong zonal field, but it also simplifies and makes plausible other aspects of the analysis, as we shall see below. In view of the relatively close alignment of the geomagnetic dipole with the rotation axis, it is appropriate to decompose velocity and magnetic fields into axisymmetric and non-axisymmetric parts. The axisymmetric magnetic field may be further divided into zonal and meridional parts. The former is produced from the latter by differential rotation, while the meridional field relies for its regeneration on the mean electromotive force (e.m.f.) produced by the non-axisymmetric fields. So much is part of conventional dynamo theory. What we have been able to do in addition is to construct a theory for the mean e.m.f. We suppose that for a given differential rotation and zonál field, the dominant contribution to the nonaxisymmetric field comes from the most unstable convective mode. As Rm is large, it is plausible to calculate the motions in the weakly nonlinear regime, when their spatial forms are essentially those

of the linear eigensolutions. Thus an e.m.f. can be constructed and fed into the mean field dynamo equations to find a steady toroidal field profile that will in general differ from the one previously selected. The process is then iterated to produce a self-consistent zonal (and hence meridional) field. It is not clear whether the converged solution is determined uniquely by the axisymmetric velocity field chosen since we have not yet investigated many cases. Our solution is not a fully dynamic one since we have no theory for the axisymmetric motion (which therefore has had to be prescribed). We have, however, made a considerable advance on linear kinematic theories and we have not had to make use of any “separation of scales” argument. The asymptotic limit Rm ~ leads to a situation which may allow a simplification of the axisymmetric equations, first put forward by Braginsky (1964a, b). Since we are unable to achieve large values of Rm in the computations, we have not carried out a comparison with the asymptotic theory, but in any it does not seem that Braginsky’s theory can case be applied to this case (see section 4). This paper is organised as follows. In section 2 we give the governing equations and show how the iteration scheme is derived. Results are presented in section 3 and in a conclusion we discuss extensions and improvements to the method.

2. Governing equations and the iteration scheme Our model of the Earth’s core consists of a sphere of Boussinesq fluid of radius r0, surrounded by an insulating mantle. (A solid inner core could be incorporated easily but has been ignored for the sake of simplicity.) The sphere rotates at an angular velocity &~ flu, and the fluid has density p. permeability p., thermal diffusivity K, and magnetic diffusivity We suppose (again for simplicity) that thermal convection provides the driving force for the non-axisymmetric instabilities which occur, though the currently favoured mechanism is cornpositional convection (see Loper and Roberts, 1983). We then separate the relative velocity field U and magnetic field B into axisymmetric and non-axisymmetric parts; writing in spherical polars =

~.

80

(r, 9, 4))

sionless form as

U=U(r,9)~+U~(r,9)+ü(r,O,4))

(3)

zXu= —Vp+A[(V XB4))xb

+(v

(4) where the azimuthal average of and UP. write

4)

~.

=

u

and b vanish,

(6) (7)

B=BM[B$+R1Bp+b]

where the tildes have been dropped to denote dimensionless quantities. Averaging the induction equation over 4) gives (see Moffatt, 1978) 1

[B

(v x E)~ s (un. V ) ~

\

—

)(

+s(B~.v + V2 2 /B \ S / s 1 1 1 E 2 4, -(Up. V)(sA) + ~V -~

‘~

(8)

—

\ (9)

—

—

where E(r 9~l= R ‘.

‘

/

(u X B4)) + R~1vx (U~x b) +

V

~12

2o(13)

-

~

V X

.

U=(~j/rtj)[R~Us?p+Up+uJ

=

=

(11)

0. Since v B = 0, we can

(5) We now non-dimensionalise by scaling lengths with r 0, velocities with ~/r0, time with ~ and magnetic fields with the maximum value BM of B. The velocity and magnetic field in (3)—(4) may now be written

~B

ab/at

xb)xB~] +qROr

(2”(u 2ir ~o ‘

—~-—

m

X b~d4) /

(io~ ‘

/

ao/a:=—(u.v)T—Rm(U4s.v)o+qv ~

U

=

0

=

V b .

(14, 15)

where the temperature T= /3r0[T(r, O)+8(r, 0, 4))], and $ is a measure of the maximum temperature gradient. The new dimensionless parameters that appear in (11)—(15) are ~=~~‘b q=-~ (16) 2flic, where a is the coefficient of thermal expansion. The modified Rayleigh number R measures the strength of the buoyancy force and the number q, though small in the core if molecular values are taken for the diffusivities, is taken as equal to unity in the present paper, as is reasonable if ic and ~ are both due largely to turbulent diffusion. Several simplifications have been made in the derivation terms of (11)—(15). Firstly,consistent the inertial viscous are neglected, with and the magnetostrophic force balance. Secondly, the poloidal part of the axisyrnmetric field has been neglected in comparison with the toroidal part, consistent with the scalings of (6)—(7) for Rm>> 1. Finally, nonlinear interactions of the asymmetric variables have been suppressed, except for the calculation of E

is the mean e.m.f. due to the non-axisymmetric field and flow. We have written s = r sin 0 in (8)—(9). For large Rm~the first term on the right hand side of (8) can be neglected; so that only the differential rotation is important in generating toroidal field. In conventional mean field dynamo models E,~, is expressed as approximately equal to aB, where the function a depends on certain quadratic averages of the velocity field. We do not make this “first order smoothing” approximation and instead calculate E~directly from the equations for the non-axisymmetric fields u and b. These may be written (see Fearn and Proctor, 1983a) in dimen-

4,. Although this neglect cannot be rigorously justified, we again appeal to the fact that Rm>> 1. For E4, to be 0(1) (as is necessary for (8)—(9) to have steady solutions which obey the scaling laws), ui. hi R ~ and is therefore small. In a full solution the amplitude of u and b would be determined for R close to its critical value R by an equation of the Landau type relating lu and Ibi to R R~.In this circumstance, the form of the field and flows would be independent of R R~at leading order. Thus if we are not particularly interested in the value of R R ,, we can solve the linearised problem (11)—(15) and then adjust the amplitude of the solution as required. Methods for solving (11)—(15) for arbitrary B, ii, and T have —

—

—

81

been developed by Fearn and Proctor (1983a). Although we have no theory as yet for U, U~, and T (which therefore have to be prescribed), it is possible to use the coupled equations (8)—(15) to obtain a consistent solution with a steady axisymmetric field in which B and B~are not prescribed in advance but emerge as part of the solution. This programme is accomplished by iteration as follows: (1) U. Ui,, and T are prescribed and a zonal field profile B10t is taken as an initial estimate for B. Set N = 0; (2) the eigenvalue problem for marginally stable oscillations of (l1)—(15) is solved to yield R~N) and the corresponding eigenfunctions Ut N) and b1 N); (3) these latter are used to calculate E,~,N)using (10); (4) E,~,N)is written as R>a~(r, 0)B~’>where R~,is an arbitrary constant; (5) steady solutions ~ ~ of (8)—(9) are found with E~,N) replaced by R~7tat”t(r, 9)B1”~~ (so that B~’~>, A~”~ emerge as an eigenfunction with R~,N)as the eigenvalue); and (6) ~ is normalised in the same way as was B~°~ and if ~ is sufficiently close to B(N) then we stop. Otherwise we increment N by I and repeat steps 2—6. Thus, taking a field BIN) at step 2 produces a field B1”~~at step 6. It is clear that if B(N) ~ as N —s oc, the latter field represents a self consistent solution of the equations. In the following section we describe some results for which convergence was found using the above scheme. Convergence was slow and we found much faster convergence when the input field used to calculate B~’~’~ ~ was taken to be the average of ~I) and BIN) (rather than simply B~). An alternative scheme to that described above is to delete step 4 and use E 4, directly in (8)—(9). This has several advantages: it does not involve the assumption that the forcing is simply of the form aB, it avoids having a singular a on the

proceed. However, there is the question of whether a steady solution exists or whether it is the most unstable. This depends on the meridional circulation (Roberts, 1972) and we decided that for the preliminary investigation presented here it would be best to use the (eigenvalue) method which permitted oscillatory solutions. Further work (to be reported elsewhere) solves the inhomogeneous problem. 3. Results In this preliminary paper we are only able to report on the small number of results to hand. A full treatment will be given in a forthcoming paper. It is first necessary to specify U, U~.,and T. The last is chosen consistent with a uniform internal heating rate in the core (see Fearn and Proctor, 1983a), while for the other two, simple fields were chosen from the kinematic dynamo study of Roberts (1972). In particular, we used 2

U= —r(1 U~,= V

x

)

2

.

*

r sin9 v x(S~)

—

(17)

6 2 (18) S = 2.5 M r (1 r) (3 cos 29 + 1) where the constant M = —0.15 was chosen to minimise approximately the critical value of R,, in its first determination. Roberts (1972) noted that the presence of meridional circulation was necessary in order that the mean field dynamo should be steady, and indeed appropriate choice of U~, not only led to low values of R a and relatively simple field structures, but also to steady modes —

being the most unstable, as required.

Figure 1 shows a converged solution for A = 1, = 100. The toroidal field profile is perfectly smooth, but the evolved poloidal field differs from Rm

those of other kinematic models in the associ2 that s2)A) does ated current distribution (cz (v not vanish on the equator (0 = ?T/2). This is because the e.m.f E 5 is an even function of (0 ir/2) —

—

equator, and mathematically this inhomogeneous problem is simpler to solve than the eigenvalue problem outlined above. Separation of the time variable is not possible but since we are interested in a steady solution we can set 8/at = 0 and

*

For historical reasons a different normalisation for U was used for the thermal instability calculations; see Fearn and Proctor. 1983a, eq. (3.11).

82

(a) bb~:uu~,

(b) E~B~(c)BBp

Fig. 1. The structure of the converged solution for the case A =1, Rm = 102, q =1, and azimuthal wavenumber m 2. In (a) the solution of (11)—(15): the real (upper) and imaginary (lower) parts of the r and 0 components of b and u, and the temperature perturbation 0 are illustrated by contour plots of their strength (see Fearn and Proctor, 1983a for further details). In (b), the zonal field profile B used in calculating (a) is shown, together with the resulting e.m.f E~and the èalculated a = E’/B. Since a is singular on the equator (see text), a contour plot is not a very effective way of representing it. To try and improve the picture we used only the interior grid points but the contours remain concentrated close to the equator. In (c), the solution of (8)—(9) using the a shown in (b) is illustrated with contour plots of the toroidal field strength B (left) and the poloidal field Bp (right). There is a noticeable difference between the B plotted in (b) and that in (c). This arises because the B which emerges from the solution of (8)—(9) is represented as a double Chebyshev series before it can be used as an input to (11)—(15). At the present low level of accuracy of solution this inevitably introduces an error but we do not believe it qualitatively affects the results or conclusions of this paper.

and does not vanish in general at the equator, while the quantity aB used in kinematic dynamo models does vanish at 0 = ‘ir/2, as it is the product of two odd functions. (This means that our “a” is singular at the equator, but this fact is adequately taken account of in the boundary conditions.) The wave fields show the beginnings of the characteristic localisation that occurs for large Rm. This critical layer phenomenon, described previously by

Feam and Proctor (1983a, b), adds to the already severe resolution problems and prevents very high values of Rm being studied. In Fig. 2 we show the approach to convergence of the quantities R~, 2, and ~ the frequency of the Rd~= (RmRa)” oscillation. Although these parameters (especially R~)seem sensitive to the form of the field B, Fig. 3, which shows the starting guess Bt0~ and five subsequent iterates, exhibits rapid convergence to

83

(a) 1200

9

9

B

a

B~

9

9

(b)aBB~

B

0

B

B

S

p

(c)aBB~

___

0

5

Fig. 2. The convergence of the critical Rayleigh number R~ (0), the critical dynamo number Rd~((>), and the frequency w~(0) of the asymmetric waves for the case illustrated in Fig. 1. The horizontal axis shows the number N of iterates around the scheme discussed in section 2. The convergence shown was 2)+3B(~”’))/4 as the input field for obtained by A using (B(” iteration N. similar combination was used for the a-effect.

Fig. Some=other converged fields: (a) A Rm =4.1 (b)A 5, R,,,examples = 1 (c) A of = 1, Rm io~.

1,

Figure 4 shows converged solutions for various

the naked eye. The wave fields also converge rapidly. Thus we are convinced that a converged solution has been obtained, although we do not know whether it is unique, or whether it depends on the original guess for B.

2

4

5

10~and five Fig. 3. The zonal subsequent iterates fieldfor profile the case for theshown initial in guess Figs. B 1 and 2. Compared with Fig. 1, some extra contour lines have t,een included to emphasise the convergence.

values of A and Rm. (Although the theory is strictly valid only if RM 1 the computations are simpler for Rm = 0(1) and the qualitative features are similar.) The value of A does not seem to be critical. (Note, though that A cannot be too large otherwise the field is found to be unstable to the field gradient instability: see Fearn, 1983.) At the larger values of Rm the concentration of the instability (and thus E 1,) at a particular location clearly has its effect on the field structure, and the example shown for Rm iO~may not be properly resolved. There is a suggestion that an asymptotic theory may be efficacious here but we have not made any progress on this yet. ~.

4. Conclusions In this paper we have derived equations that describe a self consistent steady dynamo driven by hydromagnetic waves. The toroidal field B is not prescribed, but emerges as part of the final solution. We plan to go on to compute similar solu-

84

tions for a variety of imposed zonal and mendional flow patterns. One obvious shortcoming of the theory is that axisymmetric velocities must be prescribed. To give a theory for them would mean solving the mean momentum equation in the geostrophic limit, a procedure that is known to be fraught with difficulty (Proctor, 1977; lerley, 1982; Soward and Jones, 1983). However, we do plan to address this question in the future. Finally, a comment should be made on the relevance of Braginsky’s (1964a, b) nearly-axisymmetric theory which would appear a priori to be applicable to the present problem in the limit Rm —s 00. In fact, as we have seen, the Braginsky ansatz E,~= aB does not hold near the equator. This is because B vanishes there and the assumptions of the theory cease to apply. Braginsky (1964a) noted that his expansion scheme breaks down at any location where the phase speed of the asymmetric motions matches the zonal flow. Such points always seem to be present when Rm is large for the instability modes we have studied (see Fearn and Proctor, 1983a, b). Braginsky (1984a) gave a modified theory that took account of the singularity but his analysis does not apply to the present problem since it assumes that the poloidal flow has no critical-layer structure, whereas the actual flows which arise out of thermally driven instabilities do have such a structure. Thus, even if Rm could be taken large enough for an asymptotic regime to be obtained we would not necessarily expect agreement with Braginsky’s theory.

Acknowledgement This work is supported by the Science and Engineering Research Council of Great Britain.

References Braginsky. S.l., 1964a. Self excitation of a magnetic field during the motion of a highly conducting fluid. JETP 47: 1084—1095. EngI. trans. Soy. Phys., JETP 20: 726—735 (1965). Braginsky, S.l.. 1964b. Theory of the hydromagnetic dynamo. JETP, 47: 2178—2193. EngI. trans. Soy. Phys., JETP, 20:

1462~1471(1965).

Busse, F.H., 1975. A model of the geodynamo. Geophys. J. R. Astron. Soc., 42: 437—459. Busse, FF1., 1976. Generation of planetary magnetism by convection. Phys. Earth Planet. Inter., 12: 350—358. Fearn, DR., 1983. Hydromagnetic waves in a differentially rotating annulus I. A test of local stability analysis. Geophys. Astrophys. FluidM.R.E., Dynam., 27: 137-162. Fearn, D.R. and Proctor, l983a. Hydromagnetic waves in a differentially rotating sphere. J. Fluid Mech.. 128:

1-20. Fearn, D.R. and Proctor, M.R.E., 1983b. The stabilising role of differential rotation on hydromagnetic waves. J. Fluid Mech., 128: 21—36. lerley, G.R., 1982. Macrodynamics of a-dynamos. Ph.D. thesis, MIT. Loper, D.E. and Roberts, P.H., 1983. Compositional tion and the gravitationally powered dynamo. In: A.M. Soward (Editor), Stellar and Planetary Magnetism. Gordon and Breach, London, pp. 297—327. Moffatt, H.K., 1978. Magnetic Field Generation in Electrically Conducting Fluids. Cambridge University Press. Proctor, M.R.E., 1977. Numerical solution of the nonlinear a-effect dynamo equations. J. Fluid Mech., 80: 769—784. Roberts, P.H., 1972. Kinematic dynamo models. Philos. Trans. R. Soc. London, Ser A: 272: 663—703. Soward, AM., 1979. Convection driven dynamos. Phys. Earth Planet. Inter., 20: 134—151. Soward, AM., 1983. Convection driven dynamos. In: A.M. Soward (Editor), Stellar and Planetary Magnetism. Gordon and Breach, London, pp. 237—244, 2-dynamos and Taylor’s Soward, A.M. and Jones, C.A., 1983. a constraint. Geophys. Astrophys. Fluid Dynam., 27: 87—122.