# Core-flow constraints on extreme archeomagnetic intensity changes

Philip W. Livermore, Alexandre Fournier, Yves Gallet

► To cite this version:

Philip W. Livermore, Alexandre Fournier, Yves Gallet. Core-flow constraints on extreme archeomagnetic intensity changes. Earth and Planetary Science Letters, 2014, 387, pp.145-156. 10.1016/j.epsl.2013.11.020. insu-03581150

HAL Id: insu-03581150

https://insu.hal.science/insu-03581150

Submitted on 18 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Core-flow constraints on extreme archeomagnetic intensity changes

Philip W. Livermore<sup>a,*</sup>, Alexandre Fournier<sup>b</sup>, Yves Gallet<sup>b</sup>

<figure><img src="image_2.png" /></figure>

CrossMark

## ARTICLE INFO

**Article history:**

Received 23 July 2013

Received in revised form 8 November 2013

Accepted 11 November 2013

Available online 7 December 2013

Editor: L. Stixrude

**Keywords:**

archeomagnetism

intensity spike

geodynamo

core-flow

upper bounds

## ABSTRACT

Recent studies (Ben-Yosef et al., 2009; Shaar et al., 2011) propose extreme archeomagnetic intensity changes (termed spikes) in the range ~4–5 µT/year c.a. 1000 BC in the Near East, around 40 to 50 times larger than values typical of the present-day. In order to investigate whether such extreme changes are consistent with a model of the source region of the magnetic field, namely the fluid flow at the surface of Earth's core, we construct upper bounds for instantaneous magnetic intensity change at an arbitrary site on the Earth's surface. These bounds are constrained by the amount of kinetic energy available to sustain the change, taken here to be a prescribed value for the root-mean-squared surface velocity of 13 km/yr as inferred from the current state of the core. Further, we focus attention on two end-members of optimised core surface flow structure: unrestricted and purely-toroidal. As the derivation of the bounds demands complete knowledge of the geomagnetic field at the core surface, we model the unknown field by means of a Monte Carlo approach, extending to high degree the CHAOS-4 (epoch 2010 AD) and CALS10k1.b (epoch 1000 BC) geomagnetic field models.

Using 2000 realisations for each family of stochastic field models, we find that optimised core flows are always large-scale and that they tend to generate a non-dipole, quadrupole-dominated secular variation at the Earth's surface. The dependence of the upper bounds as a function of site location reflects the large-scale structure of the intensity itself: stronger field permits more rapid change. For the site in the Near East, purely-toroidal flows have upper bounds of approximately 0.62 ± 0.02 µT/year, whereas unrestricted flows increase this bound to 1.20 ± 0.02 µT/year. We favour the former as more geophysically sound, on the account of a large body of previous results from core surface flow inversions and consistency with the existence of a stratified layer at the top of the core. Even if we allow for a generous threefold increase in the prescribed rms velocity (and a concomitant threefold increase in the bound), we conclude that the reported occurrences of extreme intensity changes as suggested in the Near East are not compatible with the commonly accepted structure of core-surface flow. However, it may be that an explanation for spikes lies beyond our current perception of core-dynamics and future work would be further motivated by seeking corroborative evidence of rapid intensity change from sites elsewhere on Earth's surface; we therefore also discuss the form that the secular variation would take in the case of simultanenous archeomagnetic spikes.

© 2013 The Authors. Published by Elsevier B.V. Open access under CC BY license.

## 1. Introduction

Over the last few years there has been considerable effort focused on studies of geomagnetic field intensity (F) variations over the past several millennia, mainly analysing the thermoremanent magnetization carried by archeological artifacts heated at the time of their manufacture or use. These studies led to the emergence of increasingly detailed composite intensity variation curves for European regions (e.g. De Marco et al., 2008; Kovacheva et al., 2009; Gómez-Paccard et al., 2012; Genevey et al., 2013; Hervé et al.,

2013), as well as for the Middle East, where many opportunities of sampling recently allowed the recovery of the geomagnetic intensity fluctuations over the past millennia BC (e.g., Thébault and Gallet, 2010; Gallet et al., in press). The results highlighted a series of regional centennial-scale intensity maxima. In particular, Genevey et al. (2009, 2013) showed the existence in Western Europe of three intensity peaks during the last millennium (during the 12th century, the second half of the fourteenth century and around 1600 AD) that the reconstructed time-varying global archeomagnetic field models are unable to capture, because of their still limited temporal resolution (Korte et al., 2011; Licht et al., 2013). The peaks observed in Western Europe are associated with intensity variation rates of $dF/dt \sim 5–10$ µT per century (or 50–100 nT/yr), comparable to variation rates presently

<footnote><sup>a</sup> School of Earth and Environment, University of Leeds, Leeds, LS2 9JT, UK</footnote>

<footnote><sup>b</sup> Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Université Paris Diderot, UMR 7154 CNRS, F-75005 Paris, France</footnote>

<footnote>* Corresponding author.
E-mail address: p.w.livermore@leeds.ac.uk (P.W. Livermore).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Contours of (a) *F* and (b) *dF/dt* from CHAOS-4 at epoch 2010 to degree 13 on the Earth's surface. The maximum rate of change of intensity is currently ~100 nT/yr.</figcaption></figure>

prevailing in many regions (as shown in Fig. 1). This rate from Genevey et al. (2009, 2013) is determined from an averaged curve and the existence of episodes with larger intensity variation rates cannot be excluded. Such a possibility was proposed by Gómez-Paccard et al. (2012) for Western Europe from rough comparisons between pairs of intensity data, although the derived rates did not take into account uncertainties in the intensity or age. These authors suggested variation rates larger than ~10 µT per century (100 nT/yr) during the Medieval era with values up to 50–80 µT per century (500–800 nT/yr) — such high rates obviously not persisting for the entire period.<fnref n="1" /> While variations close to ~100 nT/yr presently exist only in the Indian ocean and the Americas, values larger than ~150 nT/yr appear very unusual.

The best evidence for very rapid intensity variations at a rate much larger than ~100 nT/yr was most probably provided by archeointensity data obtained from the Near East. Recent archeomagnetic studies conducted in Anatolian, Syrian and Levantine regions made it possible to assemble the main intensity variation patterns over the last three millennia BC: see Fig. 2(a); (Genevey et al., 2003; Gallet and Le Goff, 2006; Gallet et al., 2006, 2008, in press; Gallet and Al-Maqdissi, 2010; Ben-Yosef et al., 2008, 2009; Shaar et al., 2011; Ertepinar et al., 2012). These changes are marked by several relative intensity maxima, among which three occurred during the Bronze Age (~3000–1200 BC), with variation rates again of ~50–100 nT/yr (e.g. Thébault and Gallet, 2010; Gallet et al., in press). But the latter appear very minor compared to the rates reported for the very beginning of the first millennium BC (beginning of the Iron Age; Fig. 2(b)). These data were principally obtained from paleomagnetic analyses of metalurgical residues (slag) in the Levant (Ben-Yosef et al., 2008, 2009; Shaar et al., 2011), and more recently from baked-clay fragments collected from a kiln discovered at Arslantepe (Turkey; Ertepinar et al., 2012). Furthermore, a selection of the available data over a wider geographic area (see Fig. 1 in Shaar et al., 2011, with data selected between 20° N and 50° N in latitude and between 10° E and 50° E in longitude) does not change the main message conveyed by Fig. 2: the beginning of the first millennium BC was characterised in the Near East (and elsewhere) by the highest geomagnetic field intensities ever observed during the Holocene and beyond. There was a suggestion of this phenomenon in the analysis of global and regional compilations of geomagnetic field intensity data (e.g. Yang et al., 2000; Genevey et al., 2008; Knudsen et al., 2008), but not with such a high intensity level. The key information obtained by Ben-Yosef et al. (2009) and Shaar et al. (2011) is that these extreme geomagnetic intensities, larger than 100 µT, were reached in a very short time in-

terval (~10–20 years), which led to the concept of a “geomagnetic intensity spike” (Ben-Yosef et al., 2009). Two such events have been proposed so far, at ~980 BC and ~890 BC (Fig. 2(b); from Fig. 7 by Shaar et al., 2011). Fig. 2(b) shows that regional intensity variations as extreme as ~30 µT (Shaar et al., 2011, considering solely the data from archeological site Timna-30), even possibly ~50 µT (Ben-Yosef et al., 2009, further considering the data from site Khirbat en-Nahas, abbreviated in what follows to “KEN”), beyond the previously assumed intensity level (~75 µT) would have occurred in approximately a decade. Figs. 2(c) and 2(d) show smooth fits<fnref n="2" /> through the data (orange curve: Timna-30 and KEN datasets, red: Timna-30 only); we present separate analyses as only the intensity recorded at Timna-30 is independently verified by a measurement of comparable magnitude nearby (Ertepinar et al., 2012). In both cases Fig. 2(d) shows similar rates of change of intensity of up to 4–5 µT/yr during the two proposed spike-events.

These archeomagnetic spikes, which have not been found elsewhere, in particular in eastern Europe (for instance neither in the Bulgarian nor Greek data sets; De Marco et al., 2008; Kovacheva et al., 2009; Tema and Kondopoulou, 2011), pose a number of questions, in particular on the precise geographical extension and time-scale of the spike events, on the relative contribution of the dipolar and non-dipolar field components to these features, as well as on the reliability of the extreme intensity data, although the latter have met strict quality criteria, or on the dating precision of these data. The kiln studied by Ertepinar et al. (2012) at Arslantepe yielded a spike-like intensity value (~100 µT), thereby confirming the extremely high intensity level during the early Iron Age. However, Ertepinar et al. (2012) emphasized the fact that the dating uncertainties for this kiln are 300 yr. It may be that the date at which the oven was fired happened to coincide with the proposed 20–30 yr spike, although this may be regarded as too much of a coincidence. Should the dates not coincide, this would favour not a spike but a centennial-scale intensity peak, albeit extreme, like other maxima so far observed in the Near East and in Europe.

Considering the importance of intensity spikes, and more generally of episodes characterised by rapid intensity changes for deciphering the geomagnetic field behaviour during the Holocene, we decided in the present study to approach the topic from a different perspective: that of asking what intensity variation rates in the

<footnote n="1">¹ Although we have reported intensity variation rates expressed per century to be consistent with the cited literature, in view of the short time-scales under consideration in this study, hereafter we will only refer to rates expressed per year.</footnote>

<footnote n="2">² We generated 100 realisations (thin blue lines) of data bootstrapped from Shaar et al. (2011) and Ben-Yosef et al. (2009), based on a uniform distribution of ages (within the bounds given) and normally distributed intensity (with the mean and standard deviation as given). We fitted uniformly-weighted cubic splines through the data using the Matlab routine SPAPS, where the tradeoff between the fit to the data and the smoothness was chosen so that the maximum $|dF/dt|$ was at most 10 µT/yr, a factor of 100 higher than typical archeomagnetic intensity change. The solid orange line shows the average spline, the red shows the equivalent average restricted to the Timna-30 dataset.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Centennial to millennial-scale geomagnetic field intensity variations in the Near East as recovered from recent archeointensity data sets. (a) Temporal distribution of archeointensity data covering the last three millennia BC obtained from Turkey (red circles, Ertepinar et al., 2012), Syria (blue circles, Genevey et al., 2003; Gallet and Le Goff, 2006; Gallet et al., 2006, 2008, in press; Gallet and Al-Maqdissi, 2010) and from the Levantine region (green squares, Ben-Yosef et al., 2008; green open circles, Ben-Yosef et al., 2009) with age calibration according to Shaar et al., 2011; green solid circles, Shaar et al., 2011). All data were reduced to the latitude of Ebla/Tell Mardikh (~35.8° N). (b) Composite high resolution geomagnetic field intensity record between 1050 and 850 BC reconstructed by Shaar et al. (2011) from archeointensity results obtained from two early Iron Age archeological sites, Kharbat en-Nahas (KEN) (open circles, Ben-Yosef et al., 2009) and Timna-30 (closed circles, Shaar et al., 2011). All data were reduced to the latitude of site Timna-30 (~29.8° N). The age calibration is provided by ¹⁴C dates and stratigraphic constraints. Two spikes in geomagnetic intensity were proposed by Shaar et al. (2011) around 980 BC and 890 BC (blue shaded zones). (c) as (b) but over-plotted with cubic-spline fits to 100 realisations of this dataset (thin blue lines), each constrained by a maximum intensity change of 10 µT/yr; orange shows the average curve, the red is the equivalent average restricted to the Timna-30 dataset. (d) The values of dF/dt for each curve in (c) as a function of time. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Near East (specifically at the site Timna-30 used in Shaar et al., 2011) are consistent with the process of magnetic field generation in Earth's core.

We adopt an optimisation approach, and ask, “what is the largest rate of change of intensity, $dF/dt$, at a site on the Earth’s surface consistent with a reasonable description of the physics of magnetic field generation?”. The physics to which we refer is that of magnetic induction on the core–mantle-boundary (CMB), the principle ingredient of which on short time-scales is the action of the core-flow, advecting and shearing magnetic field-lines. We constrain the flow amplitude to be broadly consistent with our current understanding of the geodynamo, which we take to be the condition that the root-mean-squared flow speed on the CMB is 13 km/yr (Holme, 2007). Such a constraint is vital since, as faster flows generate more rapid secular variation, the maximum value of $dF/dt$ would be formally unbounded if the flow was unconstrained (i.e. allowed to become infinitely large). Subject to this rms constraint we optimise over all possible configurations of flow on the CMB, encompassing every possible state of the geodynamo consistent with this assumption. Within this well-defined framework, the upper bounds on $dF/dt$ will allow us to critically assess the likelihood of the suggested extreme intensity changes.

We remark briefly that an alternative method is to run geodynamo models, recording the intensity changes with time and noting the largest and fastest variations. These models (see e.g. Christensen and Wicht, 2007, for a review) are now arguably able

to reproduce the main features of the observed geomagnetic field, provided the time-scales of rotation, advection by fluid flow, and magnetic diffusion, are in approximate geophysical proportions (Christensen et al., 2010). However, although the models are expected to provide a reasonably good approximation of the secular variation on time-scales of the order of a century and longer, they cannot resolve the interannual to decadal time-scales, the very features of interest in this study (consult Finlay et al., 2010, for further reading on this topic).

A full description of our method is laid out in Section 2, following which we provide a pedagogical example in Section 3. In Section 4 we describe our Monte Carlo models of magnetic field, and in Section 5 we present our main results with concluding remarks in Section 6. We end with a discussion in Section 7.

## 2. Methodology

We consider an observation site S on the Earth's surface (radius $a = 6371$ km), at colatitude and longitude ($\theta_s, \phi_s$). Changes in this magnetic field structure everywhere outside the core (assuming an electrically insulating mantle), including at S, are slaved to changes in the magnetic field on the CMB (Jackson and Finlay, 2007). Thus, for the purposes of constructing upper bounds we do not require a complete model of magnetic field generation inside the core, but only a description of its changes on the edge of the region where the geodynamo operates.
The optimisation methodology that we describe below requires complete knowledge of the magnetic field, **B**, on the CMB: this point is discussed further within the framework of Monte Carlo methods in Section 4. Further, we neglect radial diffusion of the magnetic field on the CMB as it would also require knowledge of gradients of **B** which are observationally unconstrained; we discuss this limitation in Section 7. We begin by describing the method in the limit of frozen-flux (Roberts and Scott, 1965; Jackson and Finlay, 2007), in which magnetic diffusion is neglected on the CMB; adding in horizontal diffusion is then a simple extension.

At S, the intensity is given by

$$F = |\mathbf{B}| = \sqrt{B_r^2 + B_{\theta}^2 + B_{\phi}^2}, \qquad (1)$$

so that

$$\frac{dF}{dt} = \frac{1}{F} \left( B_r \frac{dB_r}{dt} + B_{\theta} \frac{dB_{\theta}}{dt} + B_{\phi} \frac{dB_{\phi}}{dt} \right) = \frac{1}{F} \left( \mathbf{B} \cdot \frac{d\mathbf{B}}{dt} \right). \qquad (2)$$

Ignoring diffusion, the radial-component of the instantaneous rate of change of magnetic field on the CMB ($r=c$, $c=3485$ km) is

$$\frac{\partial B_r}{\partial t} = -\nabla_H \cdot (\mathbf{u}_H B_r), \qquad (3)$$

where $\mathbf{u}_H$ denotes the horizontal flow and $\nabla_H = \nabla - \hat{\mathbf{r}} \frac{\partial}{\partial t}$ is the horizontal part of the gradient operator. Modelling **B** and hence $\partial \mathbf{B} / \partial t$ as a potential field in the mantle $r \ge c$, we adopt the usual representation

$$\frac{\partial \mathbf{B}}{\partial t} = -a \nabla \sum_{l,m} \left(\frac{a}{r}\right)^{l+1} Y_l^m(\theta, \phi) \dot{\beta}_l^m \qquad (4)$$

where $\dot{\beta}_l^m = (\dot{g}_l^m, \dot{h}_l^m)$ is the vector of Gauss coefficients of the secular variation, and $Y_l^m$ are Schmidt quasi-normalised spherical harmonics. This equation is simply the description of the upwards-continuation of the magnetic field (or rather its rate of change) from $r=c$ to $r=a$. By (4) and (3), each $\dot{\beta}_l^m$ depends linearly on $\mathbf{u}$, from which it follows that all three components of $\partial \mathbf{B} / \partial t$ at $r=a$, and indeed $\frac{dF}{dt}$, depend linearly on $\mathbf{u}$ at $r=c$. The relation is trivially also homogeneous as if $\mathbf{u}=\mathbf{0}$ there can be no secular variation.

We can exploit this elementary relation to express $dF/dt$ as a linear homogeneous function of the flow,

$$\frac{dF}{dt} = \mathbf{G}^T \mathbf{q}, \qquad (5)$$

for some column vector **G** (which needs to be computed), and where the flow is written as an expansion over a set of modes, $\mathbf{u}_k$, whose coefficients are stored in the vector **q**:

$$\mathbf{u} = \sum_{k} q_k \mathbf{u}_k. \qquad (6)$$

In simple terms, element $G_k$ gives the contribution to the rate of change of intensity at S from flow mode $\mathbf{u}_k$, and encapsulates information about the background field, the coordinates of S and also the structure of the flow mode. All flows we will consider are divergence-free, and within this class we will consider the most general set which we term unrestricted, comprised of both toroidal and poloidal modes, respectively

$$\nabla \times t_l^m Y_l^m(\theta, \phi) \mathbf{r}, \quad \nabla_H [r s_l^m Y_l^m(\theta, \phi)], \qquad (7)$$

up to a fixed spherical harmonic degree $L_U$, where $\mathbf{u}_k$ is defined by the vector of coefficients ($t_l^m, s_l^m$) of zeros except for a one in the

kth position. We will also consider flows which are also purely-toroidal, describing flows that are everywhere horizontal.

For a given site location and prescribed magnetic field (to degree $L_B$), the vector **G** is then straightforward to assemble, one element at a time. Each mode of flow was taken in sequence and we calculated (i) the spherical harmonic spectrum to degree $L_B + L_U$ of $\partial B_r / \partial t$ at $r=c$ using a standard transform methodology based on Gauss-Legendre quadrature and the fast-Fourier transform, and then (ii) using (4) the contribution to $dF/dt$ at S. All elements of **G** are exact to machine precision.

The statement of the problem we need to solve then is rather succinct: maximise $\mathbf{G}^T \mathbf{q}$ subject to constraints on the flow. One feature of core-flows on which inverse studies broadly agree is the root-mean-squared (rms) velocity:

$$T_0 = \sqrt{\frac{1}{4\pi} \int_{r=c} |\mathbf{u}|^2 d\Omega} \qquad (8)$$

where $d\Omega$ is the element of solid angle, for which we take $T_0 = 13$ km/yr (Holme, 2007) as a (typical) target value. In our representation of flow, this constraint may be written

$$\mathbf{q}^T E \mathbf{q} = \frac{1}{4\pi} \sum_l (l(l+1)([t_l^m]^2 + [s_l^m]^2)) \frac{4\pi}{2l+1} = T_0^2 \qquad (9)$$

where E is a diagonal matrix with elements $l(l+1)(2l+1)^{-1}$. Expressing the rms constraint using a Lagrange multiplier $\lambda$, we therefore seek the constrained maximum

$$\left. \frac{dF}{dt} \right|_{\max} = \max_{\mathbf{q}} [\mathbf{G}^T \mathbf{q} - \lambda(\mathbf{q}^T E \mathbf{q} - T_0^2)]. \qquad (10)$$

At a local maximum, the gradient with respect to each component of the vector **q** is zero, giving

$$\mathbf{G} - 2\lambda E \mathbf{q} = \mathbf{0}, \qquad (11)$$

and so

$$\mathbf{q} = \frac{1}{2\lambda} E^{-1} \mathbf{G}, \qquad (12)$$

where $\lambda > 0$ is found by scaling the flow to the target rms. Note that changing the sign of $\lambda$ corresponds to the flow $-\mathbf{u}$ which minimises (rather than maximises) $dF/dt$. Finally, it is worth remarking that although the intensity depends nonlinearly on **B**, here we only need optimise a linear quantity subject to a quadratic constraint, whose (unique) solution is given above. More complex constraints on the flow may not lead to such a straightforward solution.

The inclusion of magnetic diffusion is a simple extension of the above methodology. Eq. (3) now becomes

$$\frac{\partial B_r}{\partial t} = -\nabla_H \cdot (\mathbf{u}_H B_r) + \eta \hat{\mathbf{r}} \cdot \nabla^2 \mathbf{B}, \qquad (13)$$

in which $\eta$ denotes the magnetic diffusivity and $\hat{\mathbf{r}}$ the unit position vector. The last term is independent of the flow and can be split into two parts, stemming from either horizontal or radial diffusion of the magnetic field. We consider only horizontal diffusion, which, at $r=c$, may be written

$$\nabla_H^2 B_r = - \sum_{l,m} \frac{l(l+1)}{c^2} [B_r]_{(l,m)} Y_l^m \qquad (14)$$

where $[B_r]_{(l,m)}$ denotes the $Y_l^m$ component of $B_r$. Since **B** is assumed known, (14) provides an additional term, independent of **u**, to (10) and has no bearing on the optimising flow. We take $\eta = 0.7$ m$^2$/s, based on a core electrical conductivity of
<figure><img src="image_5.png" /><figcaption>(b) Purely-toroidal flows</figcaption></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

**Fig. 3.** Assuming a geomagnetic field of an axial dipole with $g_1^0 = -30000$ nT, contours of $dF/dt$ maximised for an observation site at 52° N, 0° E (labelled by the star) with the optimal flow pattern for (a) unrestricted flows (b) and purely-toroidal flows. Each flow is normalised to have a rms value of 13 km/yr; the local flow speed is shown both by the coloured contours and the arrow size relative to the scale shown.

$\sigma = 1.1 \times 10^6 \Omega^{-1} m^{-1}$ (Pozzo et al., 2013). This value, the smallest such compatible with current theory, provides us with a highest estimate of the upper bounds and is discussed further in Section 4.2.

Finally, it is worth remarking that the geometric attenuation of the magnetic field through the mantle, far from being an obstacle in our analysis, actually makes feasible the optimisation procedure as described above. A theoretical observation site on $r=c$ would have an optimising flow, and associated secular variation, non-zero only at this point and consequently infinite in size. In contrast, sites on the Earth's surface are sensitive primarily to large-scale secular variation (since the smallest scales are strongly attenuated in the mantle), which we find generated by large-scale flows. Thus we anticipate (and find) that the optimisation procedure is well-posed, the optima of $dF/dt$ are finite and the resulting flows large-scale.

**3. A pedagogical example: optimal change in an axial dipole field**

We now illustrate the analysis as described above by considering the idealised example of a purely axial dipolar geomagnetic field that is described by the Gauss coefficient $g_1^0 = -30000$ nT. The flow is expanded to degree $L_U = 25$, which ensures convergence.

Consider an observation point in London (latitude 52° N, longitude 0°), at which we optimise $dF/dt$. Fig. 3 shows a contour map of $dF/dt$ on the Earth's surface along with the optimising flow at the core surface $r=c$, for both (a) unrestricted and (b) purely-toroidal flow cases.

A striking feature of the plots is the localisation of $dF/dt$, focused into quasi-circular contours, taking maximum values of (a) 887 nT/yr and (b) 171 nT/yr at the observation site. Away from the observation site, much of the Earth's surface has small and
negative values of $dF/dt$ showing that the rapid increase in intensity at $S$ comes at the expense of slower decreases elsewhere. Another point of note is that the maximum value of $dF/dt$ does not coincide with the site location. This is not an error, but simply reflects the fact that with the flow optimised for site $S$, the (different) structure of the magnetic field elsewhere is such that it produces a greater value of $dF/dt$ than at $S$. In (a), the background field is stronger at higher latitudes, which explains why the maximum of $dF/dt$ is shifted to the north of London. However, in (b) it is the horizontal gradient of $B_r$ that is important: $dF/dt = -u_H \cdot \nabla_H B_r$ (from Eq. (3)); since the largest horizontal gradients of $B_r$ are at the equator, the high is to the south of London. The values of the upper bounds themselves evaluated at London are really only of pedagogical value, as it turns out that the upper bounds change as we introduce more complex field structures. However, one persistent feature is the relative disparity of the optimum bounds from the two types of flow. Since purely-toroidal flows are a subset of the unrestricted flows, it follows that their associated maxima must be smaller. Yet, this simple statement tells us nothing about the actual difference: here, they vary by a factor of 5; as we will see later, in the more general case a factor of 2 is more typical.

The structure of the flows is very instructive, and illustrative of the more general case. Both flows, and their associated secular variation, are large-scale (see also supplementary Figs. 11(a-b)). The unrestricted flow is dominated by its high horizontal-convergence very close to, but not coincident with, the radial-projection of the observation site. Horizontally-convergent flow (with an associated down-welling) is a very efficient method of concentrating and amplifying magnetic flux, and thus results in a rapid change of intensity. By contrast, because it is everywhere horizontal the purely-toroidal case has to rely on a much less optimal method of creating intensity changes, that of advection of the poloidal field (in this case, due south beneath $S$, by means of a pair of counter-rotating vortices). Lastly, it is interesting that although both flows have rms values of 13 km/yr, they have maximum speeds of approximately 50 km/yr. Although in line with estimates of the maximum speed in the core (Finlay and Amit, 2011), for more complex magnetic fields the associated optimising flows are also more spatially complex and have greater maximum speeds than we see here. The results presented here contain horizontal diffusion, although its contribution to $dF/dt$ turns out to be very small: $-0.18$ nT/yr, around 0.1% of the total.

# 4. Stochastic magnetic fields

## 4.1. Magnetic models

To address the question of intensity variability at the beginning of the first millennium BC, our framework formally requires complete knowledge of the main geomagnetic field at that epoch. This information is unfortunately not available, mainly because of the crustal concealing of the small-scale core field (approximately from harmonic degree 14 and beyond). In addition, our ability to describe accurately the theoretically visible large-scale core field at that epoch is further hampered by the poor geographical and temporal distribution of the existing data sets for this period (e.g., Donadini et al., 2009; Korte et al., 2011; Licht et al., 2013). We therefore study two types of model that attempt to characterise the structure of the field. Model ST_CALS (STochaotic CALS10k1.b) is based upon the CALS10k1.b model (Korte et al., 2011), which provides a set of spherical harmonic coefficients to degree 10 (with standard deviations) over the entire Holocene. We use only the more robustly determined coefficients up to degree 4 and extend to higher degree using a stochastic extrapolation technique as described below. The second model, ST_CHAOS, is based upon the

<figure><img src="image_9.png" /><figcaption>Fig. 4. A comparison of 10 realisations of the Mauersberger-Lowes energy spectrum on the CMB as a function of degree $l$, for the two models ST_CHAOS and ST_CALS. The extrapolation is based on the expected spectra of the form $A(2l + 1)^{-1}$.</figcaption></figure>

geomagnetic model CHAOS-4 (Olsen et al., 2010), providing the magnetic field at epoch 2010 (see Fig. 1). While the magnetic field has undoubtedly changed in the intervening 3010 years, the basic structure is much the same: a dominant dipole and high-latitude flux lobes (Korte et al., 2011). We exploit the fact that CHAOS-4 is known to a much greater accuracy since it is constructed from high-quality satellite data of almost uniform coverage, and use it as a tightly constrained proxy for the magnetic field in 1000 BC. In this model, we assume that the magnetic field is given exactly to degree 13, above which crustal contamination forces us to adopt a stochastic extrapolation.

We assume that the unknown Gauss coefficients of degrees 5 to $L_B$ in ST_CALS and degrees 14 to $L_B$ in ST_CHAOS are distributed normally with mean zero and degree-dependent standard deviation $\sigma_l$. Our treatment of the low degrees is model-dependent: in ST_CHAOS, we assume that degrees 1–13 are known exactly, whereas in ST_CALS degrees 1–4 are treated as random variables, with mean and standard deviation taken from CALS10k1.b. Choosing $\sigma_l$ appropriately, and by generating pseudo-random values from this distribution, we then are able to construct multiple Monte Carlo realisations of the magnetic field. We statistically analyse all the optimised solutions, each calculated using a single realisation. We determine $\sigma_l$ by using the fact that the Mauersberger-Lowes energy spectrum (Backus et al., 1996) of the magnetic field, where $R_l$ denotes the contribution from degree $l$ to $\int B^2 d\Omega$ over the CMB, has the expected value

$$ 
\begin{aligned}
E(R_l) &= E \left[ \left( \frac{a}{c} \right)^{2l+4} (l+1) \sum_m \left[ (g_l^m)^2 + (h_l^m)^2 \right] \right] \\
&= \left( \frac{a}{c} \right)^{2l+4} (l+1)(2l+1)\sigma_l^2.
\end{aligned} 
\quad (15) $$

We then assume that this follows the profile $A(2l + 1)^{-1}$ (McLeod, 1996), where $A$ is calculated by a least-squares fit to the nondipolar spectra of the base models (i.e. to the spectrum of degrees 2–13 for ST_CHAOS and the expected spectrum for degrees 2–4 for ST_CALS). Fig. 4 illustrates the magnetic-energy spectra of ten realisations from each family of stochastic models. A significant point of note is that although the dipole component is stronger in CALS10k1.b than in CHAOS-4, degrees 2–4 are weaker, a possible consequence of the regularization imposed on the archeomagnetic field models, owing to the limited spatial coverage of the data sets. The lower energy available in the non-dipole spectrum has an important influence on the maximum $dF/dt$, as we shall see below.
<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Dependence of upper bounds (at Timna-30) on the truncation of the magnetic field extrapolated from ST_CHAOS. In green are shown the first 13 (unique) bounds derived from the (assumed exact) CHAOS-4 geomagnetic field model to degree 13; 2000 Monte Carlo realisations (red dots) are shown for magnetic fields stochastically extrapolated to degrees 20, 50 and 135. In each case, the mean and unbiased standard deviation are marked by the blue square with black error bars. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 4.2. Truncation

In both stochastic models for the magnetic field, we built up a suite of computations using ever larger truncations $L_B$ on the maximum spherical harmonic degree. From this set of upper bounds, how do we choose which is most relevant for the Earth? There are two issues here. Firstly, the spectral form $A(2l + 1)^{-1}$ cannot describe the geodynamo for arbitrarily high degree as the associated total (expected) magnetic energy is divergent (i.e. $\sum_l R_l$ is infinite), so there must be a threshold beyond which the spectra decreases more rapidly. Secondly, as $L_B$ increases, the upper bounds of $dF/dt$ at the site increase because ever smaller scales (and ever larger gradients) in the magnetic field make it easier for a flow to create secular variation. This behaviour is shown in Fig. 5, for a variety of truncations $L_B$ (20, 50, and 135). As we populate the concealed part of the geomagnetic spectrum, we observe a monotonic (albeit slow, quasi-logarithmic) increase in the maximum attainable $dF/dt$.

We exploit the fact that the upper bounds increase with $L_B$ by choosing the maximum truncation that is consistent with the duration of the archeomagnetic spikes of ~30 yr (Shear et al., 2011) – that is, we assume that the magnetic field structure does not diffuse within the lifetime of the spike. This maximum $L_B$ then gives us a highest estimate of the upper bounds. If the small-scale magnetic field changes only through horizontal diffusion, we constrain the e-folding time for magnetic field harmonics of degree $l$, by $c^2(l(l+1)\eta)^{-1} \le 30$ yr. Taking $\eta = 0.7 \text{ m}^2/\text{s}$, this condition is satisfied for every $l \le L_B = 135$. Note that larger values of $\eta$ (consistent with historical, rather than very recently published, values of electrical conductivity of the core) give smaller values of $L_B$ and so $L_B = 135$ provides a “greatest” upper bound over all past, present, and hopefully future, estimates of core conductivity. In addition, the slow, logarithmic increase displayed in Fig. 5 is such that multiplying this truncation by a factor of (say) 2 is likely to induce only a minor change in the estimated upper bounds.

Having noted above that the magnetic field diverges as $L_B$ increases, we check that $L_B = 135$ still gives geophysically reasonable values for the expected field. Associated with this truncation is the

<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Normalised histograms giving the discrete probability density function of the upper bound of $dF/dt$ at Timna-30 with models ST_CHAOS and ST_CALS, for both unrestricted and purely-toroidal flows. In all cases 2000 realisations of the magnetic field extrapolated to degree 135 are used. The purely-toroidal assumption leads to a considerably lower mean upper bound compared to unrestricted flows: 620 nT/yr compared with 1200 nT/yr for ST_CHAOS; 460 nT/yr compared with 950 nT/yr for ST_CALS. The bounds derived from ST_CALS have significantly more variability than those derived from ST_CHAOS, with unbiased standard deviation of 46 nT/yr compared to 17 nT/yr for unrestricted flows, and 88 nT/yr compared to 16 nT/yr for purely-toroidal flows.</figcaption></figure>

(expected) rms value of $B_r = 0.39$ mT and an (expected) minimum Ohmic dissipation of 1.5 GW (see Jackson and Livermore, 2009), all within geophysical bounds. Lastly, we adopt a truncation for the flow of $L_U = L_B + 10$ which provides well-converged solutions, as the optimising flows we find are always predominantly large-scale (see below).

## 5. Results

Fig. 6 shows normalised histograms of upper bounds derived at the site Timna-30 using 2000 realisations of the magnetic field, for both families of field models (ST_CALS and ST_CHAOS) and for unrestricted and purely-toroidal flows. The bounds derived from ST_CHAOS are displaced to the right of those from ST_CALS by 150–300 nT/yr, owing to the stronger non-dipolar magnetic field. Thus the bounds from ST_CHAOS provide a “greatest” bound of $dF/dt$ over both models. For both ST_CHAOS and ST_CALS, the mean purely-toroidal flow bounds are about half those using unrestricted flows (620 nT/yr compared with 1200 nT/yr for ST_CHAOS; 460 nT/yr compared with 950 nT/yr for ST_CALS). These bounds are greater than (and therefore consistent with) the largest values of intensity change typical of the present-day (~0.1 µT/yr), although they are well below those suggested for the archeomagnetic spikes of 4–5 µT/yr.

Another feature of Fig. 6 is the relatively broad peaks of the bounds from ST_CALS: this stems from the fact that all the model coefficients are randomised, whereas those in ST_CHAOS are random only for degrees 14 and above. The narrow distribution of bounds for ST_CHAOS shows that, as already illustrated by Fig. 5, degrees 1 to 13 are the main contributors to an optimised $dF/dt$. In all cases, the contribution from diffusion is small, supplying approximately –0.3 nT/yr to the optima. Note that from additional analysis (supplementary Fig. 12) we concluded that 2000 realisations sufficed to produce converged statistics.

Figs. 7(a–b) show velocity power spectra of the flow for both classes of magnetic models. In each plot, 10 random representative realisations associated with purely-toroidal and unrestricted flows are shown. Figs. 7(c–d) show the energy spectrum of the associated secular variation; the inset in (c) shows a zoomed version of the figure at very low degree. In all cases, the optimal flows are
<figure><img src="image_12.png" /><figcaption>(b) ST_CALS: velocity spectra</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) ST_CHAOS: SV at CMB</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) ST_CHAOS: SV at Earth's surface</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 7.** Velocity power spectra of both unrestricted and purely-toroidal flows for 10 realisations of (a) ST_CHAOS and (b) ST_CALS magnetic models, for which $dF/dt$ is optimised for Timna-30; both the energy per degree $l$ and cumulative energy (up to degree $l$) are shown. Parts (c–d) show the energy spectrum for ST_CHAOS of the secular variation at the CMB and Earth’s surface; the inset in (c) shows the data of the main plot at very low degree.</figcaption></figure>

large-scale (converging by approximately degree 10). This observation justifies a *posteriori* the choice of 13 km/yr as a consistent constraint for the rms flow speed, a value which is based on core-flow inversions which themselves assume a large-scale flow. The optimal secular variation at the core surface increases with degree $l$ until around degree 100 where it reaches a plateau, before an exponential decrease around degree $L_B = 135$. Due to attenuation in the mantle, only the large-scales dominate the spectrum at the Earth’s surface which is around 4 times higher in the unrestricted flow case than the purely-toroidal flows case, giving the factor of two in the secular variation, or rather, $dF/dt$, itself. Interestingly the secular variation is quadrupole-dominated ($l=2$), as is the current geomagnetic field (Hulot et al., 2007).

For the observation site of Timna-30, Fig. 8 shows contours of $dF/dt$ (in both global and regional views), along with the structure of the optimising flow, for the median realisation of the ST_CHAOS model for both unrestricted and purely-toroidal flows. As seen previously in Section 3, the location of the maximum value of $dF/dt$ is close to $S$ but does not coincide exactly. Although the optimising flows were normalised to have a rms value of 13 km/yr, both are strong only in the vicinity of $S$ and have maximum speeds of approximately 130 km/yr. It is also clear just how specialised the optimal flow is to the specific structure of the magnetic field, and that although still large-scale (according to Fig. 7), it differs from the typical core-flows obtained from inversion (e.g. Holme, 2007) by the notable lack of global-scale features. On a related point, Shaar et al. (2011) remarked that it was unclear whether dipolar (by which they mean “global”) or non-dipole field (i.e. “localised”) change is responsible for the archeomagnetic intensity spikes. In all our optimised cases, the intensity change is strongly localised (and therefore “non-dipolar”), although formally since magnetic intensity is a nonlinear quantity, such a separation into dipolar and non-dipolar components is not possible.

It is also of interest to compute comparable upper bounds on $dF/dt$ for other sites. Fig. 9 shows the mean of the upper bounds for the model ST_CHAOS, for both purely-toroidal and unconstrained flows, as a function of site location on the Earth’s surface. For graphical clarity, the plots use different scales. For both flow types, the main features are not only themselves similar, but also strikingly share many features of the magnetic intensity $F$ itself (Fig. 1): pairs of high-latitude intensity patches, a strong low in the Pacific and moderately low equatorial values. This link is simply explained by the fact that high rates of change occur preferentially in regions where the field intensity is high, since the secular variation is linear and homogeneous in the magnetic field itself.

Lastly, future archeomagnetic studies may seek to find archeomagnetic spikes contemporaneous with that at Timna-30. We now ask what form the secular variation would take on the Earth’s surface in the case of a dual-spike, by considering two possible (illustrative) dual-sites in (i) Xian, China and (ii) San Lorenzo, Mexico (i.e. in Mesoamerica), where archeology is particularly rich. Mathematically, we seek the maximum value of $dF/dt$ that occurs (equally) at two sites, $S$ and $S_2$, subject to the same rms constraint on the flow as before, which becomes

$$ \max_{\mathbf{q}}[\mathbf{G}^T\mathbf{q} - \mu[\mathbf{G}^T\mathbf{q} - \mathbf{G}_2^T\mathbf{q}]] - \lambda(\mathbf{q}^T\mathbf{E}\mathbf{q} - T_0^2)], \quad (16) $$

where $\mathbf{G}_2$ is the analogue of $\mathbf{G}$ at site $S_2$, and $\lambda$ and $\mu$ are Lagrange multipliers enforcing the constraints. As before, the (unnormalised) optimum flow is

$$ \mathbf{q} = \frac{1}{2\lambda} E^{-1}[(1 - \mu)\mathbf{G} + \mu\mathbf{G}_2] \quad (17) $$
<figure><img src="image_16.png" /><figcaption>(b) Purely-toroidal flows</figcaption></figure>

<figure><img src="image_17.png" /></figure>

**Fig. 8.** Contours of $dF/dt$ (in nT/yr) on the Earth's surface (left), a regional view (middle) and contours of $|u|$ with indicative arrows showing the structure of the optimal flow (right) for the realisation of ST_CHAOS giving the median value of upper bounds for (a) unrestricted (top row) and (b) purely-toroidal (bottom row) flows, optimised for a site at Timna-30 (shown by the stars).

<figure><img src="image_18.png" /></figure>

**Fig. 9.** Contour plot of the mean (over 2000 realisations) of the upper bound of $dF/dt$ as a function of position of the observation site on the Earth's surface, for the magnetic field model ST_CHAOS for (a) unrestricted and (b) purely-toroidal flows. Note that the plots have different scales.
(i)(a) Timna-30 and Xian;
Unrestricted flow

<figure><img src="image_19.png" /></figure>

(ii)(a) Timna-30 and San Lorenzo;
Unrestricted flow

<figure><img src="image_20.png" /></figure>

(i)(b) Timna-30 and Xian;
Purely-toroidal flow

<figure><img src="image_21.png" /></figure>

(ii)(b) Timna-30 and San Lorenzo;
Purely-toroidal flow

<figure><img src="image_22.png" /><figcaption>**Fig. 10.** Contour plot of $dF/dt$ over the Earth's surface, for the dual-site calculations: (i) Timna-30 and Xian, China; (ii) Timna-30 and San Lorenzo, Mexico. The magnetic field is a single realisation of model ST_CHAOS. The flows (unrestricted and purely-toroidal) are chosen to give a maximal and equal value of $dF/dt$ at both sites.</figcaption></figure>

where the required value of the weighting, $\mu$, ensuring equal values of $dF/dt$ at both sites, is found using a simple bisection algorithm.

Fig. 10 shows contours of $dF/dt$ for both dual-site problems with unrestricted and purely-toroidal flows, for a single (arbitrarily chosen) realisation of the magnetic field model ST_CHAOS. In case (i), where the two sites are relatively close, there is a single peak of $dF/dt$ which is in between the sites. By contrast, in case (ii) where the two sites are of sufficient distance apart, the flows and corresponding intensity changes are quite separate. Thus, in the event of case (i), the same archeointensity spike would necessarily be realised over much of Asia but in (ii) the spikes may not be realised except in the vicinity of the two sites. Note that this analysis does not preclude potential localised spikes in (i) for flow structures different to the optimal cases derived here. In each case, the maximum values of $dF/dt$ are about 70% of those for the single-site case.

Lastly, we commented earlier that, just as the present-day magnetic field, the one-site case had a secular variation dominated by $l=2$. Here, of the four cases presented, only case (i)(a) has a secular variation on the Earth's surface dominated by degree 2, all other cases have a dipole-dominated spectrum.

**6. Concluding remarks on the results**

In this study we have described a methodology for computing the maximum value of intensity change, $dF/dt$, at any site on the Earth's surface subject to knowledge of the magnetic field (which we handled using statistical methods) and a rms constraint on the core-surface flow of 13 km/yr, a consensual figure from the core-flow inversion literature. For the site of Timna-30,

of the two magnetic field models used, that based on CHAOS-4 had a higher mean upper bound (~1.2 µT/yr) than that based on CALS10k.1b (0.95 µT/yr). Restricting attention to purely-toroidal flows, the mean upper bounds are lowered by about a factor of 2 to 0.62 µT/yr and 0.46 µT/yr. None of these bounds are consistent with (i.e. greater than) the values of $dF/dt$ required by the archeomagnetic spikes of ~4–5 µT/yr suggested for the Near East (recall Fig. 2). A further point of note is that the largest, but poorly constrained, variation rates of up to 0.8 µT/yr in European regions during the Medieval period suggested by Gómez-Paccard et al. (2012) are either barely consistent (in the unrestricted case) or inconsistent (in the toroidal case) with our upper bounds.

The extent of the discrepancy between the upper bounds and the severity of the archeomagnetic spikes depends on which flow assumption we take to be the most physically meaningful. The two choices of flow structure we studied were deliberately chosen as end-members of a spectrum of permissible poloidal flow structures in the core: at one end they are unrestricted, at the other constrained to be zero. It is well-known that a purely-toroidal flow can adequately explain the observed secular variation within the observational uncertainties. But due to the inherent nonuniqueness of the core-flow problem, many flows with a non-zero poloidal component provide as good a fit (see e.g. Holme, 2007). However, a large body of work in the core-flow literature points to a predominantly toroidal large-scale flow structure, even when the toroidal constraint is relaxed during the inversion. Indeed, in many such studies the ratio of the toroidal to poloidal flow energy is reported to be at least ten (Eymin and Hulot, 2005; Holme and Olsen, 2006; Pais and Jault, 2008; Gillet et al., 2009; Lesur et al., 2010). Furthermore, very recent work on the physical properties of the core itself suggests that the outermost part
of the liquid core is stratified (see Hirose et al., 2013, for a re-
cent review), thereby admitting only toroidal flows. This motivates
us to adopt the upper bounds derived from purely-toroidal flows,
as the most physically compelling. The archeomagnetic spikes then
require a rate of change of intensity of around 6 to 8 times larger
than what our most favourable toroidal flows can account for.

A generous threefold increase (up to a flow rms of 40 km/yr) allows the toroidal bound to reach 1.5 µT/yr, but this still falls short of 4–5 µT/yr. Furthermore, and from a geometrical standpoint, this issue is compounded by the fact the optimising flows are very much more spatially focused than the global-flows inferred for the present-era. It would be rather fortuitous indeed if the only strong flow were to occur locally beneath Timna-30 at 1000 BC. Indeed, it is likely that such localised flows would never arise in the core, since due to the important role of the Coriolis force (and associated quasi-geostrophic flow), any rapid dynamics would be expected to be equatorially symmetric (Jault, 2008; Gillet et al., 2011), a property that our specialised optimal flows do not have.

Regarding the possibility of contemporaneous spikes of the
same severity at Timna-30 and at another site (e.g. Xian or San
Lorenzo), the upper bound is about 70% of that only for Timna-30.
This reduction in the bound is simply a reflection of the fact that
the kinetic energy available is spread more thinly over the CMB,
and less is then available to drive rapid change at any one site.
Thus simultaneous dual spikes of ~5 µT/yr are even more unlikely
than a spatially localised spike in the Near East. Furthermore, con-
firmation of the Near-East spikes could likely only be attempted in
the immediate vicinity (see Fig. 8). It is however interesting to note
that a spike in Xian, China, around 1000 BC would require either
an even more intense contemporaneous spike in an intermediate
region between the Near East and China, or for the two events to
be completely independent and asynchronous in time (for which
the more generous single-site bounds, rather than the lower dual-
site bounds, would apply). The site in San Lorenzo is sufficiently
far away from Timna-30 that two simultaneous spikes could occur
without any signal in intermediate regions.

7. Discussion

We finish by briefly discussing how the discrepancy in our com-
putations and the archeomagnetic intensity spikes might be recti-
fied.

A simple resolution is to reduce the rate of change of archeo-
magnetic intensity found in the data, either by reassessing (and
then reducing) the intensity itself or by making the duration of
the intensity spike longer. Using our preferred upper bound of
0.62 µT/yr, this would necessitate a spike lasting several centuries
rather than a few decades. The gap in the archeointensity data ob-
tained from Syria admits such a possibility of an extreme intensity
peak between 1100 BC and 800 BC. Moreover a longer duration
would be more consistent with the fact that a spike-like feature is
also observed at Arslantepe from a kiln with age uncertainties of
300 yr. However, increasing the duration of the spike (which would
then become a "simple" archeointensity maximum) would severely
contradict the stratigraphy and age models proposed by Ben-Yosef
et al. (2009) and Shaar et al. (2011) in their archeological sites. It is
further surprising that spike-like events have not been recorded in
the Greek and Bulgarian data sets, although dF/dt there should be
~0.55 µT/yr according to the results from ST_CHAOS (Fig. 8(b)).
However, it is possible that this is due solely to the poor documen-
tation of such an era, rather than the non-occurrence of a spike:
this time interval in Eastern Europe and Near East corresponds to
the ancient Dark Age (~1100–800 BC) after the sudden demise of
the flourishing Late Bronze societies, and thus a scarcity of arche-
ological artifacts.

Alternatively, since a faster flow leads to more rapid intensity
change, our preferred purely-toroidal bounds could become consis-
tent with the archeomagnetic intensity changes if we increased our
rms flow by a factor of 8 to ~100 km/yr. However, such a flow is
far away from what we consider to be the current dynamics of the
core, requiring significant changes in the core-flow since 1000 BC.
Furthermore, such a strong flow is not consistent with the under-
standing gained over the past few years from scaling laws: the
rms flow speed inside the core is controlled by the available con-
ective power (e.g. Christensen, 2010), which itself varies on the
time-scales of the secular cooling of the Earth, and can be consid-
ered constant throughout the Holocene. Although fluctuations to
this constant time-averaged flow are possible, deviations of up to
8 times the mean would be required here, which we deem not
possible.

As a further alternative, it is conceivable that the magnetic field
in 1000 BC was similar to none of our Monte Carlo realisations,
and that our computed upper bounds significantly underestimate
dF/dt for this epoch. One possibility is that the geomagnetic field
models based on CALS10k1.b and CHAOS-4 are simply too weak:
the values of dF/dt (linear and homogeneous in the magnitude
of B) therefore being too small. However, a magnetic field stronger
by a factor of 8 (as required in the purely-toroidal case) is not con-
sistent with records of the axial dipole moment over the Holocene
(Knudsen et al., 2008). Another possibility is a strong small-scale
magnetic field, but this surely would have been realised, at the
very least, in the ST_CALS stochastic model. Although it is true
that ST_CALS showed a higher variance than ST_CHAOS, 5 µT/yr is
around 100 standard deviations above the most optimistic bound
of 0.95 ± 0.046 µT/yr (with unrestricted flows). The most likely
shortcoming of our CHAOS-4 models is the fact that the axial
dipole field moment was stronger at the beginning of the first mil-
lennium BC than in more recent times (e.g. Knudsen et al., 2008).
For this reason, we performed extra calculations assuming that the
axial dipole coefficient in CHAOS-4 was twice that of its nomi-
nal value. This change led to an increase of the upper bounds
by approximately 35% and 10% for the unrestricted and toroidal
core-flow patterns, respectively. These rather moderate variations,
in particular in the case of our preferred toroidal flows, are clearly
not sufficient to reconcile the spike data with our computations.

Lastly, we concede that an explanation for the spikes may lie
outside the scope of this study. With this in mind, we briefly ex-
amine the only physics of magnetic induction that our treatment
neglects: that of radial diffusion of field. Although it is tempt-
ing to argue that diffusion in the radial and horizontal directions
should have a comparable effect, this may be too simplistic. Mag-
netic features with short radial length scale would diffuse rapidly
through the CMB and therefore, plausibly, could be important. The
dynamic origin of these localised features is far from obvious how-
ever. One possibility are upwellings that carry strong, deep field
to the surface, the concentrated flux then being expelled through
the CMB by diffusion. Such a mechanism was for instance consid-
ered by Chulliat et al. (2010) to explain the recent acceleration of
the north magnetic pole drift during the 1990s. However, a flux-
expulsion event (Bloxham, 1986) not only requires a non-stratified
outer core, but it remains to be seen whether it can deliver the re-
quired change in intensity over ~one decade, very short compared
to a core-overturn time-scale of ~one century. Modelling of these
dynamics and whether or not such a mechanism can explain the
spikes is a future extension of our analysis. Looking more broadly,
further data from the Near East would of course be highly benefi-
cial in order to add to (and perhaps to refute or corroborate) the
existing studies of Shaar et al. (2011) and Ben-Yosef et al. (2009).
Moreover, it would be very useful to extend the analysis to di-
rectional data: the datasets that motivated this study contain only
records of archeomagnetic intensity and not directional data, since
the slag materials are, in most cases, displaced from their cooling position. Should the observation of extreme rapid intensity change stand up to scrutiny, its implications would be far reaching, and may require consideration of yet poorly known core processes.

## Acknowledgements

We would like to acknowledge useful discussions with Steve Hunter and Chris Davies; comments by the deep-Earth research group at Leeds helped improve the manuscript. This study benefited from the support of the International Space Science Institute (ISSI), through the work of team 241. PWL was supported by NERC grant NE/G014043/1. The contribution of AF and YG is IPGP contribution 3455. The numerical calculations were performed on the super-computer ARC1 at the University of Leeds.

## Appendix A. Supplementary material

Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.epsl.2013.11.020.

## References

Backus, G., Parker, R., Constable, C., 1996. Foundations of Geomagnetism. CUP.

Ben-Yosef, E., Ron, H., Tauxe, L., Agnon, A., Genevey, A., Levy, T.E., Avner, U., Najjar, M., 2008. Application of copper slag in geomagnetic archaeointensity research. J. Geophys. Res. 113.

Ben-Yosef, E., Tauxe, L., Levy, T.E., Shaar, R., Ron, H., Najjar, M., 2009. Geomagnetic intensity spike recorded in high resolution slag deposit in Southern Jordan. Earth Planet. Sci. Lett. 287, 529–539.

Bloxham, J., 1986. The expulsion of magnetic flux from the Earth's core. Geophys. J. R. Astron. Soc. 87, 669–678.

Christensen, U.R., 2010. Dynamo scaling laws and applications to the planets. Space Sci. Rev. 152, 565–590.

Christensen, U.R., Wicht, J., 2007. Numerical dynamo simulations. In: Olson, P., Schubert, G. (Eds.), Core Dynamics. In: Treatise on Geophysics, vol. 8. Elsevier, Oxford, pp. 245–282. Chapter 8.

Christensen, U.R., Aubert, J., Hulot, G., 2010. Conditions for Earth-like geodynamo models. Earth Planet. Sci. Lett. 296, 487–496.

Chulliat, A., Hulot, G., Newitt, L.R., 2010. Magnetic flux expulsion from the core as a possible cause of the unusually large acceleration of the north magnetic pole during the 1990s. J. Geophys. Res. 115.

De Marco, E., Spatharas, V., Gomez-Paccard, M., Chauvin, A., Kondopoulou, D., 2008. New archaeointensity results from archaeological sites and variation of the geomagnetic field intensity for the last 7 millennia in Greece. Phys. Chem. Earth 33, 578–595.

Donadini, F., Korte, M., Constable, C.G., 2009. Geomagnetic field for 0–3 ka: 1. New data sets for global modeling. Geochem. Geophys. Geosyst. 10, Q06007.

Ertepinar, P., Langereis, C.G., Biggin, A.J., Frangipane, M., Matney, T., Ökse, T., Engin, A., 2012. Archaeomagnetic study of five mounds from upper Mesopotamia between 2500 and 700 BCE: Further evidence for an extremely strong geomagnetic field ca. 3000 years ago. Earth Planet. Sci. Lett. 357, 84–98.

Eymin, C., Hulot, G., 2005. On core surface flows inferred from satellite magnetic data. Phys. Earth Planet. Inter. 152, 200–220.

Finlay, C., Amit, H., 2011. On flow magnitude and field-flow alignment at Earth's core surface. Geophys. J. Int. 186, 175–192.

Finlay, C.C., Dumberry, M., Chulliat, A., Pais, M.A., 2010. Short timescale core dynamics: Theory and observations. Space Sci. Rev. 155, 177–218.

Gallet, Y., Al-Maqdissi, M., 2010. Archéomagnétisme à Mishirfeh-Qatna: Nouvelles données sur l'évolution de l'intensité du champ magnétique terrestre au Moyen-Orient durant les derniers millénaires. Akkadica 131, 29–46.

Gallet, Y., Le Goff, M., 2006. High-temperature archeointensity measurements from Mesopotamia. Earth Planet. Sci. Lett. 241, 159–173.

Gallet, Y., Genevey, A., Le Goff, M., Fluteau, F., Ali Eshraghi, S., 2006. Possible impact of the Earth's magnetic field on the history of ancient civilizations. Earth Planet. Sci. Lett. 246, 17–26.

Gallet, Y., Le Goff, M., Genevey, A., Margueron, J., Matthiae, P., 2008. Geomagnetic field intensity behavior in the Middle East between ~3000 BC and ~1500 BC. Geophys. Res. Lett. 35, 2307.

Gallet, Y., D'Andrea, M., Genevey, A., Pinnock, F., Le Goff, M., Matthiae, P., in press. Archeomagnetism at Ebla (Tell Mardikh, Syria). New data on geomagnetic field intensity variations in the Near East during the Bronze Age. J. Archaeol. Sci. http://dx.doi.org/10.1016/j.jas.2013.11.007.

Genevey, A., Gallet, Y., Margueron, J.C., 2003. Eight thousand years of geomagnetic field intensity variations in the Eastern Mediterranean. J. Geophys. Res. 108, 2228.

Genevey, A., Gallet, Y., Constable, C.G., Korte, M., Hulot, G., 2008. Archeoint: An upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment. Geochem. Geophys. Geosyst. 9, Q04038.

Genevey, A., Gallet, Y., Rosen, J., Le Goff, M., 2009. Evidence for rapid geomagnetic field intensity variations in western Europe over the past 800 years from new French archeointensity data. Earth Planet. Sci. Lett. 284, 132–143.

Genevey, A., Gallet, Y., Thébault, E., Jesset, S., Le Goff, M., 2013. Geomagnetic field intensity variations in western Europe over the past 1100 years. Geochem. Geophys. Geosyst. 14 (8). http://dx.doi.org/10.1002/ggge.20165.

Gillet, N., Pais, A., Jault, D., 2009. Ensemble inversion of time-dependent core flow models. Geochem. Geophys. Geosyst. 10.

Gillet, N., Schaeffer, N., Jault, D., 2011. Rationale and geophysical evidence for quasi-geostrophic rapid dynamics within the Earth's outer core. Phys. Earth Planet. Inter. 187, 380–390.

Gómez-Paccard, M., Chauvin, A., Lanos, P., Dufresne, P., Kovacheva, M., Hill, M.J., Bea-mud, E., Blain, S., Bouvier, A., Guibert, P., Archaeological Working Team, 2012. Improving our knowledge of rapid geomagnetic field intensity changes observed in Europe between 200 and 1400 AD. Earth Planet. Sci. Lett. 355, 131–143.

Hervé, G., Chauvin, A., Lanos, P., 2013. Geomagnetic field variations in western Europe from 1500 BC to 2000 AD. Part II: New intensity secular variation curve. Phys. Earth Planet. Inter. 218, 51–65.

Hirose, K., Labrosse, S., Hernlund, J., 2013. Composition and state of the core. Annu. Rev. Earth Planet. Sci. 41, 657–691.

Holme, R., 2007. Large-scale flow in the core. In: Olson, P. (Ed.), In: Treatise on Geophysics, vol. 8. Elsevier, pp. 107–130.

Holme, R., Olsen, N., 2006. Core surface flow modelling from high-resolution secular variation. Geophys. J. Int. 166, 518–528.

Hulot, G., Sabaka, T.J., Olsen, N., 2007. The present field. In: Kono, M., Schubert, G. (Eds.), Geomagnetism. In: Treatise on Geophysics, vol. 5. Elsevier, pp. 33–75. Chapter 2.

Jackson, A., Finlay, C.C., 2007. Geomagnetic secular variation and its application to the core. In: Kono, M., Schubert, G. (Eds.), Geomagnetism. In: Treatise on Geophysics, vol. 5. Elsevier, Amsterdam, pp. 148–193. Chapter 5.

Jackson, A., Livermore, P., 2009. On Ohmic heating in the Earth's core I: nutation constraints. Geophys. J. Int. 179, 923–928.

Jault, D., 2008. Axial invariance of rapidly varying diffusionless motions in the Earth's core interior. Phys. Earth Planet. Inter. 166, 67–76.

Knudsen, M.F., Riisager, P., Donadini, F., Snowball, I., Muscheler, R., Korhonen, K., Pesonen, L.J., 2008. Variations in the geomagnetic dipole moment during the Holocene and the past 50 kyr. Earth Planet. Sci. Lett. 272, 319–329.

Korte, M., Constable, C., Donadini, F., Holme, R., 2011. Reconstructing the Holocene geomagnetic field. Earth Planet. Sci. Lett. 312, 497–505.

Kovacheva, M., Boyadziev, Y., Kostadinova-Avramova, M., Jordanova, N., Donadini, F., 2009. Updated archeomagnetic data set of the past 8 millennia from the Sofia laboratory, Bulgaria. Geochem. Geophys. Geosyst. 10, Q05002.

Lesur, V., Wardinski, I., Asari, S., Minchev, B., Mandea, M., 2010. Modelling the Earth's core magnetic field under flow constraints. Earth Planets Space 62, 503–516.

Licht, A., Hulot, G., Gallet, Y., Thébault, E., 2013. Ensembles of low degree archeomagnetic field models for the past three millennia. Phys. Earth Planet. Inter. 224, 38–67.

McLeod, M.G., 1996. Spatial and temporal power spectra of the geomagnetic field. J. Geophys. Res. 101, 2745–2763.

Olsen, N., Luehr, H., Sabaka, T.J., Michaelis, I., Rauberg, J., Toffner-Clausen, L.L., 2010. CHAOS-4 – A high-resolution geomagnetic field model derived from low-altitude CHAMP data. AGU abstract GP21A-0992 http://www.spacecenter.dk/files/magnetic-models/CHAOS-4/.

Pais, M.A., Jault, D., 2008. Quasi-geostrophic flows responsible for the secular variation of the Earth's magnetic field. Geophys. J. Int. 173, 421–443.

Pozzo, M., Davies, C., Gubbins, D., Alfè, D., 2013. Transport properties for liquid silicon-oxygen-iron mixtures at Earth's core conditions. Phys. Rev. B 87, 014110.

Roberts, P.H., Scott, S., 1965. On analysis of the secular variation. 1. A hydromagnetic constraint: Theory. J. Geomagn. Geoelectr. 17, 137–151.

Shear, R., Ben-Yosef, E., Ron, H., Tauxe, L., Agnon, A., Kessel, R., 2011. Geomagnetic field intensity: How high can it get? How fast can it change? Constraints from Iron Age copper slag. Earth Planet. Sci. Lett. 301, 297–306.

Tema, E., Kondopoulou, D., 2011. Secular variation of the Earth's magnetic field in the Balkan region during the last eight millennia based on archaeomagnetic data. Geophys. J. Int. 186, 603–614.

Thébault, E., Gallet, Y., 2010. A bootstrap algorithm for deriving the archeomagnetic field intensity variation curve in the middle east over the past 4 millennia BC. Geophys. Res. Lett. 37.

Yang, S., Odah, H., Shaw, J., 2000. Variations in the geomagnetic dipole moment over the last 12 000 years. Geophys. J. Int. 140, 158–162.