Predictability of Magnetic Field Reversals

Abstract

Geomagnetic field measurements indicate that at present we may be on the brink of the Earth’s magnetic field reversal, potentially resulting in all the accompanying negative consequences for the mankind. Mathematical modelling is necessary in order to find precursors for reversals and excursions of the magnetic field. With this purpose in mind, following the Podvigina scenario for the emergence of the reversals, we have studied convective flows not far (in the parameter space) from their onset and the onset of magnetic field generation, and found a flow demonstrating reversals of polarity of some harmonics comprising the magnetic field. We discuss a simulated regime featuring patterns of behaviour that apparently indicate future reversals of certain harmonics of the magnetic field. It remains to be seen whether reversal precursors similar to the observed ones exist and might be applicable for the much more complex geomagnetic dynamo.

1. Introduction

The Earth’s magnetic field is predominantly a dipole, which is almost coaxial with the axis of rotation of the Earth (the angle between the axes is approximately ${10}^{\circ }$). It is now accepted that the field is generated by the hydromagnetic dynamo in the melted outer core [1]. The paleomagnetic analyses [2] provided evidence of its existence already as long as 4.2 billion years ago and of a period of efficient generation by convection 4.1–4.0 billion years ago.

Paleomagnetic data also revealed the dynamic nature of the geomagnetic field (see, e.g., the review [3]). The discovery of magnetic anomalies, which have a “zebra-like” stripe structure and are roughly parallel to mid-ocean ridges, and which are interpreted as the record of the Earth’s magnetic field reversals in the ocean floor minerals, turned out to be the conclusive proof leading to full acceptance of the theory of plate tectonics [4]. During the last forty million years, typically once in several hundred thousand–one million years, a reversal of the field happens, in which its polarity changes to the opposite one [1]. The duration of epochs of constant polarity varies randomly. In earlier times, its variability was higher, and the so-called superchrones, i.e., periods of stable polarity that can be up to several dozen million years long, are encountered. A polarity reversal can take up to 10–20 thousand years. As witnessed by the marine paleomagnetic record, reversals and excursions are usually accompanied by collapses of the geomagnetic dipole [5], during which the geomagnetic dipole moment decreases by a factor of 2–5 relative to its temporal mean [6]; in the course of the recent Laschamp excursion, the total field reduced by a factor of 10 [7].

The last reversal started ∼795 thousand years ago and ended 22 thousand years later [8,9]. Although on geological time scales reversals occur almost instantaneously, they are complex phenomena. For instance, the last reversal featured fluctuations of the virtual magnetic pole position, departures from the dipole configuration of the field and variation in the geomagnetic field intensity. The actual switching of polarity within this process could have happened as fast as in just 100 years [10].

The observed behaviour of the geomagnetic field may indicate that a new reversal or excursion is at present in preparation. The field intensity has been decaying by 5% per century since at least 1840, when regular measurements of the geomagnetic field began, or even since 1500 [11], while the quadrupole moment significantly increased since 1800 [1]. The first determination of the North magnetic pole position (defined as a location where the magnetic field is vertical) was carried out in 1831 in the Arctic Canada. Since then, the pole has been moving towards Siberia with the speed of 0–15 km/year, which increased between 1990 and 2005 to the present values of 50–60 km/year. By extrapolation [12], over the decade of the 2020s it will travel 390–660 km further in the direction of Siberia. The trajectory of the magnetic pole is not so straight on longer time scales: over the last 7000 years it moved chaotically around the geographic North pole [13].

No consensus has yet been achieved on how to interpret this behaviour of the magnetic pole. While the collected data indicate that a new reversal may have now started unfolding [14], a reversal or an excursion is not expected earlier than in 500–1000 years [15], and the similarity between the structure of the field at present and prior to the two most recent excursions, the Laschamp (∼41 thousand years ago) and Mono Lake (∼34 thousand years ago) excursions, suggests it will not happen in the relatively near future [16].

The presence of the geomagnetic field is necessary for the existence of life on the surface of the Earth. It prevents the solar wind from blowing off the atmosphere of the Earth and protects it from high-energy cosmic particles. It is unknown whether the magnetic field remains strong enough during the entire reversal to keep its dispersive and protective properties. Accumulated oxygen escape during intervals of an increased reversal rate could have resulted in catastrophic drops of oxygen content in the atmosphere, causing the decrease in marine biodiversity for millions of years at the levels characteristic of mass extinctions [17]. The Laschamp and Mono Lake geomagnetic excursions are synchronous with the extinction of the Neanderthals; the reduction in the shielding efficiency of the magnetosphere could yield a significant depletion of the ozone layer, causing an increase in the UV-B radiation in Europe by at least 15–20% and thus contributing to their demise [18]. By contrast, the numerical estimates of [19] suggest that during the least favourable time interval of the minimum field and domination of the quadrupole magnetic component over the dipole one during a reversal, the galactic cosmic ray flux increases at the ground by factor 3 or less, and the maximum permissible dose of radiation will not be exceeded; consequently, these authors expect a magnetic reversal to be non-fatal for the mankind and nature. A short review [20] summarises the direct influence of a reduction in the magnetic field (such as in buildings with steel reinforcement or during space flights) on living organisms, concluding that while the number of studies is clearly insufficient, some adverse effects are likely, for instance, on the development of embryons, cell signalling (by changing the contents of ions, e.g., calcium) and cardiovascular systems.

The advanced technological development of the mankind makes it vulnerable to the decrease in the geomagnetic field during magnetic pole excursions and dipole reversals. The increasing solar and cosmic radiation during the periods of unfavourable space weather are likely to adversely affect the everyday life of every human, to induce a negative impact on the economic infrastructure (e.g., causing failures of communication systems and electrical power grids), to inflict harm on agriculture and food supply, and to compromise national security. “The combined effects of losses in the satellite, electrical and agricultural infrastructures resulting from a worst-case CME [Sun’s coronal mass ejections] event hitting the Earth during a reversal could spell disaster for the nation’s [US] economy. … It is very likely the real economic loss would be at least several trillion dollars in the first year alone” [21]. The study concludes that “The nation [USA] is ill prepared to handle a disaster on the scale of a CME direct hit during a geomagnetic reversal”. Although some statements in [21] are scientifically wrong (such as “The geodynamo theory will remain unproven until enough research and scientific advancements occur to determine the origin of the magnetic field”), it is difficult to disagree with its recommendation to increase significantly the funding for geomagnetism monitoring and research. (It is notable that these issues are within the scope of the US Air Force research interests).

This implies the practical significance of studying the process of a reversal. Modelling the Earth’s dynamo proved to be a very resource-consuming numerical problem because of the extreme parameter values involved (see, e.g., the review [22] and an account of the latest developments in the geodynamo modelling [23]). Nevertheless, the occurrence of the dipole reversals was encountered in simulations. The first one was found by Glatzmaier and Roberts [24,25], and subsequently several reversals were found in [26,27,28].

A computational model of a dynamo in a rotating spherical shell driven by compositional convection was used to simulate over 6 million years of paleomagnetic time [29]. It mimicked the main features of the Earth’s magnetic field evolution: a predominantly dipole morphology of the external field exhibiting stable polarity epochs; these epochs are interrupted by dipole collapses; some of the collapses give rise to reversals or excursions. Two dipole collapses, one resulting and another one not resulting in a reversal, were analysed in detail. The dipole mode reduction was found to be several times stronger before a reversal than in its absence. Initially, the reversed field is generated inside the fluid layer due to the convective intermittency, and subsequently patches of the reversed magnetic flux appear on the boundary and move towards the poles. The behaviour of the Earth’s magnetic field since 1840 is compatible with this picture, but the high-latitude reversed flux, which is observed in the geomagnetic field, is too weak to forecast an imminent reversal. Other approaches (such as [30,31,32,33]) to the problem of the geomagnetic reversal prediction rely on the fact that a significant reduction in the axial dipole close to a reversal is observed in the paleomagnetic data and do not suggest that a reversal will occur within the next 10 thousand years.

In contrast with other approaches known in the literature, we do not aim to advance as close as possible to the range of parameters describing the conditions in the outer core of the Earth; rather, we use the parameter values that do not demand excessive numerical resources (the CPU time and memory for computing regimes with an adequate spatial resolution and of a duration sufficient for analysis) for obtaining regimes with the magnetic field reversals. In line with this philosophy, we consider the dynamo problem in a plane layer. This emphasises the fundamental properties of solutions not related to a specific geometry (a spherical layer) of the volume of the electrically conducting fluid in the outer core of the Earth and significantly simplifies the computations.

From the dynamical systems point of view, a reversal may be linked to the symmetry reversing the direction of magnetic field and to the presence of heteroclinic cycles connecting weakly unstable magnetohydrodynamic (MHD) regimes with the magnetic field of the opposite polarity. A scenario (i.e., a sequence of bifurcations) taking an MHD dynamical system to a regime involving magnetic field reversals was proposed in [34]. Briefly, on increasing a control parameter (for instance, the magnetic Prandtl number) upon the onset of magnetic field generation, this nonlinear dynamical system follows the standard sequence of bifurcations yielding a steady state, a periodic orbit, a two-frequency quasiperiodic regime and a chaotic attractor. Close (in the parameter space) to the onset of chaos, the chaotic attractors are not wide in the “magnetic direction”, but they grow in size as the control parameter is further increased. Since the system possesses the magnetic field reversal symmetry, initially there exist two separate attractors related by this symmetry (or more, due to other symmetries, which may be present in the system), but at some point they become so large that they touch each other and join into a single one, thus giving rise to a regime with magnetic field reversals.

Dynamo simulations reveal that stronger convection, slower rotation and lower electrical conductivity are beneficial for reversals [35].

Guided by these ideas, we have attempted to find convective MHD regimes involving the magnetic field reversals. The statement of the problem and the numerical methods used are described in the next section. In Section 3 we present a convective hydromagnetic dynamo that demonstrates a heteroclinic behaviour. The reversals of magnetic harmonics that it features are discussed in Section 4 together with the behaviour peculiarities of certain harmonics that seem to be of a prognostic value. Finally, the last section contains some remarks and conclusions.

2. Statement of the Problem

We consider thermal convection in a plane layer of electrically conducting incompressible fluid in the presence of a generated magnetic field. No external magnetic field is imposed. The fluid is rotating about the vertical axis. The Boussinesq approximation is assumed. The flow velocity, $\mathbf{v}$, is governed by the Navier–Stokes equation that involves the Archimedes (buoyancy), Coriolis and Lorentz forces. The magnetic field, $\mathbf{b}$, obeys the magnetic induction equation. The temperature deviation from the steady linear profile at equilibrium, $\theta =T-T_0-(T_1-T_0)x_3$, satisfies the heat transfer equation. In the reference frame co-rotating with the layer of fluid, $0\le x_3\le 1$, the nondimensionalised equations take the following form:

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{v}}{\partial t}=& \mathbf{v}\times (\nabla \times \mathbf{v})+\mathrm{P}{\nabla }^2\mathbf{v}+\mathrm{PRa}\theta {\mathbf{e}}_3+\mathrm{P}\tau \mathbf{v}\times {\mathbf{e}}_3-\nabla p-\mathbf{b}\times (\nabla \times \mathbf{b}),\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{b}}{\partial t}=& \nabla \times (\mathbf{v}\times \mathbf{b})+\frac{\mathrm{P}}{{\mathrm{P}}_m}{\nabla }^2\mathbf{b},\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill \frac{\partial \theta }{\partial t}=& -(\mathbf{v}·\nabla )\theta +v^3+{\nabla }^2\theta .\hfill \end{array}$$

(1c)

Here, p is the modified pressure; $T_0$ and $T_1$ are constant temperatures at the horizontal boundaries $x_3=0$ and $x_3=1$, respectively; and ${\mathbf{e}}_m$ are unit vectors of the employed Cartesian coordinate system; ${\mathbf{e}}_3$ is vertical; the components of vector fields in this basis are enumerated by superscripts. The usual nondimensional parameters are the Rayleigh number, Ra; the Prandtl number, P; the magnetic Prandtl number, ${\mathrm{P}}_m$; and the Taylor number, Ta ($\tau /2=\sqrt{\mathrm{Ta}}/2$ being the nondimensional angular velocity of the fluid rotation about the vertical axis). The standard procedure [36] of nondimensionalisation is assumed (which is summarised, e.g., in [37]); we stress, however, that we consider here abstract convective MHD regimes without bearing in mind any specific physical system. The equations are supplemented by the solenoidality conditions

$$\nabla ·\mathbf{v}=0,\nabla ·\mathbf{b}=0.$$

(1d)

We assume stress-free perfectly electrically conducting horizontal boundaries, where the following conditions are satisfied:

$$\begin{array}{cc}\hfill \frac{\partial v^1}{\partial x_3}{\left.\phantom{\int }\right|}_{x_3=0,1}=& \frac{\partial v^2}{\partial x_3}{\left.\phantom{\int }\right|}_{x_3=0,1}=v^3{\left.\phantom{{|}_{|}}\right|}_{x_3=0,1}=0,\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill \frac{\partial b^1}{\partial x_3}{\left.\phantom{\int }\right|}_{x_3=0,1}=& \frac{\partial b^2}{\partial x_3}{\left.\phantom{\int }\right|}_{x_3=0,1}=b^3{\left.\phantom{{|}_{|}}\right|}_{x_3=0,1}=0,\hfill \end{array}$$

(2b)

$$\begin{array}{cc}\hfill \theta {\left.\phantom{{|}_{|}}\right|}_{x_3=0,1}=& 0.\hfill \end{array}$$

(2c)

In the horizontal directions $x_1$ and $x_2$, the fields are supposed to have the periods $L_1$ and $L_2$, respectively. The boundary conditions (2) imply the time independence of the horizontal components of the mean (over the volume of a periodicity cell) flow and magnetic field; we set them to zero, meaning that no external magnetic field is imposed and the reference frame co-moving with the horizontal mean flow is used.

Solutions to the equations have been computed by the standard pseudospectral methods [38] in the form of truncated Fourier series which satisfy the boundary conditions (2):

$$\begin{array}{cc}\hfill \mathbf{v}(\mathbf{x},t)=& \sum _{n_1=-[N_1/2]}^{[N_1/2]}\sum _{n_2=-[N_2/2]}^{[N_2/2]}\sum _{n_3=0}^{N_3-1}\left(\begin{array}{c}{v}_{\mathbf{n}}^1\left(t\right)cos\left(\pi n_3{x}_3\right)\\ v_{\mathbf{n}}^2\left(t\right)cos\left(\pi n_3{x}_3\right)\\ v_{\mathbf{n}}^3\left(t\right)sin\left(\pi n_3{x}_3\right)\end{array}\right){\mathrm{e}}^{2\pi \mathrm{i}(n_1{x}_1/L_1+n_2{x}_2/L_2)},\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill \mathbf{b}(\mathbf{x},t)=& \sum _{n_1=-[N_1/2]}^{[N_1/2]}\sum _{n_2=-[N_2/2]}^{[N_2/2]}\sum _{n_3=0}^{N_3-1}\left(\begin{array}{c}{b}_{\mathbf{n}}^1\left(t\right)cos\left(\pi n_3{x}_3\right)\\ b_{\mathbf{n}}^2\left(t\right)cos\left(\pi n_3{x}_3\right)\\ b_{\mathbf{n}}^3\left(t\right)sin\left(\pi n_3{x}_3\right)\end{array}\right){\mathrm{e}}^{2\pi \mathrm{i}(n_1{x}_1/L_1+n_2{x}_2/L_2)},\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill \theta (\mathbf{x},t)=& \sum _{n_1=-[N_1/2]}^{[N_1/2]}\sum _{n_2=-[N_2/2]}^{[N_2/2]}\sum _{n_3=0}^{N_3-1}{\theta }_{\mathbf{n}}\left(t\right)sin\left(\pi n_3{x}_3\right){\mathrm{e}}^{2\pi \mathrm{i}(n_1{x}_1/L_1+n_2{x}_2/L_2)}.\hfill \end{array}$$

(3c)

Integration in time has been carried out by the third-order Runge-Kutta method ETD3RK [39] with the exponential time differencing, which was previously successfully used in [37,40,41].

3. The Convective Dynamo Exhibiting a Heteroclinic Behaviour

The dependence of convective MHD regimes on the Taylor number was investigated in [42] for the parameter values P = 1, ${\mathrm{P}}_m$ = 8, Ra = 2300 and $L_1=L_2=\sqrt{8}$ that had been previously considered in the study of convective dynamo [43]. For Ta varying from 0 to 2000, an abundance of regimes was encountered in [42], including MHD steady states (13 branches), time-periodic (8 branches), two-frequency quasiperiodic and chaotic regimes, as well as a finite Feigenbaum period-doubling sequence of bifurcations of a torus producing a chaotic regime and a torus with 1/3 of the cascade frequency. Bifurcations in which these branches emerged and disappeared and the symmetries of the regimes constituting the branches were identified.

The nonlinear stability of the space-periodic steady state for Ta = 675 from the branch ${\mathrm{S}}_8^{\mathrm{R}1}$ [42] was studied in [41] allowing for perturbations, whose period in the direction $x_i$ is $M_i$ times larger for integer $M_i\le 4$. The perturbed regime in the unitary periodicity box (for $M_1=M_2=1$) involves the flow comprised of distorted “wavy” rolls capable of magnetic field generation. The chaotic behaviour in the double periodicity box for $M_1=2,M_2=1$ is associated with a heteroclinic cycle connecting a generating double-frequency quasiperiodic regime and a non-generating steady state.

Following the ideas discussed in the Introduction, we have attempted a (rather limited so far) search for regimes with reversals in the vicinity of the above-mentioned parameter values by varying the Rayleigh number. Simulations have revealed that a similar behaviour displaying reversals is shown by some harmonics of the regime for

$$\mathrm{P}=1,{\mathrm{P}}_m=8,\mathrm{Ra}=2200,\mathrm{Ta}=675,L_1=\sqrt{32},L_2=\sqrt{8},$$

(4)

which we discuss in the sequel. For these parameter values, the time-stepping with the step length of 0.001 has been performed for $N_1=255,N_2=127,N_3=97$ —this resolution was found sufficient in computations [41].

Figure 1 presents the kinetic, $E^{\mathbf{v}}$, and magnetic, $E^{\mathbf{b}}$, energies averaged, like in [41], over the volume of the periodicity cell as functions of time. We observe that the regime consists of what we call segments; each such segment starts when the magnetic energy falls to a local minimum below ${10}^{-6}$. The durations of the segments vary chaotically, although moderately; the mean duration of a segment is 117 time units. All the segments follow the same sequence of the so-called phases, which we are now commenting on.

The details of the evolution of the MHD fields are better seen in Figure 2, a blow-up of Figure 1 for the time interval between 1000 and 1100. A new segment starts at $t\approx 1050$. The magnetic field grows exponentially, while the kinetic energy remains approximately constant (phase i, see the lower plateau of $E^{\mathbf{v}}\approx 106.9$). The end of the plateau is marked by a sharp decrease in $E^{\mathbf{v}}$ to values below 50 followed by an upsurge of the kinetic energy giving rise to a plateau at a higher level $E^{\mathbf{v}}\approx 116.2$ (phase $ii$). The coincidence of the heights of the first and third plateaux of $E^{\mathbf{v}}$ suggests that essentially this is the same flow, but nevertheless the magnetic energy exponentially decays, while the kinetic energy is constant at the first plateau level, but it grows, albeit not exponentially, during the entire time interval spanned by the third plateau. Figure 3, Figure 4 and Figure 5 display isosurfaces of the kinetic energy density (panel (a)) and isolines of the flow velocity components on the horizontal midplane and on vertical cross-sections of the periodicity box (panels (b)–(d)). The plots are the flow snapshots at times $t=1042.996$ (Figure 3, the first plateau, phase $vi$), 1052.996 (Figure 4, the second plateau, phase i) and 1056.996 (Figure 5, the third plateau, phase $ii$). They attest that in phases i, $ii$ and $vi$ the flow consists of rolls.

During the end phase $vi$ of the previous segment (the first plateau in Figure 2), the axes of the rolls become parallel to the diagonal $x_1=x_2$ of the rolls’ square of periodicity (see Figure 3). The horizontal size of each roll is $L_1\sqrt{2}/4=1$. (The length of the bottom side of the largest vertical rectangle normal to the line $x_1=x_2$ inside the original periodicity cell is $L_2\sqrt{2}=4$. It contains two flow periods in this direction, as seen in Figure 3d. In order to increase the resolution of the figure, Figure 3c shows cross-sections of a pair of adjacent rolls rotating in opposite directions in a half of this rectangle.) Surprisingly, this orientation of the rolls during phase $vi$ is the same in all segments. (It must be noted that rolls with the axes parallel and normal to the line $x_1=x_2$ are not equivalent—and they have different stability properties—because rotation by ${90}^{\circ }$ is not a symmetry of the system unless the horizontal periodicity cell of the allowed perturbations is a square.) Independent magnetic field computations for the fixed flow attest that they are capable of kinematic magnetic field generation. The fast decay of the magnetic energy $E^{\mathbf{b}}$ during phase $vi$ is therefore a transitory phenomenon: the plateaux are too short for the magnetic field to restructure itself into the configuration, where the growing magnetic mode(s) becomes dominant so that $E^{\mathbf{b}}$ starts increasing. The minimum of the magnetic energy heralds the beginning of phase i (the lower plateau), the rolls rotate by ${45}^{\circ }$, their axes become parallel to the shortest side of the main rectangle of periodicity (i.e., to the $x_2$ axis) and they become $L_1/3$-periodic (see Figure 4). (Figure 4 shows a half of the largest side of the periodicity cell, i.e., 3/2 of the flow period.) Afterwards, their axes rotate further by ${45}^{\circ }$ about the vertical axis in the same direction (see Figure 5) while the kinetic energy $E^{\mathbf{v}}$ switches from the lower to the upper plateau, marking the transition to phase $ii$. (The sizes of the rectangles shown in Figure 5 and of the rolls are now the same as in Figure 5). The two changes in the flow orientation are illustrated by plots of isolines of the flow velocity components parallel and perpendicular to the roll axes in panels (b)–(d) of Figure 3, Figure 4 and Figure 5.

Separate computations for $\mathbf{b}=0$ have revealed that none of the rolls is stable with respect to convective (amagnetic) perturbations. Examples of a trajectory following a heteroclinic connection in the Rayleigh–Bénard convection in the rotating layer of fluid and visiting rolls of different orientation were (in a purely convective setup) discussed in [44]; changes in the roll orientation were attributed to the action of the Küppers–Lortz instability constrained by the periodicity conditions in the horizontal directions. The same mechanism is at work in the present case; this explains why the horizontal sizes of the rolls remain close (being constrained by the spatial periodicity, the sizes of individual rolls change in phases $vi,i,ii$ from 1 to $\sqrt{8/9}$ and back to 1, respectively).

Upon termination of the three phases $vi$, i and $ii$ associated with the rolls, a quasiperiodic phase $iv$ sets in after a considerably long transitory phase $iii$ (spanning the time interval roughly $1059<t<1074$ in Figure 2). In phase $iv$, the regime involves two basic temporal frequencies, 0.263 and 1.767; the respective periods, 3.802 and 0.566, are well visible in Figure 2. The flow generates the magnetic field, and the magnetic energy grows exponentially, experiencing low-amplitude oscillations. The flow involves structures featuring significantly different spatial scales such as thin filaments and large vortices comparable in size to the smaller period $L_1$ (see Figure 6). This flow is stable with respect to convective (amagnetic) perturbations. When the field becomes relatively strong ($E^{\mathbf{b}}$ exceeds 0.1 at $t\approx 1015$ in Figure 2), its kinematic generation cannot be sustained any longer, and phase v of the saturated field evolution starts. The magnetic energy continues to oscillate near its highest values, surges to the maximum and engages afterwards in an exponential decay. This manifests the beginning of a new segment.

In principle, all connections between the heteroclinic nodes might be designated separate phases of the regime. However, since, unlike phases $iii$ and v, the two transitions between the rolls are abrupt, introducing separate phases for these heteroclinic connections seems unjustified. The exact locations of the boundaries between phases are not well defined. Oscillations with the two basic frequencies of phase $iv$ set in almost immediately upon termination of phase $ii$, although initially they are polluted by noise. While in phase $iii$ the trajectory approaches the invariant quasiperiodic flow in the hydrodynamic subspace of the phase space from a relatively long distance, in the subsequent phase $iv$ this approach is gradually offset by an increasing departure in the magnetic subspace.

Thus, in each segment of the evolution of the dynamo at hand, the trajectory in the phase space follows a heteroclinic cycle, joining four weakly unstable stationary convective (amagnetic) objects, all generating magnetic field kinematically: the steady states comprised of rolls parallel to the $x_2$ axis, or parallel to the diagonals of their horizontal periodicity square, and a two-frequency quasiperiodic flow. These nodes are always visited by the heteroclinic cycle in the same order. The first three segments (up to $t\approx 280$) are significantly further from the heteroclinic cycle than all the subsequent segments: we observe, for instance, that their duration is smaller than that of the subsequent segments, and their magnetic energy $E^{\mathbf{b}}$ falls to much higher (order 0.001) values towards the end of a segment than at the end of the subsequent segments (when the values $E^{\mathbf{b}}<{10}^{-6}$ are admitted).

The nodes of the heteroclinic cycle and the dominant magnetic modes that they can generate kinematically are characterised by different groups of symmetries. Here is the list of the independent symmetries that they possess:

• Diagonal rolls (phase $vi$), whose axes are parallel to the diagonal $x_1=x_2$ of the horizontal periodicity cell. The flow symmetries: $r{\gamma }_{L_2/2}^{x_2},{\gamma }_{\alpha }^{x_1}{\gamma }_{\alpha }^{x_2}$ (any shift along the roll axes), $s_2$. The $L_1/2$-periodicity of the rolls in $x_1$ is a consequence of the first two symmetries: ${\gamma }_{L_1/2}^{x_1}={\gamma }_{L_2}^{x_1}{\gamma }_{L_2}^{x_2}{\left(r{\gamma }_{L_2/2}^{x_2}\right)}^2$. The dominant magnetic mode symmetries: ${\gamma }_{L_1/2}^{x_1}$ (the mode has the periodicities of the generating rolls), $qr{\gamma }_{L_2/2}^{x_2},s_2$ (it has the Bloch structure $exp(\mathrm{i}(x_1+x_2)/L_2)\mathbf{B}(x_1-x_2,x_3)$; the growth rate is 2.15).
• Rolls along the shortest side of the horizontal periodicity box (phase i). The flow symmetries: ${\gamma }_{\alpha }^{x_2},r{\gamma }_{L_1/6}^{x_1},s_2$. The dominant magnetic mode symmetries: $q{\gamma }_{L_2/2}^{x_2}$, $qr{\gamma }_{L_1/6}^{x_1},q{s}_2$. Two dominant modes are associated with a real eigenvalue 1.0 and have the Bloch structure $exp(\pm \mathrm{i}{x}_2/L_2)\mathbf{B}(x_1,x_3)$.
• Diagonal rolls (phase $ii$), whose axes are normal to the diagonal $x_1=x_2$ of the horizontal periodicity cell. The flow symmetries are the same as those of the diagonal rolls of phase $vi$. The dominant magnetic mode symmetries: $q{\gamma }_{L_1/2}^{x_1},qr{\gamma }_{L_2/2}^{x_2},s_2$ (the mode has the Bloch structure $exp(\mathrm{i}(x_1-x_2)/L_2)\mathbf{B}(x_1+x_2,x_3)$); the growth rate is 0.48).
• The time-dependent quasiperiodic regime (phases $iii-v$). The flow and the dominant magnetic mode (its growth rate is 0.14) have the symmetry $r{\gamma }_{L_2/2}^{x_2}$.

We use the same notation for symmetries, as in [41,42]: ${\gamma }_{\alpha }^{x_n}:(x_1,x_2,x_3)\mapsto (x_1+\alpha ,x_2,x_3)$ for $n=1,2$ are shifts in the directions of the horizontal Cartesian axes $x_n$ (above, the shift $\alpha$ is arbitrary); $s_2:(x_1,x_2,x_3)\mapsto (-x_1,-x_2,x_3)$ is the reflection in horizontal planes about the origin (an appropriate horizontal shift of the origin of the coordinate system may be necessary for detecting this symmetry); and $r:(x_1,x_2,x_3)\mapsto (x_1,x_2,1-x_3)$ is the reflection about the horizontal midplane. Finally, the symmetry $q:(\mathbf{v},\theta ,\mathbf{b})\mapsto (\mathbf{v},\theta ,-\mathbf{b})$ reverses the magnetic field.

The symmetries explain why the magnetic field decays during phase $vi$: while the dominant magnetic mode generated in phases $iii$–v by the quasiperiodic flow has the symmetry $r{\gamma }_{L_2/2}^{x_2}$, in the subsequent phase $vi$ the dominant magnetic mode has the opposite symmetry $qr{\gamma }_{L_2/2}^{x_2}$; the field component having this symmetry is significantly inhibited during the preceding phases $iii$–v.

4. Polarity Reversals of Harmonics and Their Precursors

Individual harmonics constituting the magnetic field (more precisely, their coefficients in the expansions (3), but we will not make this distinction in what follows) feature significantly different patterns of behaviour shown in Figure 7: either they “burst” shortly after a segment begins, have a small amplitude and decay fast (for instance, Im$b_{1,1,2}^3$ shown in the upper panel) or are activated near the middle of a segment and oscillate, having a substantial amplitude that varies in time significantly until the field demise at the end of the segment; the oscillations are either almost symmetric about zero (e.g., Im$b_{1,1,1}^3$, the middle panel) or show a major trend towards non-zero values of a specific sign (e.g., Im$b_{4,0,0}^2$, the lower panel). The prevailing sign of the coefficient of a harmonic of the latter type varies between the segments seemingly randomly; it is natural to interpret the change in the sign as the occurrence of a polarity reversal of the respective harmonic.

For this reason, the convective dynamo regime under discussion has attracted our attention as a testbed for the investigation of the problem whether the forecasting of magnetic reversals is possible. While the study [29] suggests a positive answer, it proved impossible to find precursors of reversals in the paleomagnetic data. This question, therefore, remains so far open.

Figure 8 shows the evolution of the coefficients of two harmonics: Re $b_{0,1,0}^1$ (blue) and Re $b_{0,2,0}^1$ (black and red colours). The former one is averaged over intervals of 10 time units to obtain a clear indication of the prevailing sign of this harmonic (this operation is analogous to integrating a harmonic over a segment to find its prevailing sign, like we proceed for the harmonic Re $b_{0,2,0}^1$). It is activated in the beginning of a segment, when the rolls parallel to the $x_2$-axis persist, as a part of the magnetic field generated by these rolls, decreases almost to zero during the transition between phases i and $ii$, again grows exponentially during phase $ii$, but quickly decays and stays at the noise level until the end of the segment. The latter one, Re $b_{0,2,0}^1$, is excited when the transition to the quasiperiodic regime in phase $iii$ begins. This harmonic has a preferred polarity, clearly visible in most segments (e.g., at $760<t<900$, segment № 8, the dominant sign of the coefficient is negative; we assign successive numbers to the segments in the ascending time order), although not in all (e.g., at $500<t<610$, segment № 6, and in the next segment, № 7). For determining the polarity of a harmonic in each segment, we integrate the coefficient for this harmonic over time from the beginning of the segment of interest (for Re $b_{0,2,0}^1$ the integrals are shown in green in Figure 8), and the sign of the integral at the end of the segment is interpreted as the preferred polarity of the harmonic in this segment. A change in polarity at the turn of a segment is regarded as a reversal (in Figure 8, in segments ending with reversals, the plot of Re $b_{0,2,0}^1$ is shown in red).

We do not use the first three segments for testing the regularities potentially prognostic for reversals, since during this initial time interval the regime has not yet saturated: for instance, the minimum magnetic energy does not fall off to the standard values below ${10}^{-6}$ in the beginning of these segments (Figure 1), and the bursts of Re $b_{0,1,0}^1$ are anomalously strong (Figure 8). Out of the total of 11 full segments (the last segment № 15 is incomplete) available for the analysis, 8 segments end up in a reversal of Re $b_{0,2,0}^1$.

The following pattern is notable: a reversal of Re $b_{0,2,0}^1$ occurs at the end of segment № N if and only if in segments № N-1 and № N the spikes of Re $b_{0,1,0}^1$ are of different signs. This regularity is violated only once: in the beginning of segments № 6 and № 7 (the red segments № 3 and № 4 in Figure 8), the spikes of Re $b_{0,1,0}^1$ have the same sign, but a reversal occurs at the end of segment № 7 (this may just reflect the fact that in both segments the integrals are significantly less than their typical values in other segments and thus the polarity of Re $b_{0,2,0}^1$ is not sufficiently well defined).

This pattern of the magnetic field behaviour is apparently associated with the structure of the heteroclinic network (which is complex due to the presence of many symmetries) in the phase space of the dynamical system under consideration. It appears that the trajectory chooses early on the symmetry subspace and the heteroclinic cycle it will follow, and this choice is reflected both in what will be the sign of the harmonic Re $b_{0,1,0}^1$ and which symmetry subspace will be selected in the next segment affecting the prevailing sign of Re $b_{0,2,0}^1$. It remains to be understood, however, why this mechanism operates with a delay of two segments.

The collapse of the dipole component of the magnetic field in the numerical model [29] was preceded by the emergence of magnetic flux spots of reverse polarity on the surface of the Earth and their poleward migration. The mathematical nature of this precursor, formulated in physical terms in [29], may be similar to the one identified here: some time before the reversal, a change in sign occurs in a group of harmonics of wave numbers larger than unity, which an observer perceives as the appearance of a magnetic flux spot or spots of reverse polarity.

Similarly, necessary conditions for reversals of individual magnetic field harmonics have been identified, which are also formulated in terms of the behaviour of certain harmonics prior to the reversals. We discuss now two potentially prognostic regularities of this kind, one involving a magnetic field harmonic only and another one involving flow harmonics.

The temporal evolution of the harmonic Re $b_{0,0,1}^1$ is shown in Figure 9. Four reversals occur during the time span of the 11 segments of the computed time series that are amenable for the analysis; the segments that terminate in a reversal of this harmonic are shown in red. (The prevailing signs of the harmonic in each segment are well defined, and we do not show the plots of the sign-detecting integrals.) Figure 10 presents the plot of the increments per the computational time step h of this harmonic (i.e., the differences Re $\delta b_{0,0,1}^1\left(t\right)=$ Re $(b_{0,0,1}^1(t+h)-b_{0,0,1}^1\left(t\right))$), which coincide with the Euler approximations of the derivative of the harmonic up to the constant factor $1/h$, where h is the time step used in the computations ($h=0.001$). Figure 11, a blow-up of the plot in Figure 10, displays this quantity in detail in segments № № 8–11. We use in Figure 10 and Figure 11 the same colour coding as in Figure 9 to emphasise the segments that are followed by a harmonic polarity reversal. We find that a segment ends up in a reversal only if it involves an interval of oscillations where the envelope curve of the increments has a smooth parabolic shape. Two such intervals are seen in Figure 11, roughly $866<t<900$ and $1155<t<1170$, preceding the reversals in which the two “red” segments terminate. Similar, albeit shorter structures are also present in the two remaining “red” segments (the structure in the last such segment № 13 is slightly less recognisable). The presence of this prognostic pattern is not sufficient for a reversal to happen: similar parabolic shapes are seen in segments № № 6, 12 and 14, which finish without reversals.

We illustrate the third characteristic pattern of a possible prognostic value by Figure 12 showing the temporal behaviour of the flow harmonics Re $v_{0,0,1}^1$ and Re $v_{0,0,1}^2$. Like the prognostic pattern just discussed, it is a necessary but not sufficient condition for a reversal of the magnetic field harmonic Re $b_{0,0,1}^1$ (displayed in Figure 9). We observe that a reversal of Re $b_{0,0,1}^1$ happens at the end of segment № N (such segments are again shown in red in Figure 12) only if the flow harmonics Re $v_{0,0,1}^m$ have the opposite polarity in segments № N-1 and № N-2, and this happens for both $m=1$ and 2. (The polarity of the flow harmonic is perfectly visible for Re $v_{0,0,1}^1$; for Re $v_{0,0,1}^2$, we show the time integral from the beginning of the segment as an indicator of the polarity).

Again, both necessary conditions can be tentatively interpreted as structural features of the dynamical system. It appears that a reversal of Re $b_{0,0,1}^1$ requires that the trajectory passes along a heteroclinic connection a region of the phase space, where the evolution of the system is accompanied by the formation of a structure with the characteristic parabola-shaped envelope of oscillations of the increments Re $\delta b_{0,0,1}^1\left(t\right)$, and follows in advance two heteroclinic cycles with the opposite polarity of both flow harmonics Re $v_{0,0,1}^1$ and Re $v_{0,0,1}^2$.

5. Remarks and Conclusions

We have studied a regime of the magnetic dynamo sustained by thermal convection in a rotating layer of electrically conducting fluid. For the specific choice (4) of the governing parameters, we have obtained a regime of the dynamo, whose evolution is a succession of segments of randomly varying duration. In each segment, the simulated trajectory in the phase space visits four weakly unstable amagnetic convective regimes involving flows that are steady rolls of three different orientations of the axes and a quasiperiodic flow possessing two incommensurate temporal frequencies. All these flows are capable of kinematic magnetic field generation, the rolls are unstable with respect to convective (amagnetic) perturbations (the Küppers–Lortz instability controlling the lifespan of the flows) and the quasiperiodic flow is stable.

The harmonics constituting the magnetic field exhibit three main distinct patterns of behaviour. The prevailing, over a segment duration, polarity of some of the harmonics varies between the segments seemingly randomly. We have focused on monitoring the prevailing signs of the coefficients for individual basis functions in the expansion (3b) regarding the sign change upon the start of the next segment as a polarity reversal of the respective harmonic.

Certain parallels may be drawn between the behaviour of harmonics in the convective dynamo under consideration and those constituting the geomagnetic field. The behaviour of the main harmonic of the Earth’s magnetic field, a dipole, is not fully correlated with the behaviour of other constituent harmonics. The non-dipole components are much weaker than the dipole and evolve more rapidly (on the time scale of decades rather than of centuries for the dipole) [1]. Partially, this can be explained by the location of observers of the geomagnetic field on or over the surface of the Earth, i.e., at a distance of about half a radius of the planet from the core–mantle boundary, inside which the dynamo acts in the melted outer core: Since the mantle is largely dielectric, the field outside the core can be approximated as a gradient of a harmonic potential (see, for instance, a discussion in [45]). When the potential is expanded in spherical harmonics, the components proportional to harmonics of degree n decay at distance r from the centre as $r^{-n-1}$, implying that all non-dipole components decay faster than the dipole ones. The structure of the field on the core–mantle boundary differs significantly from the simple dipolar one (see, e.g., [16]) due to the important role the higher harmonics play. This suggests that forecasting reversals of any given individual harmonic is the underlying general problem to be addressed when studying the problem of the Earth’s magnetic field reversal. Furthermore, the geomagnetic field may consist of harmonics of different types of behaviour: like in the convective dynamo under consideration, some harmonics may experience oscillations that are symmetric about zero, and thus the symmetry reversing the magnetic field, $\mathbf{b}\to -\mathbf{b}$, which plays the crucial role in the Podvigina scenario [34] for the occurrence of reversals, does not necessarily affect the behaviour of all constituent harmonics.

We have used the dynamo regime under consideration to investigate the predictability of reversals of the prevailing polarity for some harmonics of the magnetic field. We have presented three regularities linking the reversals of some harmonics to the behaviour of other quantities (such as other magnetic harmonics, their increments, or the flow harmonics), which represent necessary conditions for the reversals (one such precursor appears to provide also a sufficient condition, but it is not faultless). It must be noted that while we use the examples of two specific harmonics, Re $b_{0,2,0}^1$ and Re $b_{0,0,1}^1$, to discuss precursors for reversals of their prevailing signs, we have identified many similar pairs of reversing and forecasting harmonics. The observed patterns demonstrating a potential prognostic value have a natural hypothetical interpretation in terms of possible specific features of the underlying dynamical system.

Thus, we have resolved the fundamental question of the predictability of reversals: Our results demonstrate that the preparation of polarity reversals of some harmonics is a long process, which can extend backwards well beyond the previous reversal. The reversals can be associated with significantly earlier reversals of other hydrodynamic or magnetic harmonics, or with a characteristic behaviour of some other quantities, and thus can be predicted.

However, our convective dynamo is distinct from the geodynamo: the parameter values are significantly different, as well as the geometries of the fluid reservoirs. We have explored the problem of reversal predictability in the abstract statement. The structure of the respective dynamical system in the phase space controls not just the occurrence of polarity reversals of some harmonics but the collective behaviour of all of them; our results indicate that the trajectories close to the attractor of the convective dynamo equations may show a peculiar behaviour of harmonics well ahead of the reversals of other harmonics, and these peculiarities can therefore forecast well in advance of the event that a reversal will follow. While in this formulation our results appear universal, by their nature they are not intended for a detailed comparison with other predictions for reversals of the Earth’s magnetic dipole—such comparisons are unlikely to reveal any similarities beyond the most general ones. For instance, a significant reduction in the dipole mode accompanying the Earth’s magnetic reversals (see [5,6,29,30,31,32,33]) is not mirrored in our dynamo. Unlike in [30], our potential precursors are not probabilistic.

Our study gives rise to a large number of further fundamental questions: The patterns that we have found, will they remain precursors when tested for a longer evolution of the dynamo or if the evolution is simulated for other values of the governing parameters (including when the parameter values are significantly modified so that no underlying heteroclinic connection exists that is responsible for the natural division of the evolution history into segments like in the present case)? If a large set of similar patterns predicting reversals is assembled, would it be possible to combine them in groups of patterns used simultaneously (using the ideas of pattern recognition, for instance) for minimising the prediction errors? It is of interest to study in detail which properties of the dynamical system representing the considered convective dynamo are responsible for the identified regularities that are apparently prognostic. How can all these ideas and results be transferred to other geometries, including for application to the problem of prediction of reversals of the dipole generated in a spherical shell mimicking the Earth’s outer core?