A STATISTICAL M O D E L 01' R E V E R S A L S IN T i l l ' C R O D Y N A M O 1 by CHRISTIAN R. S O U L L A R D B . S c . T h e University of B r i t i s h C o l u m b i a , 11)98 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L M E N T O F T H E R E Q U I R E M E N T S F O R Till-] D E C R E E O K MASTER. OF SCIENCE T H E FACULTY OF G R A D U A T E STUDIES (Department of E a r t h and Ocean Sciences) We accept, this thesis as c o n f o r m i n g to the required standard THE UNIVERSITY OK BRITISH April COLUMBIA 2001 © C h r i s t i a n R. S e u l l a u b 200-1 Library Authorization In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. OH-/ZL*=>r^HDate (dd/mm/vVvv) *" Name of Author (please print) Title of Thesis: D e 9 r e e : J- C L J ^ ^ j AUi h\Qc Department of -^ArMi J The University of British Columbia Vancouver, BC Canada ^ Year: Q ,^„. f ( ? ^ "2. • O o <<„ ; W „ . r ^ j 4 Abstract I study a simple model of a turbulent dynamo. Using a combination of analytic solutions and scaling arguments I derive a set of governing equations that describe the evolution of the magnetic field. This simplified system predicts behaviour which is qualitatively similar to that seen in the Earth's magnetic field and in numerical simulations. In particular, the model predicts multiple steady-field states, and that in certain cases the larger field states are more stable than their weaker counterparts. This provides a possible explanation for the observed periods of high stability in the geomagnetic field as well as for some of the observed behaviour in numerical simulations. ii Contents Abstract ii Table of Contents iii List of Figures v Acknowledgements vii 1 Introduction 1 2 Background 4 2.1 Governing Equations 4 2.2 Dimensionless Equations 6 2.3 Alfven's Theorem 7 2.4 Magnetic Field Decomposition 7 2.5 Mean Field Theory 8 2.6 The au Dynamo 11 3 Previous W o r k 15 4 Results 21 4.1 Model of the Velocity Field 21 4.2 Coupled Equations 22 4.3 Determination of M(UJ,V , Z 4.3.1 24 B) T Estimate of the Lifetime of Helical Disturbances iii 26 4.3.2 Solution of the Induction Equation 29 4.3.3 Calculation of the Current 37 4.4 Evolution of the Toroidal Field 40 4.5 Properties of the Potential 44 4.6 Reversals 48 5 Applications 61 6 Future Study 66 7 References 68 iv List of Figures 1 Distortion of field lines by helical fluid motion 12 2 Distortion of field lines by differential rotation 14 3 A bistable potential 20 4 The ensemble of upwelling cylinders used in the dynamo model 23 5 Contour integration for inversion of the Laplace Transform 36 6 Integration of the current 41 7 Current as a function of R 42 8 Stationary current as a function of R 42 9 Approximate current as a function of R 43 10 The function f(x) = ci(100/[l + I0x }) 48 11 The function f(x) = (1 + x ) sin(100/[l + 10x ]) 12 The potential (121) with E = 1 0 " , Ro = 10" , R 13 The potential (121) with E = 1 0 " , Ro = 10" , R w = 1, R 14 The potential (121) with E = 10~ , Ro = 10" , R u = 100, R 15 The potential (121) with E = 1 0 " , Ro = 10~ , R = 1, R a = 98 . . . . 51 16 The potential (121) with E = 10~ , Ro = 10" , = 1, R a = 100 . . . . 51 17 The potential (121) with E = 10^ , Ro = 10~ , R = 1, R 18 The potential (121) with E = 10" , Ro = 10~ , = 100, R 19 Figure (18) plotted over a wider range of BT 2 2 49 2 15 5 LJ 15 5 15 5 15 5 w 15 8 5 5 u 8 5 v = l,R = l a = 0.1 . . . a a a 49 . 50 = 0.1 . . . 50 = 1 a = 0.1 52 . . . 52 53 20 Figure (19) plotted over a wider range of B? 53 21 The potential (121) with E = 10" , Ro = 1, IL = 100, R 22 Figure (21) plotted over the range between the two large troughs 54 23 Evolution of the field in the potential of figure (18). B = 1.1, q = .02 . . 58 24 Evolution of the field in the potential of figure (18). B 0 = 0.01, q = .02 . 58 25 Evolution of the field in the potential of figure (18). B 0 = 0.01, q = 0.2 . 59 26 Evolution of the field in the potential of figure (18). B = 1.1, q = 0.2 . . 59 27 Evolution of the field in the potential of figure (18). B = 0.01, q = 0.3 60 28 Another realization of the same situation as figure (25) 60 29 Phase plot of the K u a n g and Bloxham dynamo 64 30 Evolution of the field with non-stationary noise 65 3 Q = .1 54 0 0 0 vi . Acknowledgement s Thanks to my parents and to Michelle. For all your help and support, thank you much. Many thanks are also due to Dr. Bruce Buffett, to whom I owe a great deal. vii 1 Introduction Paleomagnetic data indicates that the Earth has possessed a magnetic field for at least the last 3 billion years (Buffett, 2003). The field is dominantly dipolar in structure (Jacobs, 1994), at least in the region exterior to the Earth, but there are also measureable quadrupole and octopole components. Although a great many investigations, both experimental and theoretical, have been devoted to the subject, the origin of this field remains an outstanding problem in geophysics. Several explanations have been put forward and debated over the years. However, many of them can now be eliminated by simple arguments. For example, the possibility that the field is ferromagnetic in origin is eliminated by the fact that the temperature inside the Earth is well above the Curie temperature of its composing metals. It is also not possible that the present magnetic field is simply a remnant of one that existed when the Earth was formed. To see this, we need only do a simple order of magnitude calculation of the magnetic diffusion time in the Earth's conducting liquid outer core. The diffusion time is approximated by d /w, where d is the radius of the outer core and 2 T] (~ 2 m /s) 2 is the magnetic diffusivity of liquid iron. For the Earth, d /rj 2 is around 20000 years. Thus any field that had no regeneration mechanism would not survive very long. One interesting property of magnetic fields, not shared by other physical fields, is that they can produce their own sources. The source of a magnetic field is current, and moving a conductor in a magnetic field induces current in the conductor. It is this 1 regeneration mechanism, in which a magnetic field essentially amplifies itself, that is the focus of dynamo theory. It is fair to say that most current theoretical studies of the Earth's magnetic field are concerned with field regeneration due to the flow of liquid iron in the outer core. Perhaps the most striking feature of the geomagnetic record is the presence of magnetic reversals; changes in sign of the dipole component of the field. These events are relatively fast, geologically speaking, taking place over only a few thousand years. A typical polarity interval, on the other hand, is around 10 years (Merrill et al., 1996). A l 5 though many interpretations of the data are possible, it appears that the reversal process is one in which successive reversals are statistically independent of one another; the length of one polarity interval, or chron, does not bear on the length of the next. However, the process appears to be nonstationary, in that the average length of chrons seems to vary in time. Nowhere is this variation more apparent than in the so-called superchrons, periods in which the field remained in the same polarity for about 50 million years. Although a great deal of data has been collected on magnetic reversals, their ultimate cause and the factors that contribute to their statistical properties still remain unclear. However, as we shall see, the idea that the magnetic field is ultimately generated by dynamo action in the liquid core is consistent with all of the observed phenomena. For example, the paleomagnetic record indicates that the field has spent roughly the same amount of time in each polarity. This fact follows naturally from an induction description since the governing equations are invariant under changes of sign of the magnetic field. 2 In this work, we present a model of field generation in the core, which, though somewhat simplified, reproduces some of the qualitative features seen in the geomagnetic field. Our approach is to replace the complicated coupled equations of magnetic field induction and fluid flow with a simple system whose behaviour nonetheless captures many of the essential properties of the full problem. Chapter 2 provides an overview of the aspects of dynamo theory that are relevant for our model and Chapter 3 reviews some previous work that has been done along similar lines. Chapter 4 presents our model as well as the calculations that are used to draw predictions from it. Chapter 5 gives some possible applications of the model to explain numerical dynamo results and some observed features of the geomagnetic record. The concluding chapter is devoted to suggestions for future study. 3 2 Background 2.1 Governing Equations The theoretical study of the geodynamo is based on the coupled momentum conservation and induction equations for an incompressible flow. The first describes the influence of various forces acting within the fluid outer core: the Lorentz force and the Coriolis force, along with buoyant and viscous forces. The second describes the generation of magnetic field due to the moving conducting fluid. The Boussinesq approximation to the momentum equation is at + v •V v - — V P+ a G p r - 2f2 x v + — ( V x B ) x B + ^ V v 2 5 o pp, (1) where p is the background density, P is the pressure, a is the thermal expansion coefficient, 0 is the temperature deviation from a static background T ( r ) , p is the magnetic permeability, f i is the angular velocity vector of the outer core, go is the acceleration due to gravity, and v is the viscosity. The evolution of the magnetic field, B , is governed by the induction equation: ? = V x ( v x B ) + ,,V B 2 (2) at where rj is the magnetic diffusivity. The temperature perturbation obeys the heat equation: 89) -77- + v • V 0 + v • V T — K V 9 at 2 4 (3) where K is the thermal diffusivity and the static profile T(r) is maintained by internal heat sources. The fields also satisfy the additional constraints V-v = 0 (4) V-B = 0 (5) Equation (4) results from the incompressibility of the fluid, and (5) from the apparent lack of magnetic monopoles. Equations (l)-(5) form a complete description of a rotating conducting fluid that is influenced by, and generates, a magnetic field. However, the non-linearity of the system precludes possibility of an analytic solution, and in general we can expect the solutions to be chaotic. However, some qualitative observations can be made about the system. For example, if the field B is replaced by —B, the equations are unchanged. This suggests that the fluid does not "know" the sign of the magnetic field. In particular, if there exists a stable configuration of the field of one sign, then there should also exist a stable configuration of the opposite sign. It is not surprising then, considering the chaotic nature of the fluid flow, that systems of this type often exhibit irregular fluctuations in field polarity (Rikitake (1966) provides an example of a mechanical chaotic system that undergoes reversals). W h a t is not clear, however, is the mechanism by which such a chaotic flow produces a coherent field. In order to answer this question, we must resort to simplified models and numerical simulations. The work in this thesis will primarily deal with the so-called kinematic turbulent dynamo, in which the velocity field is a prescribed 5 random function. 2.2 Dimensionless Equations It is often convenient to express the above system in non-dimensional form. This involves a simple rescaling of the variables, though different scalings will result i n different dimensionless parameters. our timescale is A s our unit of length, we choose d, the radius of the outer core, d /rj, 2 the magnetic field scale is (2Qpbirj) ^ , 1 2 and the temperature scale is A T which is based on an average radial gradient of the imposed background over the distance d. The dimensionless momentum and induction equations are then Ro(—+v-Vv) <9B ~8t = — V P — k x v + R a 0 r + (V x B) x B + E V v (6) V x (v x B) + V B (7) 2 2 where af3g d 4 Ra = 0 8 is the Rayleigh number, Ro = 2Qd 2 (9) is the Rossby number, and E v 2Qd 2 is the E k m a n number. 6 (10) 2.3 Alfven's Theorem Consider equation ( 2 ) with 77 = 0 : dB ~dt V x (v x B ) (11) We can use a vector identity on the right hand side to obtain: — + v V B = B • Vv (12) We can recognize the left hand side as the advective derivative of B , describing its rate of change with respect to the moving fluid elements. Equation ( 1 2 ) is the same equation as that satisfied by a vector that permanently connects two fluid elements. That is to say, the magnetic field is "frozen i n " to the fluid. This result is known as Alfven's theorem. It tells us that a magnetic field can be amplified by velocity gradients (i.e. V v ) in the fluid. Of course, in a physical situation rj 7^ 0 , but the theorem still holds approximately. It then provides us with a convenient way to visualize how magnetic field lines will be stretched and distorted by the conducting fluid before the effects of diffusion set in. 2.4 Magnetic Field Decomposition Condition (5) allows to separate the magnetic field into poloidal and toroidal components defined by: Bp = VxVx B = Vx[T(r,^,t)f] T 7 [P(r, 6, cf>, t)t] (13) (14) This decomposition is convenient, not only because it reduces the number of functions to be considered, but also because it allows us to view the regeneration process as an interplay between these two separate fields. 2.5 Mean Field Theory The focus of this work is on a simplified model of a turbulent dynamo, i.e. a turbulent flow that maintains a coherent magnetic field. As such, we will be required to consider the velocity and.magnetic fields as having an essentially random component. Furthermore, we assume that the randomness in both fields occurs on a scale much smaller than macroscopic scales of interest. As such, we write: V = V +v (15) B = B+b (16) where V and B are the average parts and v and b are the fluctuating parts of the fields. Furthermore (v) = 0 (17) (b) = 0 (18) where (•) denotes a spatial average taken over regions much greater than the correlation lengths of the random fields yet still much smaller than the length scale over which the averages vary. Assuming these length scales are well-defined, random fields in spatially 8 disparate regions can be considered different realizations of the same random process. We will be considering a turbulent kinematic dynamo, meaning that the statistical properties of v are assumed given, although we will justify many of these properties by means of approximate solutions of the momentum equation. Substituting (15) and (16) into (2), we can similarly separate the induction equation into mean and fluctuating parts: — ot = (9b — = Vx(VxB) + Vx£ + (19) T]V B 2 V x ( V x b ) + V x ( v x B ) + V x G + ??V b 2 (20) where £ = (v x b) (21) G = vxb-(vxb) (22) Comparing (19) with (2), we can see that the averaged field has an additional source term, V x £, which depends on the statistical properties of the random fields. If we could express this term as a function of B , we would be in a position to solve for the large-scale field. It can be shown (Moffatt, 1978) that under the two-scale hypothesis, this relationship is of the form: dB St = OijBj + frjk-^- (23) where ctij and Pijk are pseudo-tensors and summation is implied over repeated indices. The pseudotensors are functions of the statistical properties of the velocity field. In the 9 case of homogeneous and isotropic turbulence, we have: aij Pijk = ady (24) — flCijk (25) so that £ = a B — 3V x B (26) The source in (19) is then Vx£ = a V x B + /5V B 2 (27) It is clear that the second term simply adds 3 to the magnetic diffusivity, making the total n + 3. This leads to the unsurprising conclusion that one effect of the randomness in the velocity field is to increase the dissipation of the magnetic field. The effect of the first term can be determined by examining the mean electric field given by £ { 0 ) = aB (28) Since £ ^ is just an effective electric field, it will produce a current according to Ohm's law: J = a£^ (29) where a is the conductivity. This relationship leads to: J = aaB Equation (30) (30) indicates a less obvious feature of the turbulent dynamo; the existence of an average current parallel to the field. This is not possible in a laminar theory of 10 field generation since in that context the current is given by V x B , which can never be parallel to B . It is not clear, however, how this current actually comes about in a physical situation. Parker (1970) was the first to show how this can work. He envisioned an ensemble of small-scale helical upwellings that produce local currents with a component parallel to the large-scale field. According to Alven's theorem, the lifting and twisting motion of the fluid produces the loops shown in figure 1 along with an associated current. Averaging over many of these "cyclonic events" as he called them, produces an overall current satisfying (30). A s we will see, helical upwelling is a prominent feature of rotating convection, making this mechanism a strong candidate for field generation in the Earth's core. 2.6 The ato D y n a m o In order to determine how field regeneration occurs, we need a model of how the poloidal and toroidal fields work to reinforce each other. A s mentioned earlier, differential velocities in the conducting medium can lead to field amplification through stretching of field lines. In a rotating fluid, differential rotation of the fluid is one means by which the toroidal field can be amplified by the poloidal field (figure 2). Differential rotation can arise from the geostrophic component of the flow. This type of flow results when the Coriolis and pressure forces balance each other, and is characterised by a velocity field that is constrained to the ^-direction and has no variation with z (assuming cylindrical coordinates centred at the Earth's centre). Such a flow is characterised by a set of nested 11 Figure 1 : The distortion of lines of toroidal field by helical fluid motion. This type of motion can produce small-scale currents that are a) anti-parallel and b) parallel to the large-scale toroidal field. 12 quasi-rigid rotating cylinders, similar to those shown in figure 2. This provides us with an obvious mechanism for describing the regeneration of the toroidal field. However, if the flow is strictly geostrophic then there is no way to generate the poloidal field from the toroidal field. To complete the dynamo cycle, we can appeal to via the mechanism encapsulated in (30), known as the a-effect, for the generation of poloidal field. The resulting cycle is called an au> dynamo. A model of this type is the focus of this study. One final thing to note is that the rotation of helical upwellings and the sense of the differential rotation for the Earth, both of which being directly influenced by the Earth's rotation, imply that the single field loop (figure la) leads to a current that reinforces the present field, and a double loop will create a current that serves to run down the field. Thus if a majority of the small-scale currents are anti-parallel to the toroidal field, the present polarity (whichever of the two it may be) will be reinforced. Conversely, an excess of parallel currents favours the opposite polarity and a reversal may ensue. 13 b) a) Figure 2: Schematic representation of the distortion of a field line due to differentially rotating cylinders. The initially purely radial field shown in a) is distorted to produce a (/•-component, i.e. a poloidal field adds to a toroidal field. 14 3 Previous Work Due to the complexity of the geodynamo problem, exact results are not possible. Even outside of the geophysical context, it is rare to find examples of self-sustaining dynamos whose properties are completely understood. We must resort to numerical solutions or simplified models. It is only recently that numerical simulation of the geodynamo has become realistic (Glatzmaier and Roberts, 1997; Kuang and Bloxham, 1999). Though they are set up slightly differently, both the Glatzmaier and Roberts and the Kuang and Bloxham dynamos produce self-sustaining field generation over many magnetic diffusion times. To date, only the Glatzmaier dynamo has exhibited a magnetic reversal in a long run (Glatzmaier and Roberts, 1995). However, as we shall see, the Kuang and Bloxham dynamo has also reversed but it was at the outset of an aborted run. Neither model has yet been able to produce a solution using Earth-like parameter values. The main problem seems to be the smallness of the E k m a n number (~ 10~ ) (Zhang and Schubert, 2000). This 15 parameter cannot be set to zero due to the singular nature of the equations in this limit; the solutions obtained for E = 0 are fundamentally different from those obtained with a finite E and then applying the limit E —> 0 (Foias, 2001). However, using a very small E leads to instability in the code for which no solution has yet been found. The present work is inspired by the order-of-magnitude analysis by Buffett and Bloxham (2002) of the Kuang and Bloxham numerical dynamo. The authors showed that by 15 using seemingly crude approximations in analyzing the solution, one can gain valuable insight into the behaviour of the system. In fact, they were able to clearly demonstrate the generation of toroidal field via differential rotation, and poloidal field by the alpha effect. Among many other results, they showed that an unambiguous dominant length scale appears when considering the processes of ohmic and viscous dissipation. Their approach was to individually analyze a set of length scales, indexed by the angular parameter m, and to determine the effects of the various processes at each scale. This information is easily extracted from the numerical data. For example, by using order of magnitude estimates for the rates of viscous and ohmic dissipation, they were able to extract a length L m that L over which each of these processes acts at each scale. They found has only a weak dependence on m; in dimensionless units (scaled using d, the m radius of the outer core) L m « 0.05. Upon defining a magnetic Reynolds number for each scale: where v m is a typical velocity, and a scale-dependent diffusivity r] is used because the m numerical code uses this method to model sub-grid processes (Kuang and Bloxham, 1999), they found that the efficiency of dynamo action peaks at one particular length scale; the scale at which R m ~ 1. We will make use of these results in constructing our order-of-magnitude model of magnetic reversals. Many authors have studied simplified models that are designed to capture particular features of the geodynamo. One idea, that inspired the present work, is to exploit the 16 formal similarity between magnetic field generation in a rotating convective system and a system of interacting spins. The idea is that the loops of field that produce a current in a helical upwelling have a formal correspondence with the spin states of the atoms in a ferromagnetic solid. A net magnetization is produced in such a solid when there is an excess of atoms with a particular magnetic orientation (Bergersen and Plischke, 1994; Reif, 1965). Similarly, a net dipole field is created in a rotating fluid when an excess of elemental currents point in a particular direction. The advantage of this approach is that magnetic spin systems, such as the Ising model and its variants (Montroll and Newell, 1953), have been studied extensively in statistical mechanics and many numerical methods and approximation schemes have been devised to deal with their complexities. Mazaud and Laj (1989), who were seemingly among the first authors to exploit this correspondence, devised a simple interacting spin system that reproduced the qualitative features of the geodynamo. Their model consists of an ensemble of spins that can each be in one of two states corresponding to an elemental current parallel or antiparallel to the toroidal field. The sum of these currents determines the strength of the overall dipole field. The spins are subject to a stochastic time evolution, where the state of a given spin in the next time step is random but is influenced by the interaction between its present state and the states of its neighbouring spins. W i t h this model, they were able to reproduce many of the qualitative features of the Earth's magnetic field; magnetic reversals in particular. However, their model was too idealized to allow any direct comparison with the Earth's behaviour or with the results of numerical simulations. In particular, their 17 time evolution was not based on the induction equation, and thus the parameters that appeared were not easy to associate with physical quantities. This work follows the same basic program: we have a fixed number of spins that can either enhance or hinder field growth, whose states evolve stochastically. We extend their work by allowing the loops to be in a continuum of states, and we use an order-of-magnitude approximation to the induction equation and the basic structure of convection to physically motivate the time evolution. This approach, that employs a simplified induction equation for the magnetic field evolution, is quite common in the literature. We focus attention on quasi-kinematic models, where the velocity field is imposed and externally maintained, since it is essentially a model of this type that is considered in the present work. However, for a kinematic aw dynamo, only infinite field growth or inevitable decay to zero is possible unless the parameters are exactly tuned to specific values (Zeldovich et al., 1983; Moffatt, 1978). For the Earth's magnetic field and in numerical simulations, it is quite clear that steady non-zero fields persist over the long term. This equilibration is thus likely a result of the non-linear feedback of the magnetic field on the fluid velocity. However, this interaction can be quite complicated and it is commonly modelled by inserting a term non-linear in B (often referred to as a "quenching" term) into the induction equation that hopefully captures some element of this feedback. For example, Hoyng et. al. (2001a; 2001b) studied a kinematic mean-field dynamo with a simple quenching term that described the effect of the magnetic field on the helicity (a-quenching). However, the form of this term 18 was only vaguely physically motivated and was primarily used "in want of better" (Hoyng et al., 2001a). They also introduced random fluctuations in the helicity to model the effects of turbulence and derived a set of stochastic differential equations for the evolution of the large-scale field. Although the model itself was highly simplified, they retained the full equations for the vector components and solved them numerically. They found that it was useful to analyze the resulting data in terms of a Fokker-Planck equation (Risken, 1996; V a n Kampen, 1976) for the amplitude of the dipole field. In particular, they found that their dipole amplitude satisfied the equation of a heavily damped particle in a bi-stable potential (figure 3), subject to a random driving force. A s a result, the dipole amplitude remained in the vicinity of one of its two stable states with occasional transitions between them; the system exhibits reversals. Following their lead, we also analyze the qualitative behaviour of our model in terms a damped particle in a potential. However, the result of our model is a multi-stable potential in which the dipole field has many stable states, as opposed to the bistable potential of Hoyng et. al. In a series of papers, Blanter et.al. (2002), developed a model of an a dynamo, in 2 which all components of the magnetic field are generated by turbulence. They devised a model of multiscale turbulence, in which helical motions of all scales serve to generate the overall field. The "cyclones", as they called them, are subject to a stochastic evolution, in which they can break apart into smaller cyclones, as in a normal turbulent cascade, and also join together to form larger structures, i.e. an inverse cascade. While their turbulence undergoes a complicated time evolution, they used an order-of-magnitude version of the 19 Figure 3: A bistable potential induction equation, much like ours, for the evolution of the magnetic field components. The generation term in this equation is a function of the numbers of cyclones at the different scales and both positive and negative contributions are possible. Of course, as is always the case in a model of this type, they needed an assumption about the form of the quenching. They introduced a term proportional to B ^ , 1 2 although, as they admit, their choice was arbitrary. W i t h this model, they were able to reproduce reversals as well as some features that also result from our model. In particular they found that the zero field state can serve as a stable transition between two reversing states of finite magnitude. 20 4 Results 4.1 Model of the Velocity Field M y model of the fluid flow is a simple one that resembles real convection. Numerical convection models indicate that the structure of upwelling is essentially columnar (Zhang and Schubert, 2000). F l u i d parcels rise parallel to the axis of rotation with Coriolis effects inducing helical flow. In a turbulent setting convection may be more irregular but we can still assume it has these basic properties. Our model is of an ensemble of upwelling "cells" as pictured in figure 4. The angular velocity of the 2 t h site is denoted Ui and the upwelling (z-component in cylindrical coordinates) velocity given by v . Zi velocity field at the i th site will be assumed to be of the simple form: v = UirH(r 0 - r ) 0 + v H(r Zi where H{x) is the Heaviside step function. a rigid cylinder of radius ro. Vzi = z, v The local 0 - r)z (32) Thus, the upwelling is approximated by In the present work, I will set all velocities equal, i.e. = co. The toroidal field is distorted by the rising and twisting motions of these cells, producing a small-scale current at each site (see figure 1). The largescale current, obtained by performing a spatial average over the cells, gives rise to the poloidal field. However, the poloidal field can either be constructively or destructively influenced by these convection cells, depending on the final orientation of the field loops produced by the fluid. In turn, through the Lorentz force, the magnetic field influences 21 the upwelling. It is important to realize that this effect does not always act in opposition to the flow, and in fact the presence of a magnetic field can be an aid to convection (Kuang and Bloxham, 1999). As such, we will need a model of the interaction of the magnetic field on the upwelling cells. As we will see, to solve the coupled system that can arise in even a simple model is not easy, as these equations will involve an essential non-linearity. However, we cannot simply ignore the feedback of the magnetic field onto the flow because, as will become clear in the next section, we need this feedback in order for stationary states of the field to become possible. 4.2 Coupled Equations In order to study the interplay between the poloidal and toroidal fields, I consider a simplified model that captures the essential features of the aw dynamo. To this end, I introduce the following system, which amounts to an order-of-magnitude approximation to the (dimensionless) toroidal induction equation: dB T - ^ Bp = B^BP-BT (33) = M(u,v ,B )B z T T (34) where R is the Reynolds number associated with differential rotation, and the function w M(u>,v , BT) is a function that will be derived from the mean toroidal current. The z behaviour of the system will be determined by this function and arriving at its form is one goal of this work. Equations (33) and (34) are intended to capture the spirit of the 22 Figure 4: The ensemble of upwelling cylinders used in the dynamo model. 23 acu dynamo; the toroidal field is generated by differential rotation, while the poloidal field is proportional to the mean current generated via the a-effect. It is clear from these equations that if M is independent of the magnetic field, then no steady field other than BT = Bp = 0 is possible unless the constants are tuned to exact values. 4.3 Determination of M(LU,V , BT) Z The function M describes the average current emanating from our convection cells. But in order to derive it, we need a model of how the field interacts with the flow. In particular, we would like to see how an upwelling region, initially of helical structure, evolves under the influence of viscosity and Ohmic diffusion. The governing equation for the fluid is i Ro(— + v • Vv) = - V P + E V v + B - V B - z x v (35) 2 ot We can linearize these equations by looking at small perturbations to a uniform reference state: where v 1 ; v = v + vi (36) B = Bo + Bx (37) P = Po + P i (38) 0 B i , and P i are small perturbations to the reference state. Switching to a comoving frame with velocity v , we can take v = 0 and in our case, B 0 0 0 = Pj-x. Neglecting any terms that are second order in the perturbations, the equations become: R o ^ = - V P i + E V v + Bo • V B i - z x v i 2 x 24 (39) We could then solve this equation coupled with the induction equation: — = V x (vi x B) + r?V B 2 (40) for an initially uniform B and initially helical v . Although it is possible to make a similar first-order approximation to the induction equation and solve the coupled equations, it will not produce the behaviour we wish to study. This is because the frozen-in effect is captured in the first term of (40) and is lost in any approximation. Instead, we need to retain the full induction term V x ( v i x B ) . A s we shall see in the next section, however, the linearized equations are not without value. It is likely that equations (39) and (40) can be solved numerically, but since we can guess, the properties that the solution will have, we can easily replace this system with a readily solvable one. I therefore use a simplified system that has the following properties: 1. A n initially uniform field is deformed by a helical disturbance. 2. The evolution of this disturbance is slowed by viscous and Ohmic diffusion. 3. A given helical disturbance is driven for a limited time. I model convective upwelling by (32), but introduce a finite lifetime, after which the velocity is shut off. This lifetime is calculated in the following section. 25 4.3.1 Estimate of the Lifetime of Helical Disturbances I use (39) to determine the evolution of the fluid. I begin by writing the full induction equation in the form: — = - v • V B + B •V v + V B at 2 (41) and employ the approximation (37) to obtain: ^ = Bo • V + V B 2 V l X (42) I use the linearized equation i n this context because I am only interested i n the decay rate, and not the orientation of the field loops. These equations admit wave-like solutions of the form: (B , ,P ) 1 Vl = (B ,0 ,P )e^ -^ r 1 1 1 1 (43) Substituting equations (43) into (39) and (42) yields the following constraint on the initial conditions: {-iuKo + Ek )v\ 2 = -iAk + i(B„ • k)Bi (44) The solenoidal conditions on the magnetic and velocity fields requires: Bi-k =0 (46) vi-k =0 (47) Inserting (45) into (44) gives -iavx + z x vx = — iPik 26 (48) where a = Rou + iEk - (49) 2 cu + ik 2 v Taking the cross product of k with (48) and using (47), we have —iak x v i = v i ( k • z) (50) iak -v = (k • z)k x v (51) Taking the curl again, we find 2 x x Equations (50) and (51) can be solved for v i . This leads to the dispersion relation: °2 = ( 5 2 ) % A thorough discussion of the properties of these wave solutions is given in Moffatt (1978). Their time dependence is all that interests us here, though, specifically the decay-rate. This information is found by solving (49) for LO. Rearranging (49) we have: Row + (ik Ro + iEk - a)u - iak - Ek - (B k ) 2 2 2 2 4 =0 2 T x (53) The solution is easily found by the quadratic formula: u = -^{CT-Z/C (RO + E) 2Ro 2 - y/(ik Bjo + iEk - a) + 4Ro[zrjr;; + Ek 2 2 2 2 A + {B k ) }} (54) 2 T x Simplifying the discriminant yields: u = -^{a 2Ro - ik {Ro + E ) - J[a - ik (E - Ro)] + 4 R o ( £ f c ) } 2 2 2 2 T 27 x (55) Next, we use the fact that usually R o << 1. Writing (55) as: u = -^-{a -ik (Ro + E) 2 2Ro - [a - z / c ( E - R o ) ] , 2 V + 4 k-^ (E-Ro)P 2 } ( 5 6 ) allows us to make use of the approximation: y/l + x « 1 1 (57) + - X for small x. Applying (57), we find for our two solutions: (B k ) 2 ui U 2 -?Ar = 2 — — T x — (58) a + ik (Ro - E) 2 a ik E (B k f Ro Ro a + ik (Ro-E) 2 T v 1 [ ' x 2 The lifetime of our helical disturbances is derived from the imaginary part of (59), rather than (58), since (59) yields the shorter timescale of the two. After some rearranging and using R o < < 1 we have: ^ ^ _ 2 fc E k (B k ) (Ro-E) 2 2 2 T x „ - — (60) Using (52) we finally have 1 fc£E Ro , "T" fc4(B fcx) (Ro-E) 2 (61) r kl Of course, equation (61) is only the lifetime of a particular fc-mode. In what follows, our discussion will also focus on only one mode (length-scale), which we implicitly assume to be the dominant scale of the system. It is important to note, also, that the analysis leading up to equation (61) is valid only for (Ro — E) > 0, which is the case of interest 28 for the Earth. The next task is to determine the field orientation for a given velocity and lifetime. 4.3.2 Solution of the Induction Equation Equation (40) with velocity fields similar to (32) has been studied by numerous authors (Zeldovich et al., 1983; Parker, 1966; Moffatt, 1978). Since our velocity field has no ^-dependence, the problem is 2-dimensional in the coordinate sense. The r and (f> components of B are then derivable from a potential, A , and the total field can be written z as (Moffatt, 1978): B = -z x (VA ) + BJ Z (62) Equation (40) for these components in the interior (r < r ) becomes: 0 d M dt dBi dt = ~ ^ + vV Ai = dB ~ ^ + ^ (63) 2 1 2 B l (64) To aid finding the exact solution, we assume that the exterior region (r > ro) is nonconducting, i.e. rj oo. The equations in this region are: V A'J = 0 V BJ = 0 2 2 I which are the same equations satisfied by the (65) • finite-conductivity (66) stationary solution in this region. In terms of the potential, A , the r and (p components of the field are: z * r~d6 = 29 ( 6 ? ) B, = (68) The initial field (and the field at infinity) is given by: B 0 = B cos (pr — B sin (p(p 0 0 (69) which implies AJr,(p,0) = B*(r,0,O) = l i m AJr,d>,t) -» B r sin 6 0 (70) r—>oo Um5 (r,^,<) = 0 z (71) The far-field solutions suggest that we write: Ai = <ZB f(r,ty1> (72) I?, = S S o / M ) ^ (73) = ZB f (r,ty* (74) = ZBog (r,ty* (75) 0 7 for r < ro and yl 7 7 Z?f 7 H 0 n for r > r . In all cases S denotes imaginary part. The equations for / and g are: 0 | - ; f <9t Ji dg" dt _ ( r | ; - ( - r <9r ^ <9r ^ r ^ ^ rj ^ r) d A„II dg 11 r dr (^) dr 30 > ' + ( 7 7 ) 2 u - ik9 r U 2 (79) The continuity of the magnetic field leads to the following boundary conditions at the interface: = /"(ro,*) f(r ,t) 0 (80) df". df , 1 dr dr | r o .9 (r ,t) 7 0 = (81) | r o (82) 9 (ro,t) n and integrating equation (64) over radius leads to the jump condition r 0 - r„ dr f(r , = dr (83) t) 0 rjr 0 Equations (70) and (71) lead to the far-field boundary and initial conditions: f (r,0) = Um f (r,t) = r (84) g (r,0) = l i m g (r,t) = 0 (85) n n U n This system can be solved by using the Laplace transform (Dettman, 1984; Boyce and D i P r i m a , 1992) F ( s , r ) = £ [ / ' ( £ , r)], etc. Equations (76) and (79) become: 7 d F* dF 2 . 1 2 ILO + S r - — + r-—+ ar d?F 2 2 ar A 11 r -—-+r-—+ dr dr 1 d G 2 2 dG I 1 r) r -—+r-—+ dr ,d?G" r —r^r dr 1 , s — r 2 7] 7] 3 = 86 r\ r 3 2 S T n I I - l F iui + s o2 , dr 1 J rj dF u 2 r 2 r -l)F 2 r dG , s + r — r - . + ( — r dr r] 7 = / , r7 - l G = 11 2 - ) G H = 87 r\ r 3 rj r , 88 , 3 89 rj These are Bessel equations, whose solutions are Bessel functions (Olver, 1972) of a complex argument. After solving these equations and applying the boundary conditions, the 31 result for F(r, s) is: 2wr /2j ( 1 F (r 7 = 5) ? — " 1 ± v The next step is to invert these transforms. " v " + '' (90) T Though these formulas look complicated they are relatively simple to invert. In fact, all the poles are first-order, meaning that lim ^ (s — s S j i s )F(r,s) n where s is the n pole, is finite and non-zero. This makes the th n inversion of F(r, s) simple: /(r,i) = £ l i m ( 5 S n )F(r,s)e (92) s t S—*Sn n which is derived by applying the residue theorem (Dettman, 1984) to an appropriate integration in the complex plane (figure 5). F(r, s) seems to have poles at 5 = 0, s = — iu, and s = n —o~nV/ o r — ^ where a is the n zero of th n The residue of the pole at s = 0 JQ(Z). yields the stationary solution: A t s = —iu, the residues of the two terms in (90) are equal and opposite, meaning that there is really no pole here. To calculate the residues at s = n Taylor expansion of Jo(iJ(s Jo(i\I + + iu)/r)rn) °r ) li 0 about ~ -^\(s 5 = —c^rj/r^ n x 32 iu, the following s is helpful: + iu) 2 J (a )(s - s ) n — n n (94) then l i m ^ ( s — s ) F ( r , s)e is readily calculable, which, when added to (93) produces 7 s S n st n the solution previously obtained by Parker (1966): f \ r , t) = f st) + E 2 Z4rn , „ / W l Jl(0-nZ:) , - ^ t - i u t ^ (95) 0 The solution for the z component of the field is not as simply derived. After application of the boundary conditions, the Laplace transforms of the solutions are: o s(s + r -,11/ n x 7 5 (96) 0 F (r ,s) J _ Q V G ( r , s) has poles at 0 Q iv r z 2 iLo)2j (iJ^r ) r?2 (s + G (r,.s) = iu) J$(iJ^r ) = 0, Sj = [\ + ^^-~ Q 7 i ^ r 0 —o-jw/r , 2 J Q { i ^ r Q ) - J x ( i ^ T Q (97) ) ) r — iu, and possibly at s = —iu. A s before, the residues of the two terms in (96) are equal and opposite at s = —iu, so there is actually —a rj/rl — iu 2 no pole there. The other poles at Sj = are second-order, which complicates the formulas somewhat. To invert G ( r , s), I employ the convolution theorem for Laplace 7 transforms: £{ f h(t - r)h (r)dT} 2 Jo = H^Hzis) (98) where C{h(t)} = ifi(s) (99) C{h (t)} = H (s) (100) 2 33 2 The most obvious choice is to set H^s) = H {r,s) = 2 F (r 7 0 ; S ) . (101) V JL- V W/2 (5 + w ) i/2 v J o ( ? > yffi ( 1 Q 2 ) r o ) since we have already inver ted F ( r , s). We can then immediately write: 7 0 ,•1/20 1/2 JJ T /<a ) iL e M*) = -—^ oo r , .„ Y + z4r £ — 0 (103) 2 The inversion of (102) is straightforward, the result being: h (r,t) = 2 ro J2 j , \ ^ e (104) Using (103) and (104) in (98), we get the following solution for g (r, £): 7 V^V J o ( e - T f r o ) 1 2__, —:—:— +8 £i J i ( a e n )(^? 0 + i) o 'o + 8 — E - T 7 - v ° e E w< 0 0 The behaviour of the system as t —> oo is given by the stationary solution: 9lt(r) = ^ £ / " (106) The stationary solutions (106) and (93) are identical to the ones obtained for the case of finite conductivity external to the cylinder. This can be easily seen by setting all time 34 derivatives to zero and observing that the equations satisfied by the exterior fields are the same regardless of the conductivity. The full time-dependent solution for g (r, t) is } somewhat complicated looking and it is likely that the sums over m converge to something simpler, but this possibility is not pursued any further. I do not attempt to calculate the external solutions f (r, n of (97) reveals that inverting G (r,s) n t) and g (r, n t). Examination is not a straightforward proposition since finding the residues will require a knowledge of the zeros of a certain combination of Bessel functions, rather than just those of JQ(Z). However, this is probably not much more difficult computationally because the zeros of JQ(Z) must be found numerically anyway. O n the other hand, g (r, n t) is not really important for our purposes because we only want to find the current produced in the rotating region. Actually, one encounters this nontrivial denominator issue when trying to calculate even the interior solution of g(r, t) for the case of finite conductivity everywhere. This was the reason for our electing to follow Parker (1966) and set a = 0 in the exterior region. However, as previously mentioned, the cases of finite and zero exterior conductivity have identical long-term stationary behaviour so we can be reasonably confident that, at least interior to the cylinder, the two solutions are similar. W i t h our solution for the magnetic field in hand, we are now in a position to calculate the current produced by the rising cylinders. 35 Im s ii C • • 2 2 -o r|/r 3 0 • - - -CO 2 2 2 2 -im-G r\/r -ia-o r\/r 2 0 x Q Figure 5: The contour integration used i n the inversion of the Laplace transforms. T h e contour C is completed to the left by an infinite semicircle. The dots indicate the locations of the poles. 36 4.3.3 Calculation of the Current The current is given by the Maxwell equation (Griffiths, 1989): J = - V x B (107) ii We are interested in the current that flows parallel to the initial field, ie, J x = J • x. In terms of our solution in cylindrical coordinates, we obtain: 1 dB , 1 TTT cos c/H dB z Jx = \ir o<p which only depends on B . z z . — sin0 , , (108) \x or Relevant for our purposes is the total current crossing the planes (f> = -rr/2 and (p = 3TT/2 in the rotating region (Figure 6): I* = 2 / Jo Recall that B J (r, -)dr x (109) z inside the cylinder r < ro is described by ^ (r, i ) , so the integral for I 7 z x can be written as I x = —Q% (r ,£) 7 0 + 8t + 8 5 V + — ? — - T —-—1 Recall that <r is the n th n (no) zero of Jo(-z), the zero-order Bessel function of the first kind. The behaviour of the current as a function of u for a fixed t is shown in figure 7. The 37 behaviour is quite clear: the current direction depends on the number of turns the loop . of field is given by the upwelling fluid (see figure 1). The number of turns depends on the velocity of the fluid, as well as the magnitude of the field (the timescale for the upwelling depends on the field). The overall dipole field will then be given by the sum of these currents over the ensemble of upwelling cylinders (figure 4). The stationary (t —> oo) current is shown in figure 8, where the oscillatory behaviour vanishes and only constructive currents are possible. Clearly the multiply twisted state, like that pictured in figure l b , cannot be maintained against the effects of diffusion, at least not in the case of an externally fixed constant velocity. In other words, destructive current states, which ultimately lead to reversals, are short-lived by nature. Another interesting feature of the stationary solution is the decay of the current as co becomes large. This eventual exclusion of field lines from a differentially rotating region is a phenomenon known as flux expulsion (Moffatt, 1978). It is a generic feature of induction equation solutions in which the fluid velocity has isolated regions of vorticity (Moffatt, 1992). In order to make use of (110) in the coming analysis, we will require some sort of approximation. For this, we will use the function: IJt) = W B vr T 8/i z 0 7] . smut . . (111) v ' Comparing figures 7 and 9 reveals the minimal error involved in this approximation. In particular, it is clear that the sinusoidal dependence on to is correct and since it is only 38 this relationship that is important, we will drop the factor of 8 in (111). Furthermore, the minus sign will also be dropped since we will henceforth consider currents that reinforce the present field polarity to be positive. A s mentioned previously, these constructive currents are associated with the single loops of field (figure la) and are always directed antiparallel to the field. The helical motions are short-lived, as also mentioned earlier, and we introduce a finite cutoff time, 5t, which depends on the system parameters as well as the magnetic field itself. This is the time at which we will evaluate (111) to determine the current, since after this time the current is subject only to the effects of diffusion and its orientation is no longer changed by the fluid. W i t h these considerations, our current becomes: I (5t) = —^smcoSt fj, rj (112) x In reality, a helical disturbance like the one described here would obviously not evolve in this way. It would decelerate while rising due to the effects of viscous dissipation and the Lorentz force. Our approximation to this is the present model, where the fluid rises at a constant speed and then simply stops after a finite time has elapsed. A s before, this lifetime is given by (61). Furthermore, the average poloidal (dipole) field is approximately proportional to the average current: BP = fil (113) x This leads to the following expression for the poloidal field in terms of the toroidal field: ,, d2 V To Bp = —B sin T ' — ^ — k 2 E Ro 39 fc4(B kl )2(R (114) where we have identified B 0 with the toroidal field B . T The d /rj 2 factor arises because our original calculation of the timescale was dimensionless (recall that d is the radius of the outer core and rj is the magnetic diffusivity). We now introduce the Reynolds numbers R z = v ro/rj z and R a Although we used = cor^/rj. of Ru, we use ro in the definitions R a as the length scale in the definition d and R because this choice will prove convenient in z future analysis. These definitions lead to the factor of which we will set equal to (d/ro) , 2 k (k would do as well), since ro clearly represents a dominant length scale in the model 2 2 and d/r 0 is a sort of dimensionless wavenumber. For realistic estimates of the physical properties (e.g. low viscosity), we can use the approximation R o — E « Ro. Equation (114) becomes: B P = R B sin Z T (115) :lf R k2E k Tkx)2Ro Ro + fe? Comparison of (115) and (34) indicates that: M(B , T co, v ) z = R z sin 7 fc2E Ro + ( W R Q ( 1 1 6 ) kl We now insert (115) into (33) to obtain the evolution of the toroidal field. 4.4 E v o l u t i o n of the Toroidal F i e l d Using (115) in (33) yields, in dimensionless form: a l Ro" + ~ f c i Equation (117) forms the basis of all subsequent analysis. Although this equation can be integrated numerically, its qualitative behaviour can be elucidated by analogy with a 40 Figure 6: Top view of a rotating cylinder depicting the planes over which the current J x is integrated to find I . x The integration runs over the plane cf) = 2,n/2 from r to 0 and 0 over the plane <> / = 7 r / 2 from 0 to TQ. 41 kit Figure 7: The total current, equation (110), plotted against the angular velocity of the cell. The parameters used are v r /rj z 0 = 1, tij/rl = 0.1, n m a x = 20. Figure 8: The stationary solution for the current. Note that the oscillatory behaviour with to disappears in the long-time limit. 42 /(*) cor. o n 400 Figure 9: A plot of the function f(x) = — {v r /8n) sinut for v r /n z 0 z 0 = 1, tv/rl = 0.1. Comparison with figure 7 indicates that this function provides an acceptable approximation for the current. damped particle in a potential. The equation for such a particle is: dx dx dt ^ ^ dt dV 2 (118) dx 2 where m is the mass, 7 is the damping constant, and V(x) is the potential. If the particle is heavily damped, i.e. m / 7 < < 1, then (118) reduces to: dx ldV 7 dx (119) B y analogy we can write (117) i n the form: dB dV 1 dt (120) dB T which implies that the potential is: -Rryk~, 2 2 0 " o" ^ A"; R 7 4T:> 4 43 v + ^ L B 2 ) R*R 2 E/c Zf + B ) Ro k k y 2 Rk 2 2 2 T 2 sin a Efc Ro 2 x Rofc fcg 4 2 fc (121) r>2 £>r where ci is the cosine integral function define by •°° cos t ci (122) The structure of the potential is examined in the next section. 4.5 Properties of the Potential The potential ( 1 2 1 ) is quite complicated and contains many parameters. The values of the wavenumber k set the scale that is being observed. A more complete description of the system would contain many different scales all influenced by the effect of the same toroidal field BT- This would lead to a potential that described the overall behaviour of the dipole field, and would be some combination of the potentials that describe the different scales. However, as mentioned earlier, Buffett and Bloxham (2002) found that generation of the dipole field was most efficient at one particular length scale, justifying the present approach of focusing on an individual scale. We can make a few general observations about the potential. If any of the parameters describing the dynamo action, R , P ^ , R^, is set to 0 , then the potential takes the form: z V(B ) T = \B% (123) resulting in the expected conclusion that under the sole influence of ohmic diffusion B T = 0 is the only stable state of the system. Another limiting case of interest is E —> 0 . 44 However, it is not possible to set E = 0, since this results in a term of the form sin(l/B ) T in (121) which has infinite oscillations in the vicinity of the origin. Thus the effect of reducing E is to increase the number of stationary states at smaller values of Br- The parameters that describe the convective upwelling, R , and the omega effect, FLj, occur z only in the combination R^Rz, so their effects cannot be separated. This product appears in the sinusoidal term of the potential in (121). The Reynolds number for the alpha effect R a (the rotation velocity of a convective cell) also influences the number of stationary states but it has the additional role of determining the strength of the cosine integral term compared with the sine term. In fact, the oscillatory behaviour of the cosine integral term and the sine term, as functions of BT, is quite similar. The argument for both functions can be approximated by 100/(1 + 10B ) T for plausible values of the physical parameters (see figures 10 and 11). Because of the columnar structure of convection, we can restrict attention to the longest variation wavelength, k — 1. We could also consider k — 0, i.e. no ^-dependence z z at all, but this would eliminate the sinusoidal part of equation (117) and with it the possibility of non-zero stable states. Because the small-scale currents that generate the large-scale field are themselves the result of small scale velocity features, we can focus attention on higher values of the other wave numbers. If we seek to compare our model with the Kuang and Bloxham dynamo, we run into a problem. In their simulation, they use E — Ro = 10~ . Examination of equation 5 (56) indicates that in this case, the magnetic field no longer contributes to the lifetime 45 (imaginary part) of the wave. This contribution from the field is crucial to the present model, since equilibration is not possible without it. To analyze this situation, one would likely be required to study higher-order corrections to the linearized system (39), (42). However, E = Ro = 10~ is hardly Earth-like. In fact for the case of the Earth, the 5 Ekman number is completely negligible with respect to the Rossby number; E = 10~ , 15 Ro = 10~ . Nonetheless, in the absence of anything better, we will be guided in many 5 of our approximations by the numerical results of Kuang and Bloxham. The analysis of Buffett and Bloxham (2002) found that one particular length scale, L = 0.05 in dimensionless units, stands out prominently when considering the processes of ohmic and viscous diffusion. This leads us to set k x = k = 20, and therefore k y 30, in what follows. They also found that induction is most efficient at a magnetic Reynolds number of 1. Their magnetic Reynolds number bears some similarity to both our R a so we set R z to 1 and we keep R a in the range 0.1 — 1. We leave R w and R z as an adjustable parameter, keeping in mind that only the product R^Rz is relevant. Plots of this potential using a range of parameters are shown in figures (12) to (22). Figures (12) to (16) use the Earth-like E = 1 0 ~ for R , a R, produced. z 15 and R o = 1 0 . A range of values - 5 and R^ are considered to give an idea of the kind of behaviour can be In these figures we see that stable configurations exist for extremely small fields (B ~ 10~ ), a direct result of the smallness of the E k m a n number. B y varying the 10 Reynolds numbers for the different parts of the flow, a variety of behaviour is observed. In figure (12), for example, we see that the system has a number of stable states, including 46 zero field. However, the relative stability of these cannot be gleaned directly from the potential by comparing the depths of the wells. It is true that a system in a deeper well will require a larger "kick" to escape, but as we shall see in the next section, when these kicks are provided by stochastic noise, the exact form of this noise will be of paramount importance in judging relative stability. We can observe changes in stability with changes in parameters, however. Reynolds number Figures (13) and (14) illustrate the fact that increasing the (actually, the product RzR^) amplifies the oscillating part of the potential, making the states more stable than they would be at a smaller R^. A l l other things being equal, the states in figure 14 are more stable than those in figure 13. Figures 15 and 16 indicate how a fairly small change in a dynamo parameter in this regime can cause a significant qualitative change in the system, in this case the transition of the zero state from unstable to stable. In figures 17 through 19 I use E = 10~ and Ro = ICT to bring us closer to the values 8 5 adopted in the numerical model. Figure 17 shows the transition from a system dominated by ohmic diffusion to one in which a sustained field becomes possible. Stationary states do exist but they are highly unstable and the system will eventually end up in the zero state. Figures 18 through 20 all show the same potential, one that has many stable states, but plotted over different ranges to emphasize the various features. The overall structure seems to be of two deep wells with a number of smaller ones in between. The two wells, being about 3 orders of magnitude deeper than the smaller ones, are classically much more stable. This leads us to expect that even with our field-dependent noise, these two 47 fix) Figure 1 0 : The function f(x) = ci(100/[l + 10rc ]) 2 states will reverse far less frequently than the others. In the next section, we will see that this expectation proves to be correct. The remaining figures (21 and 22) are included to show that similar behaviour to that described above can also be seen at very different parameters. 4.6 Reversals In the classical analogy we are using here, reversals are of course not possible. Once the system ends up in a stable stationary state, it remains there. However, if we realize that the fluid flow is turbulent and subject to essentially random fluctuations, we can borrow from the theory of stochastic differential equations (Risken, 1996; Van Kampen, 1976). The idea is that the system is still subject to the evolution described by the differential equation but some of the parameters are now declared to be random variables with given statistical properties. In our case, we will consider fluctuations in the upwelling velocity 48 /(*) 1 I1 1 0 .5: ° 4-1C -1 = (1 + x ) sin(100/[l + 2 A 4-id°S° 2-i- 18 2-10" 18 -Vio- 18 D1 1 5 5 other plots of this potential, k = 20, k = 1, k — 30, R z z 49 2 _ .0 Figure 12: The potential (121) with E = l f r , Ro = 10" , x Wx }) 11 T ) |2 1 0 111 V(B ) 42|l0" I 1 i-o.s Figure 11: The function f(x) |1 l|II = 1. = 1, R a = 1. A s with all 50 51 52 54 v defined as follows: z v = v + 5v (t) z z (124) z where v is the mean upwelling velocity and 5v (t) is a Gaussian random process (Grimz z mett and Stirzaker, 1992) with the following properties: {Sv ) z (5v (t)5v (t')) z z = 0 (125) = q 5{t-t!) (126) 2 where q is a parameter describing the strength (variance) of the fluctuations. Condition (125) just means that a fluctuation is as likely above as below the mean, and (126) arises because the correlation time is taken to be extremely small and is effectively zero. These properties are chosen strictly for simplicity and likely do not accurately model the properties of real turbulence (Beck, 2001; Pope, 2000). However, they serve to illustrate the basic idea of this approach. Substitution of (124) into equation (117) yields: — - (K ^sm z Ro 5v r B T — z + z z 0 fc . 0 R " s m k 2 ^y (BTfcx) ~ (R i-)&T k\ + Rk , kHB^HKo) 2 a > v ' where now R = v r /rj. 4 fc2E Ro x ( ) 127 fc| Equation (127) describes the stochastic evolution of the toroidal magnetic field (remember that the poloidal and toroidal fields are proportional). This is a system with "multiplicative noise" (Risken, 1996) which means that the strength of the fluctuations is dependent upon B . T The standard way to proceed is to realize that BT is itself now a random variable and to derive the differential equation solved 55 by its probability distribution function: the Fokker-Planck equation. However, we will content ourselves for now with only a qualitative description of the system, as presented in the previous section. To convince ourselves that our conclusions are valid, we study the related discrete-time system: B n+l = B n + MB {R K n sin z ^BJ fc2E Ro ) 2 ( R o ) - 1) + T B n (128) n kl + where A t is the time increment and, since we are only looking at the qualitative behaviour of a randomly forced particle in our potential, we absorb all the contributions to the noise strength into a single random variable, V. The variance of V is just an adjustable parameter whose value we do not physically justify. I treat V as a Gaussian random variable with the properties analogous to (125) and (126): <r ) = n (r- r ) n where S mn m = 0 (129) s q 2 mn ( 1 3 0 ) is the Kronecker delta. A s an example, I use the parameter values of the potential shown in figure (18) with A t = 0.01 for various noise strengths and initial fields, B . 0 Figures (23) and (24) depict the evolution of the field started at Bo — 1.1, i.e. in the large well, and B 0 = 0.01. The results are unsurprising in that the system simply settles down to the nearest stationary state, with small fluctuations due to the noise. However, no fluctuations are large enough in these runs to allow the system to escape either potential. Increasing the noise strength to q = 0.2 yields figure (25) for B — 0.01; 0 56 the system now reverses between two of the stable states between the large wells. A reversal is the result of a random fluctuation large enough to the send the dipole field into the opposite well. The result for B = 1.1 is shown in figure (26). The system now no longer reverses but is now subject to rather large positive fluctuations. The reason for this is the multiplicative nature of the noise. Because the field is relatively large, the fluctuations are also large. However, when there is a negative fluctuation the noise is then much smaller in the next time step and the system is likely to return to the stationary state. Thus reversals are far more rare in the deep wells than they are in their smaller field counterparts. As expected, increasing the noise strength increases the rate of reversals, as is clear from figure (27). Continuing to increase the noise obscures the stable states and the field is dominated by the noise source. A final interesting case is shown in figure (28). This system was actually started at Bo = 0.01 with q = 0.2 but a different seed integer was used in the random number generator, making this simply another realization of the same random process. There is evidently a series of fluctuations at the beginning of the run that causes the system to end up in one of the deep wells. This transition from the less stable into the more stable configuration is interesting since, as we shall see, it is a phenomenon that has actually been observed in a numerical dynamo simulation. 57 Bj. 10r 2000 4000 6000 8000 10000 Figure 23: The evolution of the dipole field in the potential of figure (18) with B = 1.1, 0 q= .02 Bj, 0.002r 0. 0 . 0 0 1 0 . 0 0 0 5 2000 4000 Figure 24: B 0 6000 8000 = 0.01, q = .02 58 1 0 0 0 0 M 0.004 0.002 Figure 25: B = 0.01 and q = 0.2 0 Bj Figure 26: B 0 = 1.1, # = 0.2. The system is started near a more stable configuration, and there are no reversals though there are large positive fluctuations of the field. 59 0.006 0.004 Figure 27: Same as figure (25) except with q = 0.3. W i t h increasing noise strength, the reversal rate clearly also increases. 2000 4000 6000 8000 10000 " Figure 28: A l l parameters are as in figure (25) but the random number generator was seeded differently. In this realization, the system experiences an initial fluctuation that lands it in the stable state. 60 5 Applications Though at this stage our model is far from being able to yield any testable quantitative predictions, many of the qualitative features we have derived are shared by the Kuang and Bloxham numerical dynamo. Some observations of the Earth's magnetic field may also be amenable to interpretation in terms of our results. The Kuang and Bloxham dynamo has not yet exhibited a spontaneous reversal after the initial transients of a new calculation decay. However, one solution did produce a reversal at the start of a new solution. The initial amplitude of the dipole decreased rapidly and established a new state with a much smaller dipole. The field reversed several times in this low-dipole state before jumping back to a state with a large dipole field, at which point reversals ceased. The end of this run is shown i n figure (29), where the final jump from the unstable "transition state" to the final non-reversing state is depicted. Our model provides a natural interpretation of this behaviour: the system was started with a small enough initial field that the weaker unstable state was briefly observable. Similarly, the absence of reversals in subsequent long runs of the simulation is interpreted, in our model, as owing to the choice of an initial field that resides in a deep well of the potential. The behaviour shown in figure (29) supports this interpretation, since presumably if the system had been started nearer the stable state the transition state would never have been seen. Another possible application for the model is in the understanding of superchrons. 61 As mentioned earlier, these are anomalously long periods of constant polarity. We could try to explain them using the present model as simply being periods in which the system jumps into one of its very stable wells, and the probability of seeing a fluctuation large enough to get back out is extremely low. In fact, our model predicts a record that has normal reversals interspersed with sudden superchrons. However, the actual geomagnetic record does not look like this. In the epochs leading up to the cretaceous superchron, for example, the reversal rate seems to have decreased over time, until it effectively fell to zero for 100 million years. Once reversals started again, their rate gradually increased. It is clear that any model that hopes to capture even the qualitative behaviour of the Earth's magnetic field will have to account for this non-stationarity. We can extend our model to include it in two ways. We can either consider the noise process to be non-stationary, or we can allow for a time-dependent potential. For example, a model in which the noise strength slowly decreases over time would reproduce some features of the geomagnetic field. But in this case, the appearance of a superchron would have nothing to do with the existence of the stable wells; the system could easily be in one of the less stable states but with reduced noise. Figure 30 provides an example of this approach. The noise is a function of time-step defined by q = 0.3i?(sin 0.015n) (chosen n for purely illustrative reasons), where H(x) is the Heaviside step function. W i t h this q , superchrons appear at regular intervals, with stochastic behaviour in between. n It has been suggested that the Earth's superchrons may appear in a similar periodic cycle (Merrill et al., 1996). If we were to consider a time-varying potential, we could explain a 62 superchron as being a period in which the depth of the frequently reversing wells became deeper, and therefore more stable, over time. Of course, a greater physical understanding of the reversal process is necessary to physically motivate either of these possibilities. 63 2 I i i i 0 I 0.2 i i i I 0.4 I i i I 0.6 i i i I 0.8 i , i I 1 Dipole Field Figure 29: A phase plot of the final transition in an aborted run of the Kuang and Bloxham dynamo. The field remained in the stable state until the end of the run. [Courtesy, Dr. Bruce Buffett] 64 Figure 30: Evolution of the field in the potential shown in figure 18 with q now a function of n: q = 0.3iJ(sin 0.015n) where H(x) is the Heaviside step function. Note the two n "superchrons" where reversals have ceased. 65 6 Future Study There are many possible extensions of the present research. For example, many different models of the interaction of the upwelling cells and the toroidal field could be employed besides the crude one used here. Obtaining a more complete description of the interaction of the rising fluid and the toroidal field would bring us closer to being able to compare the results quantitatively with those seen in the numerical models. However, in the end, this description itself will likely require numerical modelling. We also modelled the omega effect in the simplest way possible: by the assumption that differential rotation creates a toroidal field that is proportional to the poloidal field. Although this was sufficient for our purposes, a more complete description will include the effect that an increasing field has on the differential rotation. While this would modify the form of the resulting potential, it does not bear on the prediction of multiple stable states. Another interesting area of study concerns the properties of the noise. A s stated earlier, the assumptions used in this work are not motivated by anything other than the conventional wisdom that an uncorrelated Gaussian distribution is usually reasonable as a first approximation. A better distribution and correlations could be found by studying the fluctuations in the numerical simulations or perhaps even via a more sophisticated physical model than the one used here. Once this has been determined, the next logical step will be to derive the Fokker-Planck equation describing the evolution of the dipole 66 field. From there, the reversal rate, or at least its scaling with parameters, will be calculable in principle. As also mentioned before, this work only looked at one particular scale. In trying to make comparisons with the numerical models and the Earth, we have had to assume that one particular length scale stands out in the field generation process. This may in fact be the case, but it is also very likely that this length scale depends strongly on the input parameters and on the field itself. Thus, our choice of k values was essentially arbitrary. A future endeavor should see this analysis carried out at all scales and a potential derived that is a combination of effects from all of them. We will then be in a position to not only make more quantitative predictions about the behaviour of planetary dynamos, but also to explain in a natural way many of the features seen in the Earth's magnetic field. Finally, the question of how multiple stable states arise in a generic situation could be studied. In the present work, they appeared because of the sine term in our evolution equation. However, this term merely reflected the fact that different fluid velocities resulted in the same overall current. It is quite possible that a general mathematical result underlies this prediction. In particular, it may be possible to prove that for a given model of field generation, if the mapping from the velocity field to large-scale current (and thus large-scale field) is many-to-one, then multiple stable states are a natural consequence. However, the exact statement of such a theorem, let alone its proof, is far from clear at this time. 67 7 References Beck, C , 2001, Measuring nonextensivity parameters in a turbulent couette-taylor flow: Phys. Rev. E , 63, 035303(R). Bergersen, B . , and Plischke, M . , 1994, Equilibrium statistical mechanics, 2nd ed: World Scientific Publishing C o . Pte. L t d . Blanter, E . , Shnirman, M . , and Le Moiiel, J . , 2002, U p and down cascades: Three- dimensional magnetic field model: Phys. Rev. E , 65, no. 6. Boyce, W . E . , and D i P r i m a , R. C , 1992, Elementary differential equations and boundary value problems: John Wiley & Sons. Buffett, B . A . , and Bloxham, J., 2002, Energetics of numerical geodynamo models: Geophys. J . Int., 149, 221-224. Buffett, B . A . , 2003, The thermal state of the earth's core: Science, 299, 1675-1677. Dettman, J . , 1984, Applied complex variables: Dover Publications, Inc. Foias, C , 2001, Navier stokes equations and turbulence: Cambridge University Press. Glatzmaier, G . A . , and Roberts, P. FL, 1995, A three-dimensional self-consistent computer simulation of a geomagnetic field reversal: Nature, 377, 203-209. Glatzmaier, G . , and Roberts, P., 1997, A three-dimensional convective dynamo solution 68 with rotating and finitely conducting inner core and mantle: Phys. Earth Planet. Inter., 91, 63-75. Griffiths, D . , 1989, Introduction to electrodynamics: Prentice-Hall Inc, 2nd edition. Grimmett, G . , and Stirzaker, D . , 1992, Probability and random processes: Oxford University Press. Hoyng, P., Ossendrijver, M . A . J . H . , and Schmitt, D . , 2001a, The geodynamo as a bistable oscillator: Geophys. Astrophys. Fluid Dynamics, 94, 263-314. 2001b, Magnetic field reversals and secular variation in a bistable geodynamo model: Phys. Earth Planet Inter., 125, 119-124. Jacobs, J., 1994, Reversals of the earth's magnetic field: Cambridge University Press, 2nd edition. Kuang, W . , and Bloxham, J . , 1999, Numerical modelling of magnetohydrodynamic convection in a rapidly rotating spherical shell: Weak and strong field dynamo action: J . Comp. Phys., 153, 51-81. Mazaud, A . , and L a j , C , 1989, Simulation of geomagnetic polarity reversals by a model of interacting dipole sources: Earth and Planetary Science Letters, 92, 299-306. Merrill, R. K . , McElhinny, M . W . , and McFadden, P. L . , 1996, The magnetic field of the earth: Academic Press, Inc. 69 Moffatt, H . K . , 1978, Magnetic field generation in electrically conducting fluids: Cambridge University Press. Moffatt, H . K . , 1992, Topological aspects of the dynamics of fluids and plasmas: Kluwer Academic Publishers. Montroll, E . W . , and Newell, G . F . , 1953, O n the theory of the ising model of ferromagnetism: Rev. M o d . Phys., 25, no. 2, 353-389. Olver, F . W . J . , 1972, Bessel functions of integer order in Abramowitz, M . , and Stegun, I., Eds., Handbook of Mathematical Functions:, 355-433. Parker, R. L . , 1966, Reconnexion of lines of force in rotating spheres and cylinders: Proc. Roy. S o c , 17, 60-72. Parker, E . , 1970, The generation of magnetic fields in astrophysical bodies, i . the dynamo equations: A p . J., 162, 665-673. Pope, S., 2000, Turbulent flows: Cambridge University Press. Reif, R., 1965, Fundamentals of statistical and thermal physics, 2nd ed: M c G r a w - H i l l , Inc. Rikitake, R., 1966, Electromagnetism and the earth's interior: Elsevier Publishing Company. Risken, H . , 1996, The fokker-planck equation: Springer-Verlag. 70 Van Kampen, N . G . , 1976, Stochastic differential equations: Physics Reports, 2 4 , no. 3, 171-228. Zeldovich, Y . B . , Ruzmaikin, A . , and Sokoloff, D . , 1983, Magnetic fields in astrophysics: Gordon and Breach, Science Publishers. Zhang, K . , and Schubert, G . , 2000, Magnetohydrodynamics in rapidly rotating spherical systems: Annu. Rev. F l u i d Mech., 3 2 , 409-443. 71
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A statistical model of reversals in the geodynamo
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A statistical model of reversals in the geodynamo Scullard, Christian R. 2004
pdf
Page Metadata
Item Metadata
Title | A statistical model of reversals in the geodynamo |
Creator |
Scullard, Christian R. |
Date Issued | 2004 |
Description | I study a simple model of a turbulent dynamo. Using a combination of analytic solutions and scaling arguments I derive a set of governing equations that describe the evolution of the magnetic field. This simplified system predicts behaviour which is qualitatively similar to that seen in the Earth's magnetic field and in numerical simulations. In particular, the model predicts multiple steady-field states, and that in certain cases the larger field states are more stable than their weaker counterparts. This provides a possible explanation for the observed periods of high stability in the geomagnetic field as well as for some of the observed behaviour in numerical simulations. |
Extent | 2533138 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0103852 |
URI | http://hdl.handle.net/2429/15323 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2004-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2004-0296.pdf [ 2.42MB ]
- Metadata
- JSON: 831-1.0103852.json
- JSON-LD: 831-1.0103852-ld.json
- RDF/XML (Pretty): 831-1.0103852-rdf.xml
- RDF/JSON: 831-1.0103852-rdf.json
- Turtle: 831-1.0103852-turtle.txt
- N-Triples: 831-1.0103852-rdf-ntriples.txt
- Original Record: 831-1.0103852-source.json
- Full Text
- 831-1.0103852-fulltext.txt
- Citation
- 831-1.0103852.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0103852/manifest