zotero/storage/6DKHG2YF/.zotero-ft-cache

ARTICLE
https://doi.org/10.1038/s41467-020-16888-0

OPEN

Rapid geomagnetic changes inferred from Earth observations and numerical simulations

Christopher J. Davies 1✉ & Catherine G. Constable 2

Extreme variations in the direction of Earth’s magnetic ﬁeld contain important information regarding the operation of the geodynamo. Paleomagnetic studies have reported rapid directional changes reaching 1° yr−1, although the observations are controversial and their relation to physical processes in Earth’s core unknown. Here we show excellent agreement between amplitudes and latitude ranges of extreme directional changes in a suite of geodynamo simulations and a recent observational ﬁeld model spanning the past 100 kyrs. Remarkably, maximum rates of directional change reach ~10° yr−1, typically during times of decreasing ﬁeld strength, almost 100 times faster than current changes. Detailed analysis of the simulations and a simple analogue model indicate that extreme directional changes are associated with movement of reversed ﬂux across the core surface. Our results demonstrate that such rapid variations are compatible with the physics of the dynamo process and suggest that future searches for rapid directional changes should focus on low latitudes.

1234567890():,;

1 School of Earth and Environment, University of Leeds, Leeds LS2 9JT, UK. 2 Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography, University of California at San Diego, La Jolla, CA 92093-0225, USA. ✉email: c.davies@leeds.ac.uk

NATURE COMMUNICATIONS | (2020)11:3371 | https://doi.org/10.1038/s41467-020-16888-0 | www.nature.com/naturecommunications

1

ARTICLE

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16888-0

The large-scale secular variation of Earth’s internally generated magnetic ﬁeld is now reasonably well-established by
global models based on observations spanning the past two decades1 and the historical period2. Prominent features of these

models, such as the difference in activity between Atlantic and

Paciﬁc hemispheres and rapid changes at high latitudes, have

been linked, respectively, to thermal interactions between the core and mantle3,4 and accelerating jets in the core5, and hence to the

operation of the geodynamo. Over longer timescales, paleomagnetic studies6–8 have identiﬁed changes in ﬁeld intensity B in the Levantine region around 1000 BCE of ∂B/∂t ≈ 0.75–1.5 μT yr−1

that are signiﬁcantly faster than the largest values of ∂B/∂t ≈ 0.12 μT yr−1 for the modern ﬁeld9 and averages over the Holocene

ﬁeld10. Although the nature and origin of this ‘intensity spike’ has proved controversial11–13, recent estimates of its rate of change8,12 are compatible with bounds based on the kinematics of

core ﬂow14, and apparent inconsistencies in the available data

around 1000 BCE can be explained by large age uncertainties for

some samples15. Additional support for intensity spikes has come from studies in China16 and Texas17 and from numerical simu-

lations13 that reproduce similar values of ∂B/∂t. In the simula-

tions, intensity spikes reﬂect the migration of normal-polarity ﬂux patches across the core surface13.

Rapid changes in ﬁeld direction have also attracted signiﬁcant attention, particularly in the context of polarity reversals18. His-

torically, the fastest directional changes were attributed to lava

ﬂows at Steens Mountain, though these results are now thought to be untenable19. Currently, the fastest directional changes noted

are those recorded by sediments in central Italy20, where angular

fvcahaslatuneergsetoshfainn∂P^vtVha=elu∂etVs$irot1fua∂lyP^rGÀV1=e.o∂Tmt h$aegsn0e:e1triactyerPsÀo1alreefoparobtsohiuteitomna fo(adP^cetVor)nr

reach of 10 ﬁeld,

similar to the difference in rates of intensity change between the

modern ﬁeld and spikes. The reliability of the paleomagnetic

analysis of the Italian sediments is presently the subject of active

debate21,22. Overall, it has proved difﬁcult to substantiate extre-

mely rapid changes in individual paleomagnetic records either

during reversals or at other times. Furthermore, the compatibility

of these directional changes with the physics of the dynamo

process, and the dynamical origin of such variations, has not been

investigated.

Rapid directional changes have also been obtained in simple

kinematic ﬁeld models. Brown et al.23 independently reduced the

axial dipole ﬁeld coefﬁcient in the CALS7k.2 Holocene ﬁeld

model24 to simulate synthetic reversals. They found that reversal

duration varied signiﬁcantly with location, with certain localities

such as Tahiti experiencing rapid changes in the polarity of the

VGP. Fournier et al.11 constructed a synthetic model of an

intensity spike at mid-latitudes and noted that movement of the spike can change the dipole latitude by ~1° yr−1. However, it is

still not known whether these rates represent a “speed limit” for

the ﬁeld, whether such changes exhibit systematic geographical

preference, or how they relate to core dynamics.

In this study, we seek to establish how fast the local direction of

the geomagnetic ﬁeld might change in general and whether rapid

changes occur at preferred locations on Earth’s surface. Fur-

thermore, we investigate the processes at the core surface that give

rise to the most rapid directional changes. To do this, we take

advantage of recent developments in paleomagnetic ﬁeld mod-

elling and numerical dynamo simulations to document rates of

change during stable polarity times as well as excursions and

reversals. First, we use GGF100k25, the ﬁrst time-varying paleo-

magnetic ﬁeld model for the time interval 0–100 ka. Accessing

such long timescales is important because the fastest ﬁeld changes

are by deﬁnition rare events. Our results are bolstered by higher

resolution short term models spanning the Laschamp excursion. Second, we extend a set of tools that were developed to study rapid intensity changes in dynamo simulations13. This allows us to compare the whole spectrum of directional changes between simulations and observational ﬁeld models.
Results
Models and simulations of rapid directional changes. We search for extreme changes in ﬁeld direction in GGF100k and in a suite of 16 dynamo simulations that produce a range of dynamical behaviour and reproduce various features of the geomagnetic ﬁeld (see ref. 13 and Methods). Following previous work26,27, we compare simulations in terms of the magnetic Reynolds number Rm, which is used as the unique identiﬁer in both tables and ﬁgures. Our simulations span the range Rm = 100–700 (compared to the value for Earth’s core of Rm ~103), and some exhibit polarity reversals. To ensure robust results for comparing with the paleomagnetic analyses, we consider two ﬁeld properties, the local ﬁeld vector B and its transformation to the equivalent Virtual Dipole vector PV with unit vectors denoted by B^ and P^V , respectively. For both quantities, the corresponding rate of change (labelled ∂P^V =∂t and ∂B^=∂t) is calculated using the difference between values at time t and t + Δt. We ﬁnd the maximum rate of change of directions at each location on a 2° by 2° latitude–longitude (λ, ϕ) geographic grid (Methods). The largest changes anywhere on the globe are referred to as ‘extremal events’ and are denoted with a subscript ex.
Figure 1 shows the geographic variability in ∂P^V =∂t and timeseries at the locations of ð∂P^V =∂tÞex for one simulation (Rm = 386) that undergoes excursions but not reversals, one reversing simulation (Rm = 450), and the 0–100 ka time-varying ﬁeld model GGF100k25. (Time-series for all 16 simulations studied are provided in Supplementary Fig. 1). The simulations shown in Fig. 1 have been run for many dipole decay times (see Supplementary Table 1 for details) and have complete spatial and temporal coverage, thus providing a statistically representative set of extreme rates of change that occur at different times and places across the globe. The dominant event for GGF100k is around the Laschamp excursion at 41 ka. The regional focus and relatively smooth time variation in GGF100k reﬂect both the shorter time span, and the uneven temporal and spatial coverage of the underlying paleomagnetic data. Note that the temporal sampling, Δt, and size of the spatial grid are chosen to ensure that our estimates of extreme rates of change are conservative. Experimentation with denser spatial and temporal sampling and with LSMOD.2, a shorter higher resolution model28 spanning the Laschamp excursion (Supplementary Fig. 8), suggests the possibility of extreme values a factor of 2–5 larger than those shown for GGF100k. Nevertheless, both reversing and nonreversing simulations as well as GGF100k exhibit rapid directional changes with ∂P^V =∂t > 1 yrÀ1 (see also Fig. 2).
Surprisingly, most of our simulations produce extremal directional changes that are removed by more than 5 kyr from reversals (see, e.g. Fig. 1, middle panel, where for Rm = 450, the extremal change precedes the equatorial crossing of the axial dipole by at least 5 kyr). Similar behaviour occurs in GGF100k, where the extremal directional change occurs slightly after the main phase of the Laschamp excursion. Another common feature of extreme directional changes is a considerably diminished overall ﬁeld strength, but not necessarily a local minimum in B(t) (Fig. 1). It is notable that while all reversing simulations produce high rates of directional change >1° yr−1, even non-reversing simulations can still produce similarly high values associated with global or regional excursions and intensity minima.

2

NATURE COMMUNICATIONS | (2020)11:3371 | https://doi.org/10.1038/s41467-020-16888-0 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16888-0

ARTICLE

* ˆB /*t (° yr–1)

Non-reversing simulation Reversing simulation GGF100k

a
4.94 3.8 2.92 2.25 1.73 1.33 1.02 8e-01 6e-01 5e-01
b
8.53 6.89 5.56 4.49 3.63 2.93 2.37 1.91 1.54 1.25
c
2.7 1.67 1.03 6e-01 4e-01 2e-01 1e-01 9e-02 6e-02 3e-02

* ˆPV /*t (° yr–1)

* ˆPV /*t (° yr–1)

* ˆPV /*t (° yr–1)

* ˆPV /*t (° yr–1)

V, D (°)

* ˆPV /*t (° yr–1)

V, D (°)

* ˆPV /*t (° yr–1)

V, D (°)

d

B ( T)

90

40

60

35 30

30

25

0

20

–30 –60

15 10 5

–90

0

3.5

3.5

3

Rm = 386 3

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0

425 427 429 431 433 435 437 439 441 443 445

Time t (kyrs)
e

B ( T)

90

90

60

80 70

30

60

0

50 40

–30

30

–60

20 10

–90

0

3.5

3.5

3

Rm = 450 3

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0

543 545 547 549 551 553 555 557 559 561 563

Time t (kyrs)

90

35 f

60

30

B ( T)

30

25

0

20 15

–30

10

–60

5

–90

0

3.5

3.5

3

3

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0

–49 –47 –45 –43 –41 –39 –37 –35 –33 –31 –29

Time t (kyrs)

* ˆB /*t (° yr–1)

* ˆB /*t (° yr–1)

Fig. 1 Rapid directional changes in two geodynamo simulations and GGF100k. Shown are a non-reversing simulation with magnetic Reynolds number
Rm = 386 (a, d), a reversing simulation with Rm = 450 (b, e) and the observational ﬁeld model GGF100k25 (c, f). Left column shows Mollweide projections at Earth's surface of the largest change in VGP position, ∂P^V=∂t, as a function of location in ° yr−1. Red stars show the location of ð∂P^V=∂tÞex on each plot. Note the different colour scales and that values at each location may not have occurred at the same point in time. Right-hand panels show directional data at the locations of ð∂P^V=∂tÞex over a 20 kyr period with the extreme event at the midpoint. Here the top row of each panel shows the latitude λV of P^V (purple), the dipole latitude λD (black), and the ﬁeld strength B (red); the bottom row shows ∂P^V=∂t (purple) and the rate of change of the ﬁeld vector ∂B^=∂t (blue). Simulations have been run for 232 kyrs for Rm = 386 and 415 yrs for Rm = 450.

The time-series of ð∂P^V =∂tÞex and ð∂B^=∂tÞex (Fig. 1, Supplementary Fig. 1) show a clear spike at the times of extreme directional changes. However, the shapes of the corresponding P^V and B^ time-series vary across simulations: some are spike-like, others have steep increases and decreases with a relatively ﬂat
intermediate period, and occasionally the time-series does not
show both a rise and fall. We think these distinctive spike-like
patterns in the time-derivatives are best described as local or regional rapid accelerations in magnetic ﬁeld changes. This view is supported by the fact that large values of ð∂P^V =∂tÞex and ð∂B^=∂tÞex are not necessarily associated with low values of λD, the latitude of the geomagnetic dipole axis, as commonly seen in reversals or ﬁeld excursions, although in both GGF100k and the
various simulations they do generally seem to occur at times when ﬁeld strength is low enough that non-dipole ﬁeld activity
may dominate. We will just call them extremal directional
changes. Figure 2a shows the maximum rates of change in ﬁeld direction
for all simulations. Remarkably, extremal directional changes can

reach 10° yr−1 in the simulations. This is ten times faster than the
fastest variations so far reported in individual paleomagnetic records20 and about 40 times faster than the fastest changes in Holocene ﬁeld models CALS10k.2 and HFM.OL1.A110. However,
the examples in Fig. 1 are quite compatible with our conservative estimates of peak rates of change of 2.5–3.5° yr−1 found in
GGF100k and LSMOD.2 in the millennia surrounding the
Laschamp excursion when the dipole moment is very low.
The absolute latitude corresponding to the most extreme directional changes, ∣λex∣, is generally <40° in all simulations (Fig. 2b). Using either B^ or P^V to measure directional changes gives broadly similar results, though stronger latitudinal variability arises when considering B^ compared to P^V , with noticeably faster directional changes tending to occur at lower latitudes. Figure 3a shows the histogram of ∂B^=∂t values for GGF100k,
which is well approximated by a log-normal distribution both globally and at individual latitudes. The corresponding ﬁts to
cumulative distribution functions (CDFs) at different latitudes for
GGF100k (Fig. 3b) clearly show that the probability of observing

NATURE COMMUNICATIONS | (2020)11:3371 | https://doi.org/10.1038/s41467-020-16888-0 | www.nature.com/naturecommunications

3

ARTICLE

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16888-0

a

10

10

(* ˆB /*t )ex (° yr–1)

(* ˆPV /*t )ex (° yr–1)

1

1

0.1

0.1

200 400 600 800 1000

Magnetic Reynolds number Rm

ex( ˆB) (°)

ex( ˆPV) (°)

b 80

80

60

60

40

40

20

20

0

0

–20

–20

–40

–40

–60

–60

–80

–80

200 400 600 800 1000

Magnetic Reynolds number Rm

Fig. 2 Amplitudes and locations of extreme directional changes in all simulations and GGF100k. The ﬁeld vector is denoted B^ and the VGP position is denoted P^V. In b the locations of maximum rates of change of B^ and P^V are λexðB^Þ and λexðP^VÞ, respectively. Closed (open) symbols indicate non-reversing (reversing) simulations. Where shown, values for Rm = 1000 are averages from the Holocene ﬁeld model CALS10k.210 (solid symbols) and GGF100k25 (open symbols), otherwise the values
plot below the lower limit on the ordinate. All locations are given in
degrees.

faster rates of directional change is greatest at the lower latitudes.
A synthesis of the scale and shape information presented in Fig. 3b is shown in Fig. 3c, d, which show ratios of the mean (μ) and standard deviation (σ) of the log-normal distributions at 40° and
80° latitude to the values at 0° for GGF100k and all of the dynamo simulations. Normalised mean values decrease signiﬁcantly with
increasing latitude, indicating that the CDF curves at lower latitudes are shifted towards higher values of ∂B^=∂t, thus
following the behaviour of the GGF100k model (Fig. 3b). Standard deviations (reﬂecting relative variability in ∂B^=∂t)
decrease with latitude in most simulations, as in GGF100k, whereas a few simulations with Rm ≈ 200–300 show maximum values of σ at mid-to-low latitudes. Overall, the excellent
agreement between GGF100k and the dynamo simulations shows
a general preference for rapid directional changes at lower latitudes. These variations in B^ are readily linked to core
processes as will be seen in the next section, where we outline
the strategy for investigating the origin of rapid directional
changes.

Linking rapid directional changes and core processes. We ﬁrst

consider a simple model of rapid directional changes in which the

radial CMB ﬁeld consists of an isolated moving ﬂux patch

superimposed on a static axial dipole ﬁeld deﬁned by the Gauss

coefﬁcient

g

0 1

¼

80

μT.

The

spatial

form

of

the

patch

is

described

by a Fisher–von Mises probability distribution function with

amplitude A and half width σM = 15° (Methods). The patch is moved in latitude (longitude) using constant intervals of δλ (δϕ).

pdf

a
101
10–1
10–3
10–5 0

pdf

102 101 100 10–1 10–2 10–3 10–4 10–5
0.0

c
= 2.5° = 22.5° = 42.5°

0.5

1.0

1.5

*Bˆ /*t (° year–1)

1

2

*Bˆ /*t (° yr–1)

Ratio of means at latitudes X and Y , X / Y

0.9

40°/ 0°

0.8

80°/ 0°

0.7

0.6

0.5

0.4

0.3

0.2 200 400 600 800 1000 Rm

cdf

b
1.0 0.8 0.6 0.4 0.2 0.0
10–3

= –37.5° = –17.5° = 2.5° = 22.5° = 42.5° = 62.5° = 82.5°

10–2

10–1

100

*Bˆ /*t (° yr–1)

d 1.10
1.05 1.00

40°/ 0° 80°/ 0°

Ratio of at latitudes X and Y , X / Y

0.95

0.90

0.85

0.80

0.75 200 400 600 800 1000 Rm

Fig. 3 Statistics of extreme directional magnetic ﬁeld changes. a Normalised histogram for rates of change of the ﬁeld vector B^ in the GGF100k model25 sampled on a 2° × 2° latitude–longitude grid at 100 year time intervals. The black line shows the ﬁt of a log-normal probability density function (PDF) to the data with mean μ = 0.009 and standard deviation σ = 0.011. The inset shows histograms at latitudes of 2.5° (red), 22.5° (green) and 42.5° (blue) with PDFs shown as solid lines. b Cumulative distribution functions (CDFs) of ∂B^=∂t at different latitudes (shown by colours) from the GGF100k model, showing higher rates of change at lower latitudes. c The ratio of the mean values of PDFs at two different latitudes (denoted by subscripts) for GGF100k and all
dynamo simulations as a function of the magnetic Reynolds number Rm. GGF100k is shown at Rm = 1000 as in Fig. 2. Panel d is the same as c but showing ratios of standard deviations.

4

NATURE COMMUNICATIONS | (2020)11:3371 | https://doi.org/10.1038/s41467-020-16888-0 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16888-0

ARTICLE

a 200.00
159.18 118.37 77.55 36.73 –4.08 –44.90 –85.71 –126.53 –167.35
16.73 12.65 8.57 4.49 0.41 –3.67 –7.76 –11.84 –15.92 –20.00
83.97 80.88 77.80 74.71 71.63 68.54 65.45 62.37 59.28 56.20

Z ( T)

Z ( T)

VGP latitude

200.00 159.18 118.37 77.55 36.73 –4.08 –44.90 –85.71 –126.53 –167.35
16.73 12.65 8.57 4.49 0.41 –3.67 –7.76 –11.84 –15.92 –20.00
83.18 80.02 76.86 73.70 70.54 67.38 64.22 61.06 57.90 54.74

Z ( T)

Z ( T)
* ˆPV /*t (°)

b
20

15

10


 = 10°

5


 = 20°

A = 10 (R)

A = –10 (N) 0
0 10 20 30 40 50 60

ex (Pˆ V )

c
20

= –10° (Equatorward) = 5° (Poleward)

15

VGP latitude
* ˆPV /*t (°)

11.51 10.23 8.96 7.68 6.40 5.12 3.84 2.56 1.29 0.01

* ˆPV /*t

9.92

10

8.82

7.72

* ˆPV /*t

6.61

5.51

5

4.41

3.31

2.21

1.11

0

0.01

–20 –10 0 10 20 30 40

ex (Pˆ V )

Fig. 4 Directional changes produced by a simple model of a patch of radial ﬁeld moving across the core surface. a The effect of a reversed ﬂux patch of
amplitude A = 10 μT with an initial latitude of 60° moving southward to 50° (left) and then to 40° (right) (δλ = −10°). Rows show, respectively, the vertical magnetic ﬁeld at the CMB and surface, Z, VGP latitude at the surface, and the rate of change of VGP position ∂P^V=∂t. The centre of the patch on the CMB is denoted by the star. b, c The amplitudes, ∂P^V=∂t, and locations, λexðP^VÞ, of the fastest directional changes for normal (A < 0, denoted N) and reversed (A > 0, denoted R) polarity patches. In b, patches initially at 5°, 20°, 40° and 60° latitude are moved longitudinally in increments of δϕ = 10°
(squares) and 20° (circles). In c patches initially placed at 40° longitude are moved equatorward and poleward in increments of δλ = −10° and δλ = 5°, respectively. In all cases in b, c ∂P^V=∂t is largest for reversed patches.

Time increments are arbitrarily set to Δt = 1 yr and so this model
cannot provide an absolute rate of directional change; nevertheless it allows us to systematically compare the inﬂuence of
isolated normal (A < 0) and reversed (A > 0) CMB patches on
directional changes at the surface without the complication of
multiple patches as arises in the simulations.
Figure 4 shows an example evolution of the CMB and surface ﬁelds and a synthesis of models with normal and reversed patches
moving in various directions. The main conclusion is that
reversed patches produce faster directional changes than normal
patches, all other factors being equal (compare red, reverse, and
green, normal, points in Fig. 4b, c). Interestingly, normal-polarity patches moving towards the pole can induce a change in ∂P^V =∂t in the opposite hemisphere that is almost as large as the change
induced by the reversed patches in the same hemisphere (Fig. 4c, circles). However, in a more realistic scenario where the ﬂux
patch executes some combination of the motions shown in Fig. 4
the fastest changes, which are the main focus here, would clearly arise from reversed ﬂux migration. In this model the extremal events arise in regions where the amplitude of the reversed ﬂux balances that of the dipole ﬁeld at the surface (Fig. 4a),
corresponding to a region of low inclination. The non-dipole ﬁeld, which is weaker and more variable than the dipole ﬁeld,
then controls the rate of directional change. This simple model suggests that reversed ﬂux patches could be
relevant for understanding rapid directional changes in the simulations. Figure 5 shows snapshots of λV, the latitude of P^V , together with the vertical magnetic ﬁeld at the core-mantle

boundary (CMB) and at Earth’s surface for a representative
simulation (see Supplementary Information for other cases). The extremal event arises as a reversed ﬂux patch moves southwards
beneath the northeast Asian region. This patch, which is clearly visible on the CMB and also apparent in the surface ﬁeld, creates
a local inclination anomaly of opposite sign compared to the surrounding ﬁeld, which is reﬂected in λV.
To investigate this behaviour in more detail, we have developed a tool for locally reducing the ﬁeld strength in regions
surrounding extreme events (Methods). We consider four
30° × 30° quadrants with a common vertex at the location of the extremal event and halve the CMB ﬁeld strength in each
quadrant in turn. We then compute directional changes at the surface based on these CMB ﬁelds as before. Results from three simulations are shown in Fig. 6 (see Supplementary Figs. 2–6 for
more examples). Figure 6a shows the result of performing this
analysis on the simulation in Fig. 5, which reveals that the northeast quadrant containing the large reversed ﬂux patch
makes the biggest contribution to the observed rate of change. In Fig. 6b, reversed ﬂux patches in the northwest and southwest
quadrants that lie either side of the equator are the cause of the
rapid change. In Fig. 6c, the largest contributions to the directional change arise from a reversed ﬂux patch moving
towards the observation point from the northeast. This analysis is somewhat simpliﬁed because it does not
account for the changing direction and strength of different
patches. Possible extensions to our method include tracking individual patches29 and using time-dependent factors for

NATURE COMMUNICATIONS | (2020)11:3371 | https://doi.org/10.1038/s41467-020-16888-0 | www.nature.com/naturecommunications

5

ARTICLE

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16888-0

VGP latitude

Z ( T)

66.9 50.6 34.3 18.0 1.6 –14.7 –31.0 –47.3 –63.7 –80.0
0.02 0.01 0.01 0.00 0.00 –0.00 –0.01 –0.01 –0.02 –0.02
0.42 0.32 0.21 0.11 0.01 –0.09 –0.19 –0.30 –0.40 –0.50

Z ( T)

Z ( T)

VGP latitude

66.9 50.6 34.3 18.0 1.6 –14.7 –31.0 –47.3 –63.7 –80.0
0.02 0.01 0.01 0.00 0.00 –0.00 –0.01 –0.01 –0.02 –0.02
0.42 0.32 0.21 0.11 0.01 –0.09 –0.19 –0.30 –0.40 –0.50

Z ( T)

Z ( T)

VGP latitude

66.9 50.6 34.3 18.0 1.6 –14.7 –31.0 –47.3 –63.7 –80.0
0.02 0.01 0.01 0.00 0.00 –0.00 –0.01 –0.01 –0.02 –0.02
0.42 0.32 0.21 0.11 0.01 –0.09 –0.19 –0.30 –0.40 –0.50

Z ( T)

Z ( T)

0.13

0.13

0.09

0.09

0.06 0.03

Reversed

0.06 0.03

Z ( T)

0.00 flux patch

0.00

–0.03

–0.03

–0.06

–0.06

–0.09

–0.09

–0.12

–0.12

–0.15

–0.15

Z ( T)

0.13 0.09 0.06 0.03 0.00 –0.03 –0.06 –0.09 –0.12 –0.15

Fig. 5 Evolution of an extreme change in the VGP position, P^V , for the simulation with E = 5 × 10−4, Pm = 5, Ra = 400 and Rm = 257. Rows show the VGP latitude λV, vertical component of the magnetic ﬁeld Z at Earth's surface, and at the CMB, and a local Mercator projection of the CMB ﬁeld around the location of maximum change (white star). Columns show times just before (left), during (centre) and just after (right) the extreme change. The extremal
event arises as a reversed ﬂux patch (blue) migrates from the northeast towards the point of maximum change.

reducing the local ﬁeld strength as the simulation evolves, though both options have their drawbacks. Nevertheless, we believe that the results in Fig. 6, combined with the results from the simple patch model (Fig. 4), present a strong link between rapid directional changes at Earth’s surface and the movement of reversed ﬂux patches across the core surface.
Discussion Our suite of geodynamo simulations span a wide range of physical parameters and access reversing and non-reversing dynamo regimes; however, as with all current models, they still do not reach the low values of the Ekman number E and magnetic Prandtl number Pm that characterise Earth’s core. However, all exhibit some Earth-like features, and we ﬁnd consistent results for the locations and amplitudes of extremal events and no obvious dependence on Rm. As discussed in detail in our previous work on intensity spikes13, there is no reason to believe that the rates of change in these simulations should be greater or smaller than those that would be obtained at more extreme parameters. Indeed, we ﬁnd that different values of E and Pm produce similar rates of change in the same dynamo regime (e.g. stable dipolar or reversing).
The GGF100k model is limited by the spatial and temporal distribution of the data and the accuracy of the chronological constraints at different locations. This produces a smoothed record in time that can only resolve large-scale spatial features. We expect these limitations to result in estimates of ﬁeld changes that are lower than the actual variations. This is conﬁrmed by the fact that higher rates of change are evident in LSMOD.1 and LSMOD.2, two recent higher resolution global time-dependent

ﬁeld models spanning 50–30 ka28,30 that incorporate updated chronological information. The use of the longer record provided by GGF100k reveals that the most extreme values occur at times of low ﬁeld strength and the approximately log-normal distributions of rates of change further demonstrate that extreme values are relatively rare over the past 100 kyr. From these considerations and the excellent consistency between the locations and amplitudes of extreme directional changes in both observational and numerical models, which are entirely independent, we believe our results provide strong evidence that the geodynamo can produce much faster directional changes than have previously been considered viable.
The fastest directional changes are of comparable magnitude in both simulations and paleoﬁeld models, reaching ~10° yr−1 in simulations and rising to 4.8° and 22.5° yr−1 in GGF100k and LSMOD.2, respectively, when these models are more densely sampled at an interval of 10 years. This is faster than the latest published local paleomagnetic estimates20 of 1° yr−1, and suggests that such rapid directional changes must be fully compatible with the physics of the dynamo process. Of course, such events are exceedingly rare: in the 100 kyr spanned by GGF100k, we found only one time interval with directional changes of >2° yr−1 from over 5 million samples of the ﬁeld, whereas 95% of the samples have rates of change below 0.05° yr−1 (Fig. 3). However, faster rates of directional change tend to occur at lower latitudes where the ﬁeld is weaker. Our results suggest that the most rapid directional changes occur at latitudes <40°, which could be helpful in guiding future paleomagnetic acquisitions.
The best evidence for rapid directional changes in time-varying paleoﬁeld models comes from times that are close to known excursions. Although rapid changes are not necessarily associated

6

NATURE COMMUNICATIONS | (2020)11:3371 | https://doi.org/10.1038/s41467-020-16888-0 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16888-0

ARTICLE

Z (µT)

60°E

90°E

120°E

60°N

30°N

0.13 0.09 0.06 0.03 0.00 –0.03 –0.06 –0.09 –0.12 –0.15

Z (µT)

0° 30°S

90°W

60°W

0.17 0.13 0.09 0.04 0.00 –0.04 –0.08 –0.12 –0.16 –0.20

Z (µT)

0° 30°S 60°S

90°E

120°E

150°E

0.13 0.09 0.06 0.03 0.00 –0.03 –0.06 –0.09 –0.12 –0.15

Z (µT)

60°E

90°E

120°E

0.13

0.09

60°N

0.06

0.03

0.00

–0.03

–0.06

30°N

–0.09 –0.12

60°E

90°E

120°E

–0.15

Reversed flux

NW NE SW SE Original

4
a 3.5
3 2.5 2

* ˆPV /*t (° yr–1)

1.5

1

0.5

0 688.7 688.8 688.9 689 689.1 689.2
t (kyrs)

Z (µT)

90°W

60°W

0°

30°S

90°W

60°W

Z (µT)

90°E
0° 0.17
0.13
0.09
0.04 30°S
0.00
–0.04
–0.08
–0.12
–0.16 –0.20 60°S
90°E
Reversed flux

120°E 120°E

150°E 150°E

0.13 0.09 0.06 0.03 0.00 –0.03 –0.06 –0.09 –0.12 –0.15

* ˆPV /*t (° yr–1)

6
b
5
4
3
2
1
0 435 435.2 435.4 435.6 435.8 436
t (kyrs)

* ˆPV /*t (° yr–1)

552.8

9

NW NE

c8

SW

7

SE

6

Original

5

4

3

2

1

0 552.85 552.9 552.95 553

t (kyrs)

Fig. 6 Contributions of CMB ﬁeld features to rates of directional change at Earth’s surface. The three columns show extreme changes arising in the northern hemisphere (left), equatorial region (centre) and southern hemisphere (right) in simulations with E = 5 × 10−4, Ra = 400, Pm = 5, Rm = 257 (a); E = 5 × 10−4, Ra = 250, Pm = 10, Rm = 386 (b); E = 5 × 10−4, Ra = 350, Pm = 10, Rm = 450 (c). Rows show local Mercator projections of the radial CMB ﬁeld just before (top) and during (middle) an extreme directional change located at the star. The bottom row shows the rate of VGP change, ∂P^V=∂t, at the location of the star based on the original CMB ﬁeld and four modiﬁed versions of the CMB ﬁeld wherein the ﬁeld in the denoted quadrant has been
halved in amplitude. In this projection the four quadrants denoted northwest (NW), northeast (NE), southwest (SW) and southeast (SE) tessellate the
regions shown in the ﬁrst two rows and share a common vertex centred at the star.

with large deviations of the dipole axis, they are often associated with a decrease in ﬁeld strength that accompanies excursions25,30. However, we do not yet rule out the possibility that large directional changes might occur during stable polarity times due to rapid growth of a reverse ﬂux patch that might not lead to global large-scale decrease in intensity. Such a scenario is seen in a number of our dynamo simulations (e.g. Supplementary Fig. 1e, g, h, k) and may correspond to a so-called regional rather than global excursion. The quality of the current paleomagnetic record is not adequate to map such regional behaviour in detail, but recent progress in global ﬁeld modelling suggests it can be achieved in the future.
In the simulations, extremal directional changes arise from the migration of reversed ﬂux patches towards the observation site, which produces a sharp change in direction. Reversed ﬂux patches tend to occur more frequently in the equatorial region, which explains the observation that rapid directional changes occur at low latitudes. In our simple model of an isolated CMB

ﬂux patch the fastest directional changes are always associated with migration of reversed ﬂux on the CMB; however, normal ﬂux patches seem likely to have a general role in producing slower
variations. Examining the relative contributions of normal and reversed CMB ﬂux across the statistical distribution of directional changes is a signiﬁcant undertaking owing to the detailed analysis
required for each event, but could be approached based on the
tools developed here to study extremal events.
We can compare the characteristics of directional changes to those of intensity variations, which were previously analysed13 in
9 of the 16 simulations presented here. Rapid and localised
changes in intensity comparable to the rates inferred for the
Levantine spike were found to occur predominantly at latitudes >50° and were produced by rapid migration and intensiﬁcation of normal-polarity ﬂux patches on the core surface. These results suggest that extremal intensity and directional variations reﬂect
different physical processes at the top of the core. It is of course
well-known that the local surface inclination and declination

NATURE COMMUNICATIONS | (2020)11:3371 | https://doi.org/10.1038/s41467-020-16888-0 | www.nature.com/naturecommunications

7

ARTICLE

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16888-0

sample a different part of the CMB from the intensity. For an axial dipole ﬁeld inclination sampling at the CMB is biased toward lower latitudes (except at the equator), declination preferentially samples east and west of the sampling location, and local surface intensity sampling peaks poleward of the sampling site31. Thus spatial sampling bias may further reinforce the ability to detect latitudinal differences between extremal events in direction and intensity. Further work is needed to assess whether under some circumstances accounting for this sampling bias may allow any information about intensity spikes to be extracted from local surface directional changes.
New global time-dependent ﬁeld models for periods preceding the Holocene25,28,30 are revealing geomagnetic ﬁeld behaviour that is not represented in the relatively short historical record. The strong agreement between rates and locations of rapid directional change in these models and numerical simulations means that it is now possible to probe the dynamical origin of these extreme events. It is plausible that the characteristics of extreme events contain information about the structure and dynamics of the CMB region. A stably stratiﬁed layer in the uppermost core, which has been inferred from seismology32 and may reﬂect thermal33 or chemical34,35 interactions with the mantle, may suppress radial motions and weaken the ﬁeld strength at the core surface36. On the one hand a weaker ﬁeld is generally associated with increased rates of change in our simulations; however, we also ﬁnd that extremal directional changes arise from horizontal motion of ﬂux patches, which is not hindered by stratiﬁcation. Lateral variations in CMB heat ﬂow can induce longitudinal variations in the structure of the magnetic ﬁeld37 and secular variation38 and may even produce regional stratiﬁcation39. Although we observe no clear differences between the 5 simulations that include CMB heat ﬂow variations and the 11 that do not, this sample is not large enough to allow robust conclusions on this point. Future work linking numerical simulations to global models of the paleomagnetic ﬁeld can help resolve these issues.
Methods
Geodynamo simulations. Full details of the numerical simulations are given in our previous papers13,27 and so only a brief description is provided here. Our computer code solves the equations that describe conservation of mass, momentum, energy and magnetic ﬂux in a rotating spherical shell. In the dimensionless governing equations the input parameters are the Ekman number E, measuring the ratio of viscous to Coriolis effects; the magnetic Prandtl number Pm, the ratio of viscous and magnetic diffusivities; the Rayleigh number Ra, measuring the vigour of convection; the Prandtl number Pr = 1, the ratio of viscous and thermal diffusivities; the ratio of inner to outer boundary radii, ri/ro = 0.35; and the amplitude of boundary heat ﬂow heterogeneity q⋆ (=0 for homogeneous boundaries). At ro the boundary conditions are no-slip, electrically insulating and ﬁxed heat ﬂux (FF), whereas at ri the boundary conditions are no-slip, electrically conducting or insulating, and either ﬁxed temperature (FT) or ﬁxed ﬂux. Simulations with a laterally varying heat ﬂux at ro use the pattern from ref. 40. The magnetic Reynolds number Rm, the ratio of advection and diffusion of magnetic ﬁeld, serves to uniquely identify the simulations.
We consider a total of 16 simulations, which are listed in Supplementary Table 1. The simulations with E = 5 × 10−4 and 1.2 × 10−4 were originally published in ref. 27, the simulations with E = 10−5 were published in ref. 41 and the simulations with E = 10−3 were published in ref. 42. These simulations were chosen because they have been shown to reproduce various features of Earth’s magnetic ﬁeld; however, they have never been used to study local directional changes.
Analysis of extreme directional changes. The geodynamo simulations use dimensionless variables (denoted by stars below), which must be converted to physical units in order to be compared with geomagnetic observations. Following our previous work13, time t is scaled using the advection time, t = (Rmm/RmE) t⋆d2/η, where d = 2264 km is the shell thickness, η = 1 m2 s−1 is the magnetic diffusivity and Rmm and RmE are the magnetic Reynolds number of the model and Earth, respectively. Intensity is scaled based on the polar value of the time-averaged ﬁeld strength B, which gives similar results to other choices13.

At each time the three components X, Y and Z of the local magnetic ﬁeld vector B are calculated on a 2° by 2° latitude–longitude grid at Earth’s surface as described
ipn ﬃrﬃﬃeﬃﬃfﬃ.ﬃﬃﬃ1ﬃﬃ3ﬃﬃ.ﬃﬃFﬃﬃﬃrﬃﬃoﬃﬃmﬃﬃﬃﬃﬃtﬃﬃhﬃﬃese components, the strength of the local ﬁeld vector B ¼ ðX2 þ Y2 þ Z2Þ and the inclination I and declination D are obtained. Deﬁning the instantaneous units vectors

X^ðtÞ ¼ XðtÞ=BðtÞ; Y^ðtÞ ¼ YðtÞ=BðtÞ; Z^ðtÞ ¼ ZðtÞ=BðtÞ;

ð1Þ

the rate of change of B^ between times t and t1 = t + Δt is given by

∂B^ ∂t

¼

arccos½X^ ðtÞX^ðt1Þ

þ

Y^ ðtÞY^ðt1Þ Δt

þ

Z^ ðt ÞZ^ ðt 1 Þ

ð2Þ

and the rate of change of B is

∂B ∂t

¼

Bðt1Þ À Δt

BðtÞ

:

ð3Þ

Variations of the virtual dipole vector PV are calculated in the same manner as variations in B. The virtual dipole position P^V with latitude λV and longitude, ϕV is calculated from I and D and converted to unit vectors according to

X^ V ¼ cosðλV Þ cosðϕV Þ; Y^ V ¼ cosðλV Þ sinðϕV Þ; Z^V ¼ sinðλV Þ:

ð4Þ

The rate of change of P^V between times t and t1 is then given by

∂P^V ∂t

¼

arccos½X^V ðtÞX^ V ðt1Þ

þ Y^V ðtÞY^V ðt1Þ þ Δt

Z^V ðtÞZ^V ðt1Þ :

ð5Þ

The amplitude of the virtual dipole vector PV is the virtual dipole moment, which is calculated according to43

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

PV ¼ 4πa3

ð1 þ 3cos2ðI ð2μ0 Þ

ÞÞB

;

ð6Þ

where a = 6371 km is the radius of the Earth, μ0 is the permeability of free space and B is in micro Tesla. The rate of change of PV is then obtained from

∂PV ∂t

¼

PV

ðt1

ÞÀ Δt

PV

ðtÞ

:

ð7Þ

A simple model of rapid directional changes. We generalise the model in ref. 15 in which an isolated patch is inserted into the radial magnetic ﬁeld Br at the CMB, radius r = c = 3485 km, centred at the point θc, ϕc in spherical polar coordinates (r, θ, ϕ). The functional form is based on the Fisher–Von Mises probability density
function and is given by

Bsr ðc;

θ;

ϕÞ

¼

Aκ

expκ cos θ 4π sinh κ

À

A 4π

;

ð8Þ

where

cos θ ¼ sin θ cos ϕ sin θc cos ϕc þ sin θ sin ϕ sin θc sin ϕc þ cos θ cos θc ð9Þ and

κ

¼

6561

σ

2 M

;

ð10Þ

where σM is the standard deviation in both latitude and longitude. The factor A/4π ensures that the patch ﬁeld deﬁned by Eq. (8) remains divergence-free, i.e.
Z

BsrdSc ¼ 0;

ð11Þ

where Sc denotes the core surface. The patch ﬁeld Bsr can be expanded in a convergent spherical harmonic series
on the unit sphere. It is then straightforward to add to Bsr a background ﬁeld Bbr such that the total radial ﬁeld Br ¼ Bsr þ Bbr . Here we set Bbr equal to an axial dipole ﬁeld. The total ﬁeld on the CMB can then be upward continued to the surface
where magnetic elements X, Y, Z, intensity F, inclination I, declination D can be
calculated at all points on a global grid. The aforementioned model produces a stationary ﬁeld on the CMB. Here we
generalise this model to allow the patch to move across the core surface. This is accomplished by introducing two additional parameters that deﬁne the latitudinal and longitudinal increments to the central patch location, δλc and δϕc, respectively. Over a sequence of n = 1, …, N iterations the centre of the patch on the CMB is moved from its original location (θc, ϕc) to the locations (θc + nδλc, ϕc + nδϕc). Of course, there is no physical signiﬁcance to the increments used in the calculations
and hence absolute rates of change cannot be calculated by this approach. Therefore, we simply set Δt = 1 in Eqs. (2) and (5). We calculate ∂B^=∂t and ∂P^V =∂t at all locations on a 2° × 2° latitude–longitude grid at Earth’s surface and compare values for imposed normal and reversed CMB patches.
The generalised model of an isolated CMB ﬂux patch is completely speciﬁed by the patch amplitude A, its central location and increments, θc, δλc, ϕc, δϕc, and width σM. We ﬁx σM = 15°, which produces a broad patch on both the CMB and surface that is easy to identify in visualisations. We set a background axial dipole ﬁeld with an amplitude of 80 μT at the CMB, similar to the present geomagnetic

8

NATURE COMMUNICATIONS | (2020)11:3371 | https://doi.org/10.1038/s41467-020-16888-0 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16888-0

ARTICLE

ﬁeld, and choose A = −10 μT and A = 10 μT to simulate the two cases of normal
and reversed ﬂux patches, respectively. Values of θc, δλc, ϕc and δϕc are varied for the two cases. For the chosen background ﬁeld, the solutions are identical for
positive and negative values of δϕc and so we focus on the former. We consider both positive and negative values of δλc to simulate patches moving towards and
away from the pole, respectively.

Masking features of the CMB ﬁeld. The goal here is to modify the local strength of the radial CMB ﬁeld in order to evaluate the effect of speciﬁc ﬁeld features on

directional changes at the surface. We refer to this process as “masking”. Local

manipulation of the ﬁeld must be done in real space; however, the ﬁeld must also

be represented in spherical harmonics in order to perform upward/downward

continuation between the surface and CMB. Therefore, any local change to the ﬁeld

in real space will alter the structure globally when represented in spherical har-

monics, leading to the possibility of aliasing. To mitigate this effect, we retain all

spherical harmonic coefﬁcients in the expansion of the radial ﬁeld up to the truncation used in the original simulations, usually degree L = 128 − 256. For each

simulation, sets of high-resolution Gauss coefﬁcients are produced for the 5 time

points surrounding an extreme change in the VGP.

The Schmidt-normalised Gauss coefﬁcients gml of degree l and order m at Earth’s surface are related to the coefﬁcients cml for the radial ﬁeld at the CMB by

cml

¼

ðl

þ

1Þg

m l

alþ1 c

:

ð12Þ

The cml are transformed to and from real space using Shtools44 with Gauss–Legendre quadrature and (3/2)m latitude points to ensure an exact transform. In real space, a spatial region is identiﬁed and the ﬁeld there is

multiplied by a factor f. After experimentation with values of 0.01 ≤ f ≤ 0.9 the value

f = 0.5 was chosen as this is approximately the fractional difference between the

maximum and mean absolute ﬁeld strength in many simulations and also ensures

that the ﬁeld structure remains smooth with a convergent spherical harmonic

representation.

In each simulation, we mask four different regions in turn, chosen to be the four

30° × 30° quadrants with a common vertex at the location of maximum change.

The quadrants will be referred to as NW (northwest), NE (northeast), SW

(southwest) and SE (southeast). For each of the 5 sets of Gauss coefﬁcients (the

original plus one each for the 4 quadrants), we compute directional changes as

described above, creating maps of X, Y and Z on Earth’s surface and using these to

extract rates of change at all locations.

An example calculation (Supplementary Fig. 7) shows that the method is

effective at masking speciﬁc regions, as can be clearly seen for the NE quadrant in

the middle row. The global inﬂuence of masking can be estimated from the power

spectrum of the Gauss coefﬁcients at the CMB (Supplementary Fig. 7). This shows

little inﬂuence from the masking procedure, which is conﬁrmed by visual inspection of global plots of the CMB ﬁeld (not shown).

Data availability
All data produced for this manuscript is available from the authors on reasonable request.

Code availability
The dynamo code used to produce the simulations has been described and benchmarked in ref. 45. The postprocessing steps that produce grids of magnetic elements and VGPs are standard and are described in refs. 46,13. The code used to produce rates of change is described in ref. 13 and is available from the authors on reasonable request.

Received: 6 September 2019; Accepted: 26 May 2020;

References
1. Finlay, C., Olsen, N., Kotsiaros, S., Gillet, N. & Tøffner-Clausen, L. Recent geomagnetic secular variation from Swarm and ground observatories as estimated in the CHAOS-6 geomagnetic ﬁeld model. Earth Planets Space 68, 112 (2016).
2. Jackson, A., Jonkers, A. & Walker, M. Four centuries of geomagnetic secular variation from historical records. Phil. Trans. R. Soc. Lond. A 358, 957–990 (2000).
3. Gubbins, D. & Gibbons, S. J. In Timescales of the Paleomagnetic Field, Geophysical Monograph Series Vol. 145, 279–286 (eds Channell, J., Kent, D., Lowrie, W. & Meert, J.) (American Geophysical Union, 2004).
4. Aubert, J., Finlay, C. & Fournier, A. Bottom-up control of geomagnetic secular variation by the Earthas inner core. Nature 502, 219–223 (2013).
5. Livermore, P., Hollerbach, R. & Finlay, C. An accelerating high-latitude jet in Earth’s core. Nat. Geosci. 10, 62 (2017).

6. Ben-Yosef, E. et al. Geomagnetic intensity spike recorded in high resolution slag deposit in Southern Jordan. Earth Planet. Sci. Lett. 287, 529–539 (2009).
7. Shaar, R. et al. Large geomagnetic ﬁeld anomalies revealed in Bronze to Iron Age archeomagnetic data from Tel Megiddo and Tel Hazor, Israel. Earth Planet. Sci. Lett. 442, 173–185 (2016).
8. Ben-Yosef, E., Millman, M., Shaar, R., Tauxe, L. & Lipschits, O. Six centuries of geomagnetic intensity variations recorded by royal Judean stamped jar handles. Proc. Natl. Acad. Sci. USA 114, 2160–2165 (2017).
9. Thébault, E. et al. International geomagnetic reference ﬁeld: the 12th generation. Earth Planets Space 67, 79 (2015).
10. Constable, C., Korte, M. & Panovska, S. Persistent high paleosecular variation activity in southern hemisphere for at least 10 000 years. Earth Planet. Sci. Lett. 453, 78–86 (2016).
11. Fournier, A., Gallet, Y., Usoskin, I., Livermore, P. W. & Kovaltsov, G. A. The impact of geomagnetic spikes on the production rates of cosmogenic 14C and 10Be in the Earth’s atmosphere. Geophys. Res. Lett. 42, 2759–2766 (2015).
12. Korte, M. & Constable, C. Archeomagnetic intensity spikes: Global or regional geomagnetic ﬁeld features? Front. Earth Sci. 6, 17 (2018).
13. Davies, C. & Constable, C. Searching for geomagnetic spikes in numerical dynamo simulations. Earth Planet. Sci. Lett. 504, 72–83 (2018).
14. Livermore, P. W., Fournier, A. & Gallet, Y. Core-ﬂow constraints on extreme archeomagnetic intensity changes. Earth Planet. Sci. Lett. 387, 145–156 (2014).
15. Davies, C. & Constable, C. Geomagnetic spikes on the core-mantle boundary. Nat. Commun. 8, 15593 (2017).
16. Cai, S. et al. Archaeointensity results spanning the past 6 kiloyears from eastern China and implications for extreme behaviors of the geomagnetic ﬁeld. Proc. Natl. Acad. Sci. USA 114, 39–44 (2017).
17. Bourne, M. D. et al. High-intensity geomagnetic ﬁeld ‘spike’ observed at ca. 3000 cal BP in Texas, USA. Earth Planet. Sci. Lett. 442, 80–92 (2016).
18. Valet, J.-P. & Fournier, A. Deciphering records of geomagnetic reversals. Rev. Geophys. 54, 410–446 (2016).
19. Coe, R., Jarboe, N., LeGoff, M. & Petersen, N. Demise of the rapid-ﬁeldchange hypothesis at Steens mountain: the crucial role of continuous thermal demagnetization. Earth Planet. Sci. Lett. 400, 302–312 (2014).
20. Sagnotti, L. et al. Extremely rapid directional change during MatuyamaBrunhes geomagnetic polarity reversal. Geophys. J. Int. 199, 1110–1124 (2014).
21. Evans, M. & Muxworthy, A. A re-appraisal of the proposed rapid MatuyamaBrunhes geomagnetic reversal in the Sulmona Basin, Italy. Geophys. J. Int. 213, 1744–1750 (2018).
22. Sagnotti, L. et al. On the reliability of the Matuyama-Brunhes record in the Sulmona Basin–Comment to ‘A reappraisal of the proposed rapid Matuyama–Brunhes geomagnetic reversal in the Sulmona Basin, Italy’ by Evans and Muxworthy (2018). Geophys. J. Int. 216, 296–301 (2018).
23. Brown, M. C., Holme, R. & Bargery, A. Exploring the inﬂuence of the nondipole ﬁeld on magnetic records for ﬁeld reversals and excursions. Geophys. J. Int. 168, 541–550 (2007).
24. Korte, M. & Constable, C. Continuous geomagnetic ﬁeld models for the past 7 millennia: 2. CALS7K. Geochem. Geophys. Geosys 6, 1–18 (2005).
25. Panovska, S., Constable, C. & Korte, M. Extending global continuous geomagnetic ﬁeld reconstructions on timescales beyond human civilization. Geochem. Geophys. Geosys. 19, 4757–4772 (2018).
26. Christensen, U., Aubert, J. & Hulot, G. Conditions for Earth-like geodynamo models. Earth Planet. Sci. Lett. 296, 487–496 (2010).
27. Davies, C. & Constable, C. Insights from geodynamo simulations into long-term geomagnetic ﬁeld behaviour. Earth Planet. Sci. Lett. 404, 238–249 (2014).
28. Korte, M., Brown, M., Panovska, S. & Wardinski, I. Robust characteristics of the Laschamp and Mono Lake geomagnetic excursions: results from global ﬁeld models. Front. Earth Sci. 7, 86 (2019).
29. Amit, H., Korte, M., Aubert, J., Constable, C. & Hulot, G. The timedependence of intense archeomagnetic ﬂux patches. J. Geophys. Res. 116, B12106 (2011).
30. Brown, M., Korte, M., Holme, R., Wardinski, I. & Gunnarson, S. Earth’s magnetic ﬁeld is probably not reversing. Proc. Natl. Acad. Sci. USA 115, 5111–5116 (2018).
31. Constable, C. & Korte, M. In Treatise on Geophysics Vol. 5, 309–341 (ed. Schubert, G.) (Elsevier, 2015).
32. Helffrich, G. & Kaneshima, S. Outer-core compositional stratiﬁcation from observed core wave speed proﬁles. Nature 468, 807–809 (2010).
33. Lister, J. & Buffett, B. Stratiﬁcation of the outer core at the core-mantle boundary. Phys. Earth Planet. Int. 105, 5–19 (1998).
34. Buffett, B. & Seagle, C. Stratiﬁcation of the top of the core due to chemical interactions with the mantle. J. Geophys. Res. 115, B04407 (2010).
35. Davies, C., Pozzo, M., Gubbins, D. & Alfè, D. Partitioning of oxygen between ferropericlase and Earth’s liquid core. Geophys. Res. Lett. 45, 6042–6050 (2018).

NATURE COMMUNICATIONS | (2020)11:3371 | https://doi.org/10.1038/s41467-020-16888-0 | www.nature.com/naturecommunications

9

ARTICLE

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16888-0

36. Christensen, U. & Aubert, J. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic ﬁelds. Geophys. J. Int. 166, 97–114 (2006).
37. Gubbins, D., Willis, A. & Sreenivasan, B. Correlation of Earth’s magnetic ﬁeld with lower mantle thermal and seismic structure. Phys. Earth Planet. Int. 162, 256–260 (2007).
38. Christensen, U. & Olson, P. Secular variation in numerical geodynamo models with lateral variations of boundary heat ﬂow. Phys. Earth Planet. Int. 138, 39–54 (2003).
39. Mound, J., Davies, C., Rost, S. & Aurnou, J. Regional stratiﬁcation at the top of Earth’s core due to core-mantle boundary heat ﬂux variations. Nat. Geosci. 12, 575–580 (2019).
40. Masters, G., Laske, G., Bolton, H. & Dziewonski, A. In Earth’s Deep Interior: Mineral Physics and Tomography from the Atomic to the Global Scale, Geophysical Monograph Series 63–87 (eds Karato, S., Forte, A., Liebermann, R., Masters, G. & Stixrude, L.) (American Geophysical Union, 2000).
41. Mound, J., Davies, C. & Silva, L. Inner core translation and the hemispheric balance of the geomagnetic ﬁeld. Earth Planet. Sci. Lett. 424, 148–157 (2015).
42. Sprain, C. J., Biggin, A. J., Davies, C. J., Bono, R. K. & Meduri, D. G. An assessment of long duration geodynamo simulations using new paleomagnetic modeling criteria (QPM). Earth Planet. Sci. Lett. 526, 115758 (2019).
43. Merrill R. T., McElhinny, M. W. & McFadden, P. L. The Magnetic Field of the Earth—Paleomagnetism, the Core, and the Deep Mantle, International Geophysics Series Vol. 63 (Academic Press, 1996).
44. Wieczorek, M. A. & Meschede, M. Shtools: tools for working with spherical harmonics. Geochem. Geophys. Geosys. 19, 2574–2592 (2018).
45. Willis, A., Sreenivasan, B. & Gubbins, D. Thermal core-mantle interaction: exploring regimes for ‘locked’ dynamo action. Phys. Earth Planet. Int. 165, 83–92 (2007).
46. Davies, C., Gubbins, D., Willis, A. & Jimack, P. Time-averaged paleomagnetic ﬁeld and secular variation: Predictions from dynamo solutions based on lower mantle seismic tomography. Phys. Earth Planet. Int. 169, 194–203 (2008).
Acknowledgements
C.J.D. acknowledges a Natural Environment Research Council personal fellowship, reference NE/L011328/1. This work was also supported by US National Science Foundation grant EAR1623786 for C.G.C.

Author contributions
C.J.D. developed the ideas, wrote and ran the computer codes, analysed the results and wrote the manuscript. C.G.C. developed the ideas, analysed the results and wrote the manuscript.
Competing interests
The authors declare no competing interests.
Additional information
Supplementary information is available for this paper at https://doi.org/10.1038/s41467020-16888-0.
Correspondence and requests for materials should be addressed to C.J.D.
Peer review information Nature Communications thanks Paul A. Bedrosian, Dominique Jault and the other anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.
Reprints and permission information is available at http://www.nature.com/reprints
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional afﬁliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.
© The Author(s) 2020

10

NATURE COMMUNICATIONS | (2020)11:3371 | https://doi.org/10.1038/s41467-020-16888-0 | www.nature.com/naturecommunications