# Continuous millennial decrease of the Earth's magnetic axial dipole

Wilbor Poletti, Andrew J. Biggin, Ricardo I. F. Trindade, Gelvam A. Hartmann, Filipe Terra-Nova

► To cite this version:

Wilbor Poletti, Andrew J. Biggin, Ricardo I. F. Trindade, Gelvam A. Hartmann, Filipe Terra-Nova. Continuous millennial decrease of the Earth's magnetic axial dipole. Physics of the Earth and Planetary Interiors, 2018, 274, pp.72-86. 10.1016/j.pepi.2017.11.005. insu-03663869

HAL Id: insu-03663869

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

Submitted on 10 May 2022

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
# Continuous millennial decrease of the Earth's magnetic axial dipole

Wilbor Poletti<sup>a,b,*</sup>, Andrew J. Biggin<sup>b</sup>, Ricardo I.F. Trindade<sup>a</sup>, Gelvam A. Hartmann<sup>c</sup>,
Filipe Terra-Nova<sup>d</sup>

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

## ARTICLE INFO

**Keywords:**

Geomagnetism

Archaeointensity

Axial dipole

Archaeomagnetic dataset

Data selection criteria

## ABSTRACT

Since the establishment of direct estimations of the Earth's magnetic field intensity in the first half of the nineteenth century, a continuous decay of the axial dipole component has been observed and variously speculated to be linked to an imminent reversal of the geomagnetic field. Furthermore, indirect estimations from anthropologically made materials and volcanic derivatives suggest that this decrease began significantly earlier than direct measurements have been available. Here, we carefully reassess the available archaeointensity dataset for the last two millennia, and show a good correspondence between direct (observatory/satellite) and indirect (archaeomagnetic) estimates of the axial dipole moment creating, in effect, a proxy to expand our analysis back in time. Our results suggest a continuous linear decay as the most parsimonious long-term description of the axial dipole variation for the last millennium. We thus suggest that a break in the symmetry of axial dipole moment advective sources occurred approximately 1100 years earlier than previously described. In addition, based on the observed dipole secular variation timescale, we speculate that the weakening of the axial dipole may end soon.

## 1. Introduction

The continuous intensity record of the Earth's magnetic field was started in 1833 CE by Carl Friedrich Gauss, enabling the precise direct recording of the full geomagnetic vector for the past 184 years (e.g. Kono (2007)). Nonetheless, ancient civilizations, when baking pottery, were also inadvertently recording the Earth's magnetic field. This archaeomagnetic record can be retrieved from ancient baked clay (and from historical lavas) using laboratory techniques developed more than one hundred years ago (Folgheraiter, 1899), that were subsequently significantly improved (Thellier and Thellier, 1959; Coe, 1967; Coe et al., 1978; Aitken et al., 1988; Shaw et al., 1996; Riisager and Riisager, 2001; Yu et al., 2004; Le Goff and Gallet, 2004). Archaeomagnetism provides information about geomagnetic field variations thousands of years before the “Gauss era” and can help in unveiling the processes operating in the Earth's core at time-scales longer than the past 184 years (e.g. Dumberry and Finlay (2007), Amit et al. (2011), Sanchez et al. (2016), Terra-Nova et al. (2015, 2016)).

Variations observed in intensity data from observatories, satellites,

volcanic lavas, and archaeological artifacts can be linked to the main component of the geomagnetic field (e.g., Jackson et al. (2000), Olson and Amit (2006), Gubbins et al. (2006), Finlay (2008), Korte et al. (2009), Korte and Constable (2011), Suttie et al. (2011), Licht et al. (2013), Nilsson et al. (2014), Pavón-Carrasco et al. (2014)), which originates from the movement of the outer core's conductive fluid and is dominated by the axial dipole component. Thanks to the continuous direct records over a wide spatial coverage during the Gauss era, it was possible to describe the geomagnetic dipole variation in detail for the past 184 years (e.g. Jackson et al. (2000), Gillet et al. (2013), Finlay et al. (2015)). For this period, the decay rate of the axial dipole is, on average, about 15 nT/yr, with decadal fluctuations (Jackson et al., 2000; Finlay et al., 2015). Prior to the Gauss era, the Earth's magnetic field record provided by archaeomagnetism is still scarce both temporally and spatially (Genevey et al., 2008, Donadini et al., 2009, Brown et al., 2015; Poletti et al., 2016). Yet, it is sufficiently robust for the description of local, rapid variations (de Groot et al., 2013; Genevey et al., 2016). It is also the only means to analyze geomagnetic axial dipole evolution on millennial timescales.

<footnote><sup>a</sup> Universidade de São Paulo, Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Departamento de Geofísica, Laboratório de Paleomagnetismo – USPMag, Rua do Matão, 1226, 05508-090 São Paulo, Brazil</footnote>

<footnote><sup>b</sup> Geomagnetism Laboratory, Department of Earth, Ocean and Ecological Sciences, School of Environmental Sciences, University of Liverpool, Liverpool, UK</footnote>

<footnote><sup>c</sup> Instituto de Geociências, Universidade Estadual de Campinas, Rua João Pandiá Calógeras, 51, 13083-870 Campinas, Brazil</footnote>

<footnote><sup>d</sup> CNRS, Université de Nantes, Nantes Atlantiques Universités, UMR CNRS 6112, Laboratoire de Planétologie et de Géodynamique, 2 rue de la Houssinière, F-44000 Nantes, France</footnote>

<footnote><sup>*</sup> Corresponding author at: Universidade de São Paulo, Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Departamento de Geofísica, Laboratório de Paleomagnetismo – USPMag, Rua do Matão, 1226, 05508-090 São Paulo, Brazil.</footnote>

<footnote>E-mail address: wilbor-poletti@iag.usp.br (W. Poletti).</footnote>
Several datasets of full vector archaeomagnetic data exist (e.g. Brown et al. (2015), Arneitz et al. (2017b)). From such datasets, different descriptions of global variations of the geomagnetic axial dipole have emerged (e.g. Valet et al. (2008), Genevey et al. (2008), Knudsen et al. (2008), Usoskin et al. (2016)). However, although such efforts have produced useful and detailed descriptions of millennial timescale variations, they have tended to avoid making links between surface observations and Earth's core process. Furthermore, due to fundamental differences between the magnetic intensity records obtained by direct and indirect measurements (e.g., spatial and temporal coverage, experimental errors) (Arneitz et al., 2017a), the geomagnetic axial dipole variations are usually described independently for two distinct periods: before and after 1840 CE and there have been few attempts to critically compare and integrate them. From 1840 to today, the axial dipole variations are robust, meanwhile for the period of 1590–1840 CE, the widely utilized historical field model gufm1 (Jackson et al., 2000) uses an arbitrary extrapolation of the axial dipole intensity from the Gauss era, whereas estimations incorporating only archaeointensity data tend to favor a rather flat decay of the axial dipole field (Gubbins et al., 2006; Finlay, 2008; but also see Suttie et al. (2011)).

In this work, we present a new description of geomagnetic axial dipole variations before the Gauss era by evaluating the Axial Dipole Moment (ADM) and Virtual Axial Dipole Moment (VADM) obtained from archaeointensity data for the entirety of the last two millennia. To do so, we accepted only high-quality archaeointensity estimates into our evaluation, and attempted to assess these data using a minimum number of linear trends. Our compilation indicates a significant shift in the trend of the axial dipole strength around the interval 550–750 CE, initiating a continuous decay in the same order of magnitude of the Gauss era up to the present. We attribute this shift to fundamental changes in geodynamo workings in the last millennium, ultimately attempting to link the archaeointensity record to dynamical processes within Earth's outer core.

## 2. Methods

### 2.1. Datasets

Thellier and Thellier (1959) defined the original double-heating protocol (TT) which today incorporates checks for alteration (Coe et al., 1978) and multi-domain effects (e.g. Riisager and Riisager (2001)), as well as corrections for the effects of magnetic anisotropy (e.g. Rogers et al. (1979)) and for the fast cooling-rates applied in the laboratory (e.g. Fox and Aitken (1980)). Other methods such as the Microwave (MW; Shaw et al., 1996) and Triaxe (TR; Le Goff and Gallet, 2004) have also been developed and their results have been systematically compared one to each other, thus increasing our confidence in the results (e.g., TT-TR: Le Goff and Gallet (2004), Gallet and Le Goff (2006), Genevey et al. (2009), Hartmann et al. (2010, 2011); TT-MW: Shaw et al. (1999), Hill et al. (2002a,b), Casas et al. (2005), Stark et al. (2010), Ertepinar et al. (2016); TT-TR-MW: Poletti et al. (2013)). A detailed historical and physical description of the Thellier-Thellier method and its modifications was put forward by Dunlop (2011).

GEOMAGIA50.v3 (Brown et al., 2015) is a comprehensive database comprising 14,645 data (declination, inclination, and intensity) from archaeological artifacts and volcanic material, obtained over the past half century. In our analysis we used a catalogue of archaeointensities from the GEOMAGIA50.v3 database, and some other data recently published that were not incorporated into the collection at the time of our analysis (Table A1). The time window investigated was the past two millennia as this period shows the best temporal and spatial coverage.

### 2.2. Selection criteria

Data selection was performed by checking if current laboratory criteria were satisfied (e.g. Poletti et al. (2016)). We considered seven

factors in our assessment of the archaeointensity data when considering the archaeological material. The factors were applied following the sequence in which they are presented below:

(i) Age uncertainty. For this study, we accepted data with age uncertainty less than or equal to 100 years ($\text{O}_{\text{age}} \le 100$). This rather strict choice was made to enable the comparison between archaeomagnetic and observatory/satellite data in the Gauss era (i.e., 184 years). Data were not filtered by the dating technique (except that archaeomagnetic dates were excluded);

(ii) The archaeointensity method used and the protocol adopted. We only accepted intensity data performed exclusively with the classical double-heating method at room-temperature (Thellier and Thellier, 1959) in one of its modified versions (TT) (Coe, 1967; Aitken et al., 1988; Yu et al., 2004), the microwave method (MW) (Shaw et al., 1996; Hill and Shaw, 1999), or the high-temperature Triaxe method (TR) (Le Goff and Gallet, 2004). Our choice was based on palaeointensity methods that perform a gradual and progressive replacement between the magnetizations acquired from the nature and laboratory. The results obtained from these three specific methods are more likely to be high-quality and concordant as highlighted by several works published in the last few decades (e.g. Hill et al. (2002a), Genevey et al. (2009), Poletti et al. (2013));

(iii) Additional steps to check alterations during the experiment. For TT and MW, we required additional steps in the laboratory protocol, referred to as pTRM checks, to monitor possible (thermo)chemical alterations during the gradual increase of temperature (TT) or power (MW) steps on the experiment (Coe et al., 1978). For TR, these additional steps are unnecessary (Le Goff and Gallet, 2004);

(iv) Evaluation of the influence of multi-domain (MD) grains. We required at least one test-type to verify possible MD grains influence (e.g. Riisager and Riisager (2001), Krása et al. (2003), Yu et al. (2004)), in order to avoid the violation of the principles of additivity and reciprocity, which are part of the backbone of the Thellier-Thellier method (Yu and Dunlop, 2003; Dunlop, 2011);

(v) Anisotropy thermoremanent magnetization (ATRM) correction. We accepted only data largely unbiased by anisotropy effects either by having the laboratory field applied in a direction within 10 degrees of the principal component of the natural remanent magnetization (NRM) (Rogers et al., 1979; Aitken et al., 1981), or by the correction of the tensor of ATRM being obtained experimentally and calculated through the formulation proposed by Veitch et al. (1984). Although there are other ways to correct the ATRM effect, for example, through the tensor obtained from measures of anhysteretic remanent magnetization (ARM) or magnetic susceptibility (MS), we restrict our analysis to results that take into account the same physical basis between anisotropy correction and Thellier-Thellier method (see ii). Data corrected by the ATRM effect using ARM or MS technique implicitly assume equivalence between the pairs of anisotropy tensors TRM-ARM or TRM-MS, which are not always true (Stephenson et al., 1986; Yu et al., 2003), although we acknowledge the need for further advances in this topic.

(vi) Cooling rate correction. We accepted only archaeointensity data that were corrected for cooling rate effects following the experimental procedure described by Chauvin et al. (2000) and Genevey and Gallet (2002) for data from TT, and Poletti et al. (2013) for data from MW, in order to avoid possible bias in the final archaeointensity result due to the difference between natural (NRM) and experimental (pTRMs imparted) cooling times (e.g. Fox and Aitken (1980); Dodson and McClelland-Brown (1980), Halgedahl et al. (1980), Biggin et al. (2013)). All results from TR were accepted without this correction, since TR routinely produces results consistent with cooling rate-corrected TT and MW estimates (e.g. Genevey et al. (2009), Hartmann et al. (2010), Hartmann et al.
(2011), Poletti et al. (2013));

(vii) Standard deviation of final archaeointensity estimates. We only accepted data with standard deviation up to 15% of the mean intensity (Paterson et al., 2014), and a minimum of three samples/specimens (N ≥ 3) per age.

In archaeointensity works, there are the designations “sample/specimen”, “fragment” and “site”, which are employed to distinguish the data acquired from laboratory measurements from the raw materials utilized, as well as for the calculation/presentation of the final archaeointensity results. Although there is a consensus about the term site, which is the final value of the magnetic intensity from a specific location for a given age (e.g., magnetic intensity for an archaeological ruin or destruction level that represents a specific period, magnetic intensity for a specific lava flow, etc.), there is no uniformity of nomenclature in relation to the other terms, which in turn are used in the calculations of the means. Thus, there are several ways in which the calculated mean from measured data in the laboratory is associated with a site (e.g., mean of several samples/specimens from a single fragment, mean of several fragments with a single sample each, mean of several fragments with several samples each, etc). We understand that the result of a site can be given by the mean value obtained from the results of at least three independent fragments; that the result of each fragment is given by the mean value from the results of at least two independent samples/specimens extracted from the same fragment in question; and that the result of each sample/specimen is given by the value obtained in the laboratory, processed, analyzed, approved by current selection criteria (e.g. Paterson et al. (2014)), and corrected for the possible anisotropy and cooling rate effects (e.g. Genevey et al. (2009), Hartmann et al. (2010, 2011), Poletti et al. (2016)). Due to the non-uniformity about the nomenclatures described, we did not distinguish between the results presented as site in our assessment. However, it is important to emphasize that our final results already have strong restrictions due to the previously applied selection criteria (i–vi). Finally, we suggest that future works should include more details about the distinction made between site, fragment and sample/specimen.

For volcanic rocks we applied the same selection criteria as for data obtained from archaeological materials with some exceptions. In criterion ii) the pseudo-Thellier method proposed for volcanic rocks (de Groot et al., 2013) (only 18 entries) were also accepted; in criterion v) an ATRM correction was not required; in criterion vi) the cooling rate correction was neglected once its effect has been shown to be very small for the assemblages of PSD and interacting SD grains that are most commonly found in lavas (Biggin et al., 2013).

## 2.3. ADMs and VADMs

From filtered archaeointensity data ($B_{\text{indirect}}$), axial dipole moments ($\text{ADMs}_{\text{indirect}}$) are calculated using the theorem of Hulot et al. (1997), following the strategy applied by previous works (Gubbins et al., 2006; Genevey et al., 2009; Hartmann et al., 2011), described as:

$$g_1^0(t)_{\text{indirect}} = g_1^0(t)_{\text{field model}} \frac{B(\lambda, \varphi, t)_{\text{indirect}}}{B(\lambda, \varphi, t)_{\text{field model}}} \quad (1)$$

where $\lambda$, $\varphi$, and $t$, represent longitude, latitude and age, respectively; the field models used are gufm1 (Jackson et al., 2000) for 1590 ≤ $t$ ≤ 1990 and CHAOS-5 (Finlay et al., 2015) for 1997 < $t$ ≤ 2015. Then, the absolute intensities of $g_{1\text{indirect}}^0$ from archaeointensity estimates are calculated by:

$$B_{g_{1\text{indirect}}^0} = \frac{g_{1\text{indirect}}^0}{r^3} a^3 (1 + 3\cos^2\theta)^{\frac{1}{2}} \quad (2)$$

where $r$ is the mean Earth radius, $\theta$ is the co-latitude, and $a$ is the mean radial distance from the Earth's center; for Earth surface estimation, we can approximate $a$ by $r$. Finally, for this case, $\text{ADM}_{\text{indirect}}$ are estimated

(in $\times 10^{22} \text{ A m}^2$) by:

$$\text{ADM}_{\text{indirect}} = \frac{4\pi r^3}{\mu_0} B_{g_{1\text{indirect}}^0} (1 + 3\cos^2\theta)^{-\frac{1}{2}} \quad (3)$$

where $\mu_0$ is the permeability of free space. Note that the insertion of Eq. (2) into Eq. (3) eliminates the dependence with the co-latitude, transforming it in a direct relation between $\text{ADM}_{\text{indirect}}$ and $g_1^0$. Virtual axial dipole moments ($\text{VADMs}_{\text{indirect}}$) are calculated (in $\times 10^{22} \text{ A m}^2$) using Eq. (3), replacing $\text{ADM}_{\text{indirect}}$ by $\text{VADM}_{\text{indirect}}$ and $B_{g_{1\text{indirect}}^0}$ by $B_{\text{indirect}}$.

## 2.4. Linear regression applied to the selected dataset

There are several statistical methods to infer the geomagnetic axial dipole variations through time from intensity data (e.g., splines, polynomials, moving averages). We decided to use linear regression in order to simplify the description of the geomagnetic dipole variations, thus providing a common solution across the longest possible period within the last two millennia; and also to correlate the variations described by direct and indirect data, considering their respective resolutions. We justify our parsimonious model on the grounds of four main (general) points: i) linear fits have proved to be sufficient to account for archaeological and historical data within their estimated errors during the historical period (Gubbins et al., 2006; Finlay, 2008); ii) a robust linear fit that describes the dataset taking into account all experimental errors will ignore any rapid, local variations; iii) at the point where the linear model no longer satisfies the dataset, it suggests that there has been a change in the general trend; iv) with a model (mathematical function) it is possible to make quantitative comparisons in relation to the physical models that may describe, for example, core features (e.g. Jackson (2003)). Therefore, if linear regression is statistically satisfied, this may provide us with insights into links between long-term (millennial) geomagnetic dipole variations recorded at Earth's surface and core physical mechanisms.

In this light, we describe the strategy employed in this work as follows. First an expression regarding the variations of the geomagnetic axial dipole as a function of time is defined as:

$$g_1^0(t) = \alpha t + \beta + f(t) \quad (4)$$

where $f(t)$ represents all nonlinear variations of $g_1^0$ as a function of time, and $t$ is the time in years (CE) defined for the interval $0 \le t \le t_{cy}$, where $t_{cy}$ is the current year. Since $f(t)$ is an unknown function and represents the manifestation of several mechanisms operating in the Earth's core, its modeling requires more sophisticated physical/mathematical approach as well as a large number of data. However, if we assumed the hypotheses that short-period variations of $f(t)$ can be minimized through average trends in restricted time windows (e.g., ~15 nT/yr for the last 150 years; Jackson et al., 2000; Finlay et al., 2015), we can represent the variation of the geomagnetic axial dipole as:

$$g_1^0(t_{sp}) \approx \alpha t_{sp} + \beta \quad (5)$$

where $t_{sp}$ represents the time for a sub period between 0 and $t_{cy}$, and $\alpha$ and $\beta$ represent the angular and linear coefficients. In addition, we have $\dot{g}_1^0(t_{sp}) = \alpha$.

From Eq. (5) we can calculate linear regressions for datasets belonging to different time windows. From an appropriated linear regression method, which provides both $\alpha$ and $\beta$ values and their respective uncertainties ($\sigma_{\alpha}$ and $\sigma_{\beta}$), it is possible to obtain a set of $\gamma$ linear solutions for each sub period $t_{sp}$, which can be written as:

$$L = \{\gamma_1, \gamma_2, \dots, \gamma_n\} \quad (6)$$

where each $\gamma \in L$ represents a solution with a specific value of $\alpha \pm \sigma_{\alpha}$ and $\beta \pm \sigma_{\beta}$, and $n$ represents the number of solutions belonging to the set ($n \to \infty$, since $\{\alpha, \beta, \sigma_{\alpha}, \sigma_{\beta}\} \in \mathbb{R}$). Fixing a sub period with a time window ($T_w$) that minimizes $f(t)$, it is possible to perform successive linear regressions in order to cover the entire period between 0 and $t_{cy}$.
Thus we can find a subset of linear solutions (S) given by the inter-section of the largest number of sets L (i.e., $S = L_1 \cap L_2 \cap ... \cap L_m$). Finally, from a priori information of the geomagnetic field, it is feasible to refine the number of linear solution $\gamma \in S$ that represents the linear variation of the geomagnetic axial dipole for the longest period between 0 and $t_{cy}$. The advantages of this approach are: i) it tends to minimize the effect of rapid variations of $f(t)$, thus restricting the scenario of signals from mechanisms that operate in the Earth's core; and ii) it is sufficiently robust to define a period in which the variations re-presented by equation 5 is valid, by using the condition $S \cap L = \emptyset$, i.e., it is sufficient to capture the period in which there is no significant influence of nonlinear variations expressed by $f(t)$ for long periods (millennial scale), further restricting the scenario of physical mechanisms responsible for geomagnetic axial dipole variations as a function of time for the last millennium. In this work the value of $g_l^0$ (2015) from CHAOS 5 (Finlay et al., 2015) will be used to estimate $\beta$ and gufm1 (Jackson et al., 2000) will be used to refine the $\alpha$ value.

To calculate the linear regressions, we employ the following strategy. Initially, for the data belonging to Gauss era (1840–2009 CE) we computed a linear regression and its respective “reduced residual” (RR). The RR is defined by $[y-y_a(x)]/\sqrt{(\sigma_y)^2 + (\alpha\sigma_x)^2}$, where $y_a(x) = \alpha\chi + \beta$ is an interpolated linear function and $\sigma_y$ and $\sigma_x$ are uncertainties of (V)ADM and age, respectively. The fitting parameters were obtained from the dataset into a fixed limit of $\pm 3$ of RR in order to refine the uncertainties of the linear model ($\sigma_\alpha$ and $\sigma_\beta$), assuming that “y” (VADM and ADM) and “x” (age of thermoremanent magnetization of the material) are variables with independent uncertainties (Bevington and Robinson, 2003). The described procedure was repeated for multiple earlier intervals spanning an arbitrary time period such that each contained the same number of data of the Gauss era (48 data). The end dates of each interval were 50 years apart (1959, 1909, etc) but the start date was determined solely by the requirement to have 48 data within the interval. Interval lengths therefore varied from 107 (between 1600 and 1707 CE) to 450 (between 0 and 450 CE) years. It is worth noting that the time window used for each linear regression and the shift by 50 years is sufficient to average typical rapid time-variations of the larger scale part of the geomagnetic field (Hulot and Le Mouël, 1994; Lhuillier et al., 2011).

# 3. Results

## 3.1. Filtered archaeointensity database

We applied the data selection criteria for the time-window between 0 and 2009 CE, comprising all results obtained from archaeological artifacts and volcanic rocks. From a total of 2532 data, only 413 entries from 62 published papers fulfill the selection criteria (Fig. 1; Table A1). These include 290 data from archaeological artifacts and 123 data from volcanic material. Regarding to the initial number of data, the selection criteria that excluded the largest number of entries was the archaeointensity method and protocols adopted (62%). In contrast, the criteria responsible to verify if there was previous investigation regarding MD grains influence did not exclude any data. The criteria about age uncertainties excluded 8% of the initial dataset; the criteria regarding thermo(chemical) alterations during the experiment excluded 3.7%; the criteria about anisotropy and cooling rate corrections excluded 4.6% and 3.2%, respectively; and the criteria that evaluated the final result excluded 2.6% of the initial catalogue (Fig. 1).

Fig. 2 shows the means and medians for the non-filtered (Fig. 2a) and filtered (Fig. 2b) datasets, as well as the age and intensity uncertainties of each data. For the original dataset (Fig. 2a), disregarding the uncertainties of individual measurements (i.e., box-and-whisker plot), there is a greater smoothing of the virtual axial dipole variations for the last 2000 years. This comes as no surprise, since the means and medians were estimated from a cloud of many values within a restricted range. Applying the selection criteria reduces the smoothness of the

<figure><img src="image_3.png" /><figcaption>Fig. 1. Summary of remaining datasets after the application of each selection criterion (see Methods).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Archaeointensity data for the last two millennia. Box-and-whisker plots of 100-years-subsets of the archaeointensity data are represented in blue and black, mean and standard deviation for the same subsets are represented in red, and age-intensity uncertainties of all dataset are represented in light green. a) represents the complete dataset, and b) the filtered dataset. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

variation prior to 1000 CE but enhances how well the trend is defined in the later part of the record whilst maintaining its shape (Fig. 2b).

Fig. 3 shows the spatial (Fig. 3a–c) and temporal (Fig. 3d–f) distribution of the archaeointensity data for the last two millennia. The representations were divided into three distinct periods that will be explored throughout the work: between 1840 and 2009 CE (Gauss era) (Fig. 3a and d), 1590 and 2009 CE (Fig. 3b and e), and 0 and 2009 CE (Fig. 3c and f). Although there is a greater concentration of data in Europe, it is important to note that, even after data selection, the same relative geographic distribution remained (Fig. 3a–c). The temporal
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Spatiotemporal distribution of the archaeointensity database for the last two millennia. Maps (a–c) show the geographic data distribution, where triangles (black) indicate all dataset, and stars (red) only selected data (more details in the main text). Histograms (d–f) show the distribution of age, where gray bars represent all dataset, and red bars the filtered dataset. (a, d) represent the distribution for the Gauss era (after 1840 CE), (b, e) for the “historical period” (after 1590 CE), and (c, f) for the last two millennia. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

distribution of the data provides a good coverage for the three averaged periods, albeit with some peaks in the number of data (e.g., 1940–1960 CE, 1600–1700 CE).

Considering the similarity of the spatial distribution and the good temporal coverage of the selected data for the three periods described above, the results were compared with the gufm1 and CHAOS 5 models (Jackson et al., 2000; Finlay et al., 2015) for the periods 1840–2009 CE (Gauss Era) and 1590–2009 CE. The main objectives of this comparison were to test the latitudinal distribution of the data, since this variation has a direct influence on the magnetic intensity estimates (e.g. Campuzano et al. (2015)), as well as to test the compatibility between the high-quality archaeointensity values and the historical field models (Fig. 4). Fig. 4a and b show the mean VADM-ADM values from archaeointensity data (each point represents the averaged intensity in a latitudinal-degree-spaced) and from the gufm1 (each point is given by the averaged intensity of 36 data in a longitudinal-eqally-spaced distribution, represented every five degrees of latitude). Although there is a greater concentration of data in the northern hemisphere (~72.9% and ~73.2% for Gauss era and 1590–2009 CE, respectively), it is important to note that they are distributed in a range of ~100 degrees (–38 to 64 degrees) of latitude. In addition, considering the uncertainties of the measurements, the difference between VADM and ADM presents values close or equal to zero. Also, almost all data (except one result) show a good correspondence with gufm1 for both temporal averages.

For the same periods, the archaeointensity results were compared to expected values from gufm1 and CHAOS 5 (simplified by $B_{gufm1}$; Jackson et al., 2000; Finlay et al., 2015) and the resulting distributions of the residuals plotted in Fig. 4c and d. In both cases the selected data

show less scatter than the original unfiltered data. A symmetric distribution within one standard deviation of zero was observed for both subsets (–0.58 ± 4.18 and –2.01 ± 5.48 μT), indicating a good concordance between measured data and the historical field models. In order to test all possible scenarios regarding to the differences between archaeointensity and gufm1 data, we repeated the same comparison taking into account experimental and age uncertainties of the indirect results, where a normal distribution with standard deviation covering the zero was also observed (Fig. S1a to S1i). In addition, we performed a Monte Carlo approach by using a homogeneous distribution for the intensity and age uncertainties, and again the averaged residual shows values close to zero for both periods (Fig. S2). This high-quality and historical-model-comparable catalogue is the main basis for our analysis of the temporal variation of the geomagnetic axial dipole.

## 3.2. The geomagnetic ADM for the last four centuries

Filtered intensity estimates for the Gauss era comprise 48 data time-geographically distributed (Fig. 3a and d; Table A1). These data were converted into $ADM_{indirect}$ values from the theorem of Hulot et al. (1997) (see Methods). It is worth noting that this theorem requires a complete geometric coverage of the magnetic field on the globe. Therefore, we computed the $ADM_{indirect}$ values from archaeointensity data using gufm1 (Jackson et al., 2000) for 1840–1990 CE and CHAOS-5 (Finlay et al., 2015) for 1997–2015 CE, since they are well-established field models available for these respective periods. Then a linear regression was performed from converted $ADM_{indirect}$ data, resulting in a set of linear solutions for this particular time interval. In our analyzes, the linear solution set is given by the uncertainties of α and β of a linear
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Latitudinal influence and intensity consistency of the archaeomagnetic dataset when compared with gufm1. (a–b) show the variation of the intensity for different latitudes, where opened circles (gray) represent data from gufm1, given by the averaged intensity of 36 data for every 10 degrees of longitude equally-spaced; and closed circles (red) represent the averaged intensity in a latitudinal-degree-spaced of the filtered archaeointensity data (see Methods). Histograms show the distribution of the difference between the absolute archaeointensity data obtained from laboratory and: (*) gufm1 (Jackson et al., 2000) for 1590–1990 CE, and CHAOS-5 (Finlay et al., 2015) for 1990–2009 CE. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

function type; in our case $y(x) = \alpha x + \beta$ with y representing (virtual) axial dipole moment and x the age of the thermoremanent magnetization of the material (Fig. 5; see Methods). ADM<sub>model</sub> decay rates from gufm1 and CHAOS-5 field models for the Gauss era fall well within our linear solution set (Fig. 6a).

Similarly, for the 1590 and 1840 CE time interval, we converted 81 high-quality archaeointensity data into $ADM_{indirect}$ values (Fig. 3b and e; Table A1), using gufm1 for 1590–1840 CE. Subsequently, six linear regressions were calculated individually using the same number of data for the Gauss era, every 50 years before 2009 CE (Fig. 6a, see Methods). All solution sets comprise the $ADM_{model}$ decay rate from gufm1. Between 1600 CE and 1800 CE the sets of linear solutions have greater uncertainties than those obtained for more recent periods. This is due to the sensitivity of linear regressions to uncertainties of age and axial dipole moment. For example, Schnepp et al. (2009) presented 25 archaeointensity results from oven floors collected in Germany. From these, 10 sites dated for 1665 ± 85 AD show intensity values over a wide range from 59.7 ± 2.9 to 44.1 ± 3.3. So, despite the excellent quality of individual archaeointensity estimates their relatively high age uncertainties and range in the intensity results strongly influence our regressions. Notwithstanding, the removal of these data reduces the uncertainties of regressions, but does not change any features or trends. For this reason, we decided to keep them into our analyses.

When we consider the whole filtered archaeointensity dataset comprising the past four centuries, the data are normally distributed, but the mean is slightly offset relative to intensity estimates of the gufm1 model (Fig. 4d). The same behavior is observed when this archaeointensity catalogue is converted into $VADM_{indirect}$ values and

compared with gufm1 (Fig. 6b). It is important to note that several studies have emphasized poor accuracy at the values calculated by gufm1, especially for the pre-Gauss era (e.g. Le Goff and Gallet (2017)). However, the comparisons performed in this work between high-quality archaeointensity data and those calculated by gufm1 present, on average, a satisfactory correspondence (Figs. 4c, d and 6b), since the mean residual covers the zero value within one standard deviation. Although there is a need to generate models with greater accuracy in the calculation of the complete vector of the Earth's magnetic field, the gufm1 remains the most robust full-vector magnetic field model for the historical period. Thus, we use the gufm1 and CHAOS-5 for the Gauss era and gufm1 for the period 1590–1990 CE in an attempt to capture a linear trend exclusively from archaeointensity data, which represents the average variation of the axial dipole for the historical period; and also minimizes the slightly offset of the mean residual between indirect data and model.

The mean linear trend given by the slope from the gufm1 and CHAOS-5 models (i.e., ~15 nT/yr) is one of several possible solutions that belong to the set of linear solutions between 1590 and 2009 CE (i.e., $\gamma_{gufm1 and CHAOS-5} \in S$) (Fig. 6a). However, there are many other slopes that satisfy the condition $\gamma \in S$ for this period. To determine a single linear solution that better represents the average decrease of the geomagnetic axial dipole, we firstly: i) fix the ADM for 2015 CE at $7.61 \times 10^{22} \text{ A m}^2$ (Finlay et al., 2015), in order to restrict the $\beta$ value in the Eq. (5); and ii) find a slope ($\gamma$) that minimizes the difference between absolute data from laboratory and model in order to obtain a normal distribution of the residual centered on zero value. To obtain (ii), we replaced the $g_1^0$ coefficient of gufm1 by those extracted from linear solutions trends, and then we recalibrated all coefficients
<figure><img src="image_7.png" /><figcaption>Fig. 5. Linear regression of the high-quality archaeointensity dataset for the Gauss era. Up graphic shows the filtered intensities converted into ADM<sub>indirect</sub> values as a function of age of the natural thermoremanent magnetization. Green shaded area represents the set of linear solutions obtained from the dataset into a limit of ± 3 RR, as represented by gray area in the below graphic (see Methods). Dashed and continuous green lines represent the average and extremes of the linear solution set, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

following the strategy showed by Whaler and Holme (2011):

$$g_{lnew}^{m} = g_{lgufm1}^{m} g_{larch}^{0} / g_{lgufm1}^{0} \qquad (7)$$

and

$$h_{lnew}^{m} = h_{lgufm1}^{m} g_{larch}^{0} / g_{lgufm1}^{0} \qquad (8)$$

where l and m represent degree and order, respectively. Finally, the best linear fit found is indicated as a blue continuous line in Fig. 6a (relative to the residual represented with the same color in Fig. 6b) and corresponds to an intensity decay rate of 12.5 nT/yr.

In order to statistically test the slope of 12.5 nT/yr obtained for the period between 1590 and 2009 CE, we performed 100 simulations in which 60% of 129 archaeointensity data were randomly selected, and then we calculated the 95% bootstrapping confidence intervals (N = 1999) (Hammer et al., 2001) for each linear regression, and also evaluated the residuals between the mean linear fit and i) the data used to compute the linear fit; and ii) the data that were not used to compute the linear fit (Fig. S3). Then we observe that the slope obtained above is statistically robust (Fig. S3a) and that both residuals show a normal distribution within one standard deviation of zero (Fig. S3b). Therefore, the linear solution obtained exclusively by archaeointensity data apparently emerges as a robust solution to describe the average ADM decay trend for the historical period.

### 3.3. The geomagnetic ADM for the last two millennia: VADM application as a proxy

Does the linear trend obtained for the past four centuries also describe the geomagnetic axial dipole further back in time? To address this question, we first need to test whether we are able to safely use VADM as a “proxy” of ADM.

<figure><img src="image_8.png" /><figcaption>Fig. 6. Continuous decrease of the geomagnetic axial dipole for the last four centuries. (a) Seven linear solutions sets for the last four centuries are represented in light green shaded area; continuous green lines represent the extremes of each linear solution set. Yellow and red lines represent the ADM of CHAOS-5 and gufm1, respectively. Blue line represents the common linear solution for the sets of linear regressions, obtained exclusively from archaeointensity data. (b) Normal distributions of the difference between VADMs<sub>indirect</sub> from archaeointensity data and VADMs<sub>model</sub> from field models: gufm1 (red bars) and gufm1 with its original axial dipole coefficient replaced (and the others coefficients recalibrated) by the linear trend proposed here (blue 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>

In order to estimate the error arising from our use of VADMs<sub>indirect</sub> in the pre-Gauss era instead of values of ADMs<sub>indirect</sub> calculated from theorem of Hulot et al. (1997), we compare the distributions of ADMs (both those calculated using Hulot et al.’s theorem from archaeointensity data and those taken directly from gufm1 and CHAOS-5) and VADMs (both those taken from archaeointensity data and from gufm1 and CHAOS-5), for the periods 1840–1990 CE, 1590–1840 CE and 1590–1990 CE (Fig. 7). In all cases, the residuals distributions between VADM and ADM are centered close to zero and well-inside of one standard deviation. We concluded that, although the (V)ADMs<sub>indirect</sub> are marginally offset to lower values than the ADMs<sub>model</sub> (Fig. 8a), they nevertheless provide a useful proxy for the Gauss era, as well as for the entire historical period.

In order to test whether the geographical distribution of the entire dataset used here is sufficient to define ADMs over the last two millennia, we segmented this into 150 years intervals (taken as the same duration as the Gauss era time interval belonging to gufm1 – 1840–1990 CE), and then for each sub-period we re-sampled VADMs<sub>model</sub> from gufm1 (1840–1990 CE) taking into account the longitudes and latitudes of each of the archaeointensity data from the preceding time intervals (Fig. S4). Afterwards we calculated the normal distribution of the difference between the re-sampled VADMs<sub>model</sub> and the ADMs<sub>model</sub> from gufm1 (Fig. S4). For all segments, the obtained
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Correlations between ADMs and VADMs from archaeointensity data and gufm1. Normal distribution of the difference between VADMs and ADMs for: (a, d, g) only archaeointensity data, (b, e, h) only gufm1 with the same time-space coverage, and (c, f, i) archaeointensity data (VADM) and gufm1. The comparisons were separated into the periods: (a, b, c) Gauss era, (d, e, f) 1590–1840 CE, and (g, h, i) 1590–1990 CE.</figcaption></figure>

peak of the VADM-ADM distributions was close to and within one standard deviation of the zero value. This suggests that, presuming variability with the Gauss era is reasonably representative of that for the last 2000 years, the geographical distribution of VADM<sub>indirect</sub> data in the intervals prior to 1840 would be sufficient to adequately describe the main trend of ADM in the post-1840 interval and that the limited geographical distribution of VADM<sub>indirect</sub> estimates prior to 1840 should not be a barrier to defining ADM variations.

Additionally, we show that the trend of (V)ADM<sub>indirect</sub> from the archaeointensities in the gufm1 time period is close to that of the ADM<sub>model</sub> taken from gufm1 and falls within one standard deviation of the VADM<sub>model</sub> values from gufm1 in a degree-spaced coverage around the globe for each year (64,800 estimations per year) (Fig. 8a). Finally, we tested the influence of the linear trend from this study (Section 3.2) in the original gufm1 (Jackson et al., 2000) by placing it into the gufm1 (Eqs. (7) and (8)), and then calculating the average of the radial component of the magnetic intensity at the core-mantle boundary (CMB) (Fig. 8b–f). The average field at the CMB shows similar geometry and maximum differences of 0.012 mT for Gauss era (Sup. Fig. 8d) and 0.016 mT for 1590–1990 CE (Fig. 8g). These comparisons converge in showing that VADM<sub>indirect</sub> is expected to be a good proxy for ADM<sub>indirect</sub> suitable for ascertaining variations on timescales longer than 150 years in the axial dipole over the last 2000 years.

We therefore proceeded back in time with the linear regressions using groups of 48 data (i.e., the same number used for the Gauss era), every 50 years, over the pre-gufm1 period (Fig. 9). The 26 sets of linear regressions include a total of 275 high-quality archaeointensity data (Table A1), belonging to the range of 0–1590 CE. Linearly extrapolating

the previous linear trend back 1590 CE we found that it is a common solution (that satisfy $\gamma \in S$) for 18 sets of solutions, being consistent until the period between 550 CE and 750 CE (Fig. 9). The period 550–750 CE marks the interval of two consecutive sets of linear solutions where the linear fit is a solution for the last time (550–800 CE) and where it fails for the first time (475–750 CE). Thus, for the last millennium, the single linear regression with slope 12.5 nT ($\sim -0.0032 \times 10^{22} \text{ A m}^2$), intercept $7.61 \times 10^{22} \text{ A m}^2$ at 2015 CE (Finlay et al., 2015) and valid for the period 750–2015 CE (hereafter called by archaeo_adm1.3k), appears to be a useful description of the long-term variations of the geomagnetic axial dipole and one that differs substantially from existing models in the interval 750–1200 CE (Fig. 10).

# 4. Discussion

## 4.1. Effects of small spatial and temporal variations

Some studies have consistently reported rapid local variations of the geomagnetic field for the past two millennia. For example, de Groot et al. (2013) studied lava flow sequences from Hawaii and reported a rapid increase in geomagnetic field strength of about 15 µT between ~850 and 925 CE (~200 nT/yr), followed by a rapid decrease at ~1150 CE (~190/yr). For Western Europe, Genevey et al. (2016) and Gómez-Paccard et al. (2016) used high-quality archaeointensity data to argue for a rapid decrease in geomagnetic field strength of ~100 nT/yr between 800 and 1050 CE. Rapid variations have also been reported for older periods, including dramatic field intensity spikes in southern
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** VADMs and ADMs from archaeointensity data and gufm1, and implications for the intensity variations at CMB. (a) Pink line and area represent averages and standard deviations, respectively, of VADMs<sub>model</sub> values produced from intensity measurements taken from gufm1 (which is constrained by intensity data only since 1840 CE), with a degree-spaced coverage around the globe for each year between 1590 and 1990 CE. Red line represents the ADM<sub>model</sub> of gufm1. Blue line represents the linear trend obtained exclusively from archaeointensity data. Note that the maximum discrepancy between this study's ADM and the gufm1 in the Gauss era (for which the latter is well-constrained) is $0.15 \times 10^{22} \text{ A m}^2$ (~1.8% of the total). (b), (c), (e) and (f) represent the average radial magnetic field intensity for 1840–1990 CE and 1590–1990 CE periods. (b) and (e) were computed using gufm1. (c) and (f) were computed using gufm1 with its original axial dipole coefficient replaced (and the others coefficients recalibrated) by the linear trend proposed here. (d) and (g) provide the residuals between the averages of radial intensity at CMB calculated from gufm1 and the modified gufm1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Israel reported by Shaar et al. (2011, 2016). These rapid variations have been attributed to local anomalies caused by dynamic processes at the CMB (e.g. Livermore et al. (2014), Davies and Constable (2017)). A complete assessment of these variations would require a more complete coverage of the globe with high-quality archaeointensity data, particularly in the southern hemisphere where the field may be more time-dependent (Constable et al., 2016) but which is underrepresented in the archaeomagnetic database. Our analysis tends to eliminate rapid local variations, thus describing only the long-term variations of the axial dipole field strength.

Our analysis also tends to average out the small amplitude variations of ADM<sub>model</sub> derived from observatory and satellite data that are taken into account in, for example, gufm1 and CHAOS-5 models (Jackson et al., 2000; Finlay et al., 2015; Figs. 6a and 9). Some studies reported that regional high-quality archaeointensity data have

sufficient resolution to suggest oscillatory behavior of the geomagnetic axial dipole (Genevey et al., 2009; Hartmann et al., 2011). However, given the inherent experimental errors and limited geographical coverage, we suggest that the current global archaeointensity dataset, on average, cannot reproduce small fluctuations in the geomagnetic axial dipole.

## 4.2. Comparison between archaeo_adm1.3k and geomagnetic field models

When the archaeo_adm1.3k is compared with historical models, some differences can be related to data quality instead of Earth's core dynamics. Between 1590 and 1840 CE some models described the ADM evolution by a linear trend. For example, gufm1 (Jackson et al., 2000) used a linear extrapolation from the main trend of the Gauss era, Gubbins et al. (2006) proposed a linear regression from 315 non-
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Continuous decrease of the geomagnetic axial dipole for the last millennium. Twenty-six linear solutions sets for the period before 1590 CE are represented in pink and light blue shaded area; continuous pink and light blue lines represent the extremes of each linear solution set. Yellow and red lines represent the ADM<sub>model</sub> of CHAOS-5 and gufm1, respectively. Continuous blue line represents the linear trend obtained exclusively from archaeointensity data (this study); and dashed blue line represents it linear extrapolation. The period indicated by the gray area marks the interval of the last linear regression in which the linear extrapolation is a solution and the first that it fails. Just to facilitate visualization, the sets of linear regressions represented in light blue are those that contain the linear trend as the solution, and the sets represented by pink are those that do not contain it. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Comparison between the continuous linear trend decreases (archaeo_adm1.3k) with (a) geomagnetic field models and (b) VADM curves computed using temporal and spatial averaging. For the CALS3k series, A_FM-M and ASDI_FM-M dashed lines represent the original ADM<sub>model</sub>, and continuous line the smoothed trend. Dashed lines for Gubbins et al.'s model represent the uncertainties. Dashed line for archaeo_adm1.3k represents the interval that the linear trend is a solution for the last time and fails for the first time (details in the main text). The arrow indicates the possible trend before 550–750 CE (see also Fig. 9).</figcaption></figure>

filtered archaeointensity data, and Finlay (2008) suggested a linear trend given by the best fit from inversion of the same 315 archaeointensity data (Fig. 10a). The models of Gubbins et al. (2006) and Finlay (2008) suggested a shift in decay trend at 1840 CE, which is a recurrent feature of the differences between direct and indirect estimates (Suttie et al., 2011). Here, using only high-quality archaeointensity data to describe the ADM time-evolution, the shift in decay at 1840 CE is suppressed. Our estimates of ADM before 1840 CE differ from the flat evolution proposed by Gubbins et al. (2006) and Finlay (2008). Instead, we suggest an earlier start in the decay of the ADM.

Greater time-period models obtained by different data sets and modeling strategies also show significant differences regarding archaeo_adm1.3k (Fig. 10a). We compare our main result with the CALS3k series of field models (Korte et al., 2009; Korte and Constable, 2011), and with the A_FM-M and ASDI_FM-M models, which are the mean models of ensemble of time-varying archaeomagnetic field models (Licht et al., 2013). The CALS3k and the ASDI_FM-M field

models are constructed from archaeological, volcanic and sedimentary data, whereas A_FM-M only uses archaeointensity data. For better visualization, all millennial models used in our comparison had their ADM' curves smoothed using a larger smooth factor in a sequence of third-order polynomials continuous up to the second derivative (de Boor, 1978) (Fig. 10a). Besides differences in the main trend between the archaeo_adm1.3k and the mentioned models, in the specific period 550–900 CE the CALS3k models present a low of the ADM, being opposite to the peak described here. The Licht et al.,' models show intermediate values of ADM for this period. This may be a consequence of the field models inappropriately partitioning energy into higher order components and/or their inclusion of less reliable intensity data.

Other descriptions regarding the variations of the geomagnetic axial dipole were presented through VADM curves computed using temporal and spatial averaging (e.g. Valet et al. (2008), Genevey et al. (2008), Knudsen et al. (2008), Usoskin et al. (2016)) (Fig. 10b). From a non-filtered archaeointensity database, Valet et al. (2008) proposed a third-
degree polynomial function model to describe the variation of the (V)ADM<sub>indirect</sub> for the last 2000 years, from a running-window approach using time-averaged data over 100 years and shifted by 25 years. They argued that time-averaged windows of 100 years are enough to attenuate non-dipolar contributions. The result from Valet et al.'s work shows a peak of (V)ADM values at about 700 CE, followed by a decrease up to the present. Knudsen et al. (2008) adopted a running-window approach to calculate the VADM variations during the entire Holocene. They presented a least-square fit from time-averaged data using 500 year long sliding windows shifted by 100 years (between –2000 and 2000 CE). The description put forward by Knudsen et al.'s VADM curve shows a similar behavior to that described in Valet et al.'s paper, with relatively lower values beginning from a peak in the VADM at about 700 CE. Genevey et al. (2008) proposed a VADM evolution for the last 10,000 years. Their results are given by time-averaged data using 200 year long sliding windows shifted by 100 years (between –1000 and 2000 CE), and show a peak in VADM values between 300 and 400 CE, and the beginning of the VADM decrease at about 1200 CE. Similarly, Usoskin et al. (2016) presented a description of VADM from time-averaged data using 200 year long sliding windows shifted by 10 years (between –1500 and 2000 CE), using the newly updated database GEOMAGIA50.v3 (Brown et al., 2015). Usoskin et al.'s work provides a VADM curve with variations similar to that presented by Genevey et al. (2008), but with relatively lower mean values for the period 0–1440 CE. For the last nine centuries, our results, which were obtained through regressions by different time-windows with each individual endpoints regularly spaced by 50 years, are in agreement with that proposed by Valet et al.'s, Genevey et al.'s, Knudsen et al., and Usoskin et al.'s works, and emerge as a simplification of all them. Before 900 CE the archaeo_adm1.3k presents a unique trend that started between 550 and 750 CE. This reinforces that, for the last two millennia, the current models tend to lose information about the variation of the geomagnetic axial dipole, even if it is not possible to state specifically the cause.

A promising approach to combine short and long-term changes in axial dipole intensity into a complete description of the geodynamo has been put forward by Sanchez et al. (2016). These authors constructed a geomagnetic field model for the last three millennia using non-filtered archaeomagnetic data and prior information from geodynamo simulations (Aubert et al., 2013). Their geomagnetic axial dipole presents an average decay of ~7 nT/yr for the last millennium, which differs from the result obtained here (12.5 nT/yr). In view of the new methodological approach presented by Sanchez et al' model, we tentatively suggest that the use of only high-quality archaeointensity data may improve attempts to describe physical processes of the Earth's core, which drive global millennial features of the field.

### 4.3. Implications for core dynamics

The average rate of geomagnetic axial dipole decay for the last 184 years is ~15 nT/yr (Jackson et al., 2000; Finlay et al., 2015). For this period, Finlay et al. (2016) combined geomagnetic field models (Gillet et al., 2013) and equatorially symmetric core flow models (e.g. Pais and Jault (2008), Amit and Pais (2013), Aubert (2014), Gillet et al. (2015)) to attribute the axial dipole decay to symmetry breaking in advection sources in the Southern Hemisphere. They showed that the

drift of an intense normal polarity flux path equatorward, which diminishes the ADM (e.g. Olson and Amit (2006)), is unbalanced by any other significant advection source, causing the ADM decrease. According to our analysis, the intensity of the Earth's magnetic axial dipole had an average decay rate of 12.5 nT/yr from ~750 CE to present, beginning after a clear change in the trend of the geomagnetic axial dipole (Fig. 9). Therefore, we have provided evidence for a continuous linear decay comparable in order of magnitude with the average rate for the Gauss era (e.g. Jackson et al. (2000), Finlay et al. (2015), Finlay et al. (2016)). Consequently, we suggest an early break in the symmetry of the ADM advection sources in the Earth's core at about 750 CE.

According to our analysis, the ADM has been decreasing at roughly the rate of present-day for the past ~1265 years. This corresponds to estimations of axial dipole secular variation (SV) time-scales (~1000 years) recently observed by Amit et al. (in press). This SV timescale represents the reorganization time of the axial dipole (Hulot and Le Mouël, 1994). Based on the similarity between our estimate of ADM decrease period and the axial dipole SV timescale (Amit et al., in press), we speculate that the ADM decrease may reach its end soon.

## 5. Conclusions

From a careful analysis of the current archaeointensity dataset we propose a well defined linear trend that describes the variation of the geomagnetic axial dipole for the last millennium: archaeo_adm1.3k. The main conclusions of this work are:

(i) The comparison between the data obtained by direct and indirect measurements during the Gauss era allowed the analysis of the axial dipole to be extrapolated back in time for two millennia.

(ii) The shift in the trend of the geomagnetic axial dipole variation at 1840 CE described by previous studies (Gubbins et al., 2006; Finlay, 2008) is a biased feature of the difference between direct and indirect measurements;

(iii) Considering the last 2000 years, at approximately 750 CE there was a peak of intensity of the axial dipole followed by a quasi-constant decrease, which is not captured by millennial models.

(iv) If the recent decay of the Earth's magnetic axial dipole is caused by asymmetry in the advective sources then this commenced within the interval 550–750 CE.

(v) Comparable duration of the dipole decay and the ADM SV timescale suggests that this event may be reaching its end soon, and therefore tend to disfavor the hypothesis of an imminent geomagnetic reversal.

## Acknowledgements

Richard Holme and Hagay Amit are thanked for pre-submission readings and discussions that helped in improving the paper. We thank two anonymous reviewers for their helpful and constructive comments. W.P. thanks the São Paulo Research Foundation (FAPESP) (grant #2013/16382-0). A.J.B. acknowledges a NERC standard grant (NE/ P00170X/1) and funding from The Leverhulme Trust. G.A.H. thanks CAPES (grant AUXPE 2043/2014) and CNPq (grant 454609/2014-0). F.T-N. thanks The National Council for Scientific and Technological Development (CNPq/Brazil) (grant 206997/2014-0).

## Appendix A

In our analysis we used a catalogue of archaeointensities, which was approved by the data selection criteria adopted, from the GEOMAGIA50.v3 database and from some published studies that has not been inserted into it so far (Gallet et al., 2009; de Groot et al., 2013; Cai et al., 2014; Di Chiara et al., 2014; Kissel et al., 2015; Goguitchaichvili et al., 2015; Osete et al., 2015; Roperch et al., 2015; Shaar et al., 2015; Tarduno et al., 2015; Cai et al., 2016; Genevey et al., 2016; Gómez-Paccard et al., 2016; Poletti et al., 2016; Genevey et al., in press; Salnaia et al., 2017; Shaar et al., 2017) (Table A1).
<table><caption>**Table A1**
Selected archaeointensity dataset. “a” represents archaeological and “v” represents volcanic.</caption><thead><tr><th>Number of data</th><th>Age interval</th><th>Material</th><th>Reference ID</th></tr></thead><tbody><tr><td>1</td><td>0-0</td><td>a</td><td>1</td></tr><tr><td>2</td><td>0-1300</td><td>a</td><td>2</td></tr><tr><td>18</td><td>0-1700</td><td>a</td><td>3</td></tr><tr><td>4</td><td>1030-1546</td><td>a</td><td>4</td></tr><tr><td>9</td><td>1103-1696</td><td>a</td><td>5</td></tr><tr><td>6</td><td>1150-0</td><td>a</td><td>6</td></tr><tr><td>2</td><td>1160-1390</td><td>a</td><td>7</td></tr><tr><td>13</td><td>123-1320</td><td>a</td><td>8</td></tr><tr><td>35</td><td>1237-1832</td><td>a</td><td>9</td></tr><tr><td>18</td><td>1331-1665</td><td>a</td><td>10</td></tr><tr><td>5</td><td>140-1840</td><td>v</td><td>11</td></tr><tr><td>6</td><td>1440-1835</td><td>v</td><td>12</td></tr><tr><td>1</td><td>1450-0</td><td>a</td><td>13</td></tr><tr><td>4</td><td>1550-1750</td><td>a</td><td>14</td></tr><tr><td>14</td><td>1576-1910</td><td>a</td><td>15</td></tr><tr><td>9</td><td>1591-1909</td><td>a</td><td>16</td></tr><tr><td>7</td><td>1610-1850</td><td>a</td><td>17</td></tr><tr><td>8</td><td>1610-1950</td><td>v</td><td>18</td></tr><tr><td>1</td><td>1615-0</td><td>a</td><td>19</td></tr><tr><td>3</td><td>1672-1691</td><td>a</td><td>20</td></tr><tr><td>1</td><td>1706-0</td><td>v</td><td>21</td></tr><tr><td>1</td><td>1766-0</td><td>v</td><td>22</td></tr><tr><td>1</td><td>1766-0</td><td>v</td><td>23</td></tr><tr><td>1</td><td>1790-0</td><td>a</td><td>24</td></tr><tr><td>1</td><td>1835-0</td><td>a</td><td>25</td></tr><tr><td>1</td><td>1886-0</td><td>v</td><td>26</td></tr><tr><td>1</td><td>1955-0</td><td>v</td><td>27</td></tr><tr><td>1</td><td>1959-0</td><td>a</td><td>28</td></tr><tr><td>3</td><td>1960-1982</td><td>v</td><td>29</td></tr><tr><td>2</td><td>1980-2000</td><td>v</td><td>30</td></tr><tr><td>2</td><td>1993-1998</td><td>v</td><td>31</td></tr><tr><td>5</td><td>220-1138</td><td>a</td><td>32</td></tr><tr><td>15</td><td>235-1959</td><td>a</td><td>33</td></tr><tr><td>6</td><td>260-1690</td><td>v</td><td>34</td></tr><tr><td>13</td><td>27-899</td><td>a</td><td>35</td></tr><tr><td>2</td><td>316-601</td><td>a</td><td>36</td></tr><tr><td>2</td><td>330-0</td><td>a/v</td><td>37</td></tr><tr><td>14</td><td>337-1575</td><td>a</td><td>38</td></tr><tr><td>3</td><td>350-1855</td><td>a</td><td>39</td></tr><tr><td>7</td><td>360-1380</td><td>v</td><td>40</td></tr><tr><td>5</td><td>370-605</td><td>a</td><td>41</td></tr><tr><td>3</td><td>37-100</td><td>a</td><td>42</td></tr><tr><td>4</td><td>385-1975</td><td>v</td><td>43</td></tr><tr><td>10</td><td>388-424</td><td>a</td><td>44</td></tr><tr><td>5</td><td>450-810</td><td>a</td><td>45</td></tr><tr><td>9</td><td>480-1660</td><td>v</td><td>46</td></tr><tr><td>21</td><td>550-875</td><td>a</td><td>47</td></tr><tr><td>18</td><td>588-2009</td><td>v</td><td>48</td></tr><tr><td>9</td><td>590-1950</td><td>v</td><td>49</td></tr><tr><td>8</td><td>622-1180</td><td>a</td><td>50</td></tr><tr><td>4</td><td>691-1320</td><td>a</td><td>51</td></tr><tr><td>16</td><td>726-1950</td><td>v</td><td>52</td></tr><tr><td>4</td><td>765-1779</td><td>v</td><td>53</td></tr><tr><td>1</td><td>775-0</td><td>a</td><td>54</td></tr><tr><td>5</td><td>800-1563</td><td>v</td><td>55</td></tr><tr><td>12</td><td>815-1797</td><td>a</td><td>56</td></tr><tr><td>1</td><td>825-0</td><td>a</td><td>57</td></tr><tr><td>8</td><td>850-1088</td><td>a</td><td>58</td></tr><tr><td>5</td><td>852-1350</td><td>a</td><td>59</td></tr><tr><td>11</td><td>87-1770</td><td>v</td><td>60</td></tr><tr><td>2</td><td>934-1783</td><td>v</td><td>61</td></tr><tr><td>4</td><td>990-1450</td><td>a</td><td>62</td></tr></tbody></table>

Reference ID: 1 – Tema et al. (2012); 2 – De Marco et al. (2008); 3 – Chauvin et al. (2000); 4 – Tarduno et al. (2015); 5 – Salnaia et al. (2017); 6 – Gómez-Paccard et al. (2006a); 7 – Stark et al. (2010); 8 – Cai et al. (2014); 9 – Genevey et al. (2009); 10 – Schnepp et al. (2009); 11 – Tanaka and Kono (1991); 12 – Roperch et al. (2015); 13 – Shaar et al. (2017); 14 – Hartmann et al. (2009); 15 – Hartmann et al. (2010); 16 – Hartmann et al. (2011); 17 – Gallet et al. (2005); 18 – Bowles et al. (2005); 19 – Osete et al. (2015); 20 – Poletti et al. (2016); 21 – Kissel et al. (2015); 22 – Gratton et al. (2005a); 23 – Conte-Fasano et al. (2006); 24 – Goguitchaichvili et al. (2015); 25 – Gómez-Paccard et al. (2006b); 26 – Tanaka et al. (2009); 27 – Cottrell and Tarduno (1999); 28 – Catanzariti et al. (2008); 29 – Chauvin et al. (2005); 30 – Michalk et al. (2008); 31 – Carlut and Kent (2000); 32 – Genevey et al. (2003); 33 – Gómez-Paccard et al. (2008); 34 – Mankinen and Champion (1993); 35 – Mitra et al. (2013); 36 – Catanzariti et al. (2012); 37 – Böhnel et al. (2003); 38 – Genevey and Gallet (2002); 39 – Tema et al. (2013); 40 – Pressling et al. (2007); 41 – Genevey et al. (2017); 42 – Gómez-Paccard et al. (2013); 43 – Pick and Tauxe (1993); 44 – Shaar et al. (2015); 45 – Fanjat et al. (2013); 46 – Pressling et al. (2006); 47 – Genevey et al. (2016); 48 – De Groot et al. (2013); 49 – Yoshihara et al. (2003); 50 – Gómez-Paccard et al. (2012b); 51 – Cai et al. (2016); 52 – Spassov et al. (2010); 53 – Yu (2012); 54 – Donadini et al. (2008); 55 – Di Chiara et al. (2014); 56 – Genevey et al. (2013); 57 – Gallet et al. (2009); 58 – Gómez-Paccard et al. (2016); 59 – Gómez-Paccard et al. (2012a); 60 – Gratton et al. (2005b); 61 – Stanton et al. (2011); 62 – Goguitchaichvili et al. (2011).
## Appendix B. Supplementary data

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.pepi.2017.11.005.

## References

Aitken, M.J., Alcock, P.A., Bussel, G.D., Shaw, C.J., 1981. Archaeomagnetic determination of the past geomagnetic intensity using ancient ceramics: allowance for anisotropy. Archaeometry 23, 53–64.

Aitken, M.J., Allsop, A.L., Bussel, G.D., Winter, M.B., 1988. Determination of the intensity of the Earth's magnetic field during archaeological times: Reliability of the Thellier technique. Rev. Geophys. 26, 3–12.

Amit, H., Pais, M.A., 2013. Differences between tangential geostrophy and columnar flow. Geophys. J. Int. 194, 145–157.

Amit, H., Coutelier, M., Christensen, U.R., 2017; al, in press. On equatorially symmetric and antisymmetric geomagnetic secular variation timescales. Phys. Earth Planet. Inter.

Amit, H., Korte, M., Aubert, J., Constable, C., Hulot, G., 2011. The time-dependence of intense archaeomagnetic flux patches. J. Geophys. Res. 116, B12.

Arneitz, P., Egli, R., Leonhardt, R., 2017a. Unbiased analysis of geomagnetic data sets and comparison of historical data with paleomagnetic and archaeomagnetic records. Rev. Geophys. 55, 5–39.

Arneitz, P., Leonhardt, R., Schneppe, E., Heilig, B., Mayrhofer, F., Kovacs, P., Hejda, P., Valach, F., Vadasz, G., Hammerl, C., Egli, R., Fabian, K., Kompein, N., 2017b. The HISTMAG database: combining historical, archaeomagnetic and volcanic data. Geophys. J. Int. 210, 1347–1359.

Aubert, J., 2014. Earth's core internal dynamics 1840–2010 imaged by inverse geodynamo modelling. Geophys. J. Int. 197, 119–134.

Aubert, J., Finlay, C.C., Fournier, A., 2013. Bottom-up control of geomagnetic secular variation by the Earth's inner core. Nature 502, 219–223.

Bevington, P.R., Robinson, D.K., 2003. Data Reduction and Error Analysis for the Physical Sciences. McGraw-Hill.

Biggin, A.J., Badejo, S., Hodgson, E., Muxworthy, A., Shaw, J., Dekkers, M.J., 2013. The effect of cooling-rate on the intensity of thermoremanent magnetization (TRM) acquired by assemblages of pseudo-single domain, multi domain, and interacting single domain grains. Geophys. J. Int. 193, 1293–1294.

Böhnel, H., Biggin, A.J., Walton, D., Shaw, J., Share, J.A., 2003. Microwave palaeointensities from a recent Mexican lava flow, baked sediments and reheated pottery. Earth Planet. Sci. Lett. 214, 221–236.

Bowles, J., Gee, J.S., Kent, D.V., Bergmanis, E., Sinton, J., 2005. Cooling rate effects on paleointensity estimates in submarine basaltic glass and implications for dating young flows. Geochem. Geophys. Geosyst. 6, Q07002.

Brown, M.C., Donadini, F., Korte, M., Nilsson, A., Korhonen, K., Lodge, A., Lengyel, S.N., Constable, C.G., 2015. GEOMAGIA50. v3: 1. general structure and modifications to the archaeological and volcanic database. Earth Planets Space 67, 83.

Cai, S., Jin, G., Tauxe, L., Deng, C., Qin, H., Pan, Y., Zhu, R., 2016. Archaeointensity results spanning the past 6 kiloyears from eastern China and implications for extreme behaviors of the geomagnetic field. Proc. Natl. Acad. Sci. 114, 39–44.

Cai, S., Tauxe, L., Deng, C., Pan, Y., Jin, G., Zheng, J., Xie, F., Qin, H., Zhu, R., 2014. Geomagnetic intensity variations for the past 8 kyr: New archaeointensity results from Eastern China. Earth Planet. Sci. Lett. 392, 217–229.

Campuzano, S.A., Pavón-Carrasco, F.J., Osete, M.L., 2015. Non-Dipole and Regional Effects on the Geomagnetic Dipole Moment Estimation. Pure Appl. Geophys. 172, 91–107.

Carlut, J., Kent, D.V., 2000. Paleointensity record in zero-age submarine basalt glasses: Testing a new dating technique for recent MORBs. Earth Planet. Sci. Lett. 183, 389–401.

Casas, L., Shaw, J., Gich, M., Share, J.A., 2005. High-quality microwave archaeointensity determinations from an early 18th century AD English brick kiln. Geophys. J. Int. 161, 653–661.

Catanzariti, G., Gómez-Paccard, M., McIntosh, G., Pavón-Carrasco, F.J., Chauvin, A., Osete, M.L., 2012. New archaeomagnetic data recovered from the study of Roman and Visigothic remains from central Spain (3rd-7th centuries). Geophys. J. Int. 188, 979–993.

Catanzariti, G., McIntosh, G., Gómez-Paccard, M., Ruiz-Martínez, V.C., Osete, M.L., Chauvin, A., Scientific Team, A.A.R.C.H., 2008. Quality control of archaeomagnetic determination using a modern kiln with a complex NRM. Phys. Chem. Earth. 33, 427–437.

Chauvin, A., Garcia, Y., Lanos, Ph., Laubenheimer, F., 2000. Palaeointensity geomagnetic field recovered on archaeomagnetic sites from France. Phys. Earth Planet. Int. 120, 111–136.

Chauvin, A., Roperch, P., Levi, S., 2005. Reliability of geomagnetic paleointensity data: the effects of the NRM fraction and concave-up behavior on paleointensity determinations by the Thellier method. Phys. Earth Planet. Int. 150, 265–286.

Coe, R.S., 1967. The determination of paleointensities of the Earth's magnetic field with emphasis on mechanisms which could cause non-ideal behavior in Thellier's method. J. Geomag. Geoelectr. 19, 157–179.

Coe, R.S., Grommé, C.S., Mankinen, E.A., 1978. Geomagnetic paleointensities from radiocarbon dated lava flows on Hawaii and the question of the Pacific non-dipole low. J. Geophys. Res. 83, 1740–1756.

Constable, C., Korte, M., Panovska, S., 2016. Persistent high paleosecular variation activity in southern hemisphere for at least 10 000 years. Earth Planet. Sci. Lett. 453, 78–86.

Conte-Fasano, G., Urrutia-Fucugauchi, J., Goguitchaichvili, A., Morales-Contreras, J., 2006. Low-latitude paleosecular variation and the time-averaged field during the late Pliocene and Quaternary—Paleomagnetic study of the Michoacan-Guanajuato volcanic field, Central Mexico. Earth Planets Space 58, 1359–1371.

Cottrell, R.D., Tarduno, J.A., 1999. Geomagnetic paleointensity derived from single plagioclase crystals. Earth Planet. Sci. Lett. 169, 1–5.

Davies, C., Constable, C., 2017. Geomagnetic spikes on the core-mantle boundary. Nat. Commun. 8, 15593.

de Boor, C., 1978. A Practical Guide to Splines. Springer, New York.

de Groot, L.V., Biggin, A.J., Dekkers, M.J., Langereis, C.G., Herrero-Bervera, E., 2013. Rapid regional perturbations to the recent global geomagnetic decay reveled by a new Hawaiian record. Nat. Commun. 4, 2727.

De Marco, E., Spatharas, V., Gómez-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.

Di Chiara, A., Tauxe, L., Speranza, F., 2014. Paleointensity determination from São Miguel (Azores Archipelago) over the last 3ka. Phys. Earth Planet. Int. 234, 1–13.

Dodson, M.H., McClelland-Brown, E., 1980. Magnetic blocking temperatures of single-domain grains during slow cooling. J. Geophys. Res. 85 (B5), 2625–2637.

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

Donadini, F., Kovacheva, M., Kostadinova, M., Hedley, I.G., Pesonen, L.J., 2008. Palaeointensity determination on an early medieval kiln from Switzerland and the effect of cooling rate. Phys. Chem. Earth. 33, 449–457.

Dumberry, M., Finlay, C.C., 2007. Eastward and westward drift of the Earth's magnetic field for the last three millennia. Earth Planet. Sci. Lett. 254, 146–157.

Dunlop, D.J., 2011. Physical basis of the Thellier-Thellier and related paleointensity methods. Phys. Earth planet. Inter. 187, 118–138.

Ertepinar, P., Langereis, C.G., Biggin, A.J., de Groot, L.V., Kulakoğlu, F., Omura, S., Süel, A., 2016. Full vector archaeomagnetic records from Anatolia between 2400 and 1350 BCE: Implications for geomagnetic field models and the dating of fires in antiquity. Earth Planet. Sci. Lett. 434, 171–186.

Fanjat, G., Camps, P., Valdivia, L.A., Sougrati, M.T., Cuevas-Garcia, M., Perrin, M., 2013. First archeointensity determinations on Maya incense burners from Palenque temples, Mexico: New data to constrain the Mesoamerica secular variation curve. Earth Planet. Sci. Lett. 363, 168–180.

Finlay, C.C., Aubert, J., Gillet, N., 2016. Gyre-driven decay of the Earth's magnetic dipole. Nat. Commun. 7, 10422.

Finlay, C.C., 2008. Historical variation of the geomagnetic axial dipole. Phys. Earth Planet. Int. 170, 1–14.

Finlay, C.C., Olsen, N., Toeffner-Clausen, L., 2015. DTU candidate field models for IGRF-12 and the CHAOS-5 geomagnetic field model. Earth Planets Space 67.

Folgheraimer, M., 1899. Sur les variations séculaires de l'inclination magnétique dans l'antiquité. J. Phys. 5, 600–667.

Fox, J.M.W., Aitken, M.J., 1980. Cooling-rate dependence of thermoremanent magnetization. Nature. 283, 462–463.

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

Gallet, Y., Genevey, A., Fluteau, F., 2005. Does Earth's magnetic field secular variation control centennial climate change? Earth Planet. Sci. Lett. 236, 339–347.

Gallet, Y., Genevey, A., Le Goff, M., Warmé, N., Gran-Aymerich, J., Lefèvre, A., 2009. On the use of archeology in geomagnetism, and vice-versa: recent developments in archaeomagnetism. C. R. Phys. 10, 630–648.

Genevey, A., Gallet, Y., 2002. Intensity of the geomagnetic field in western Europe over the past 2000 years: New data from ancient French pottery. J. Geophys. Res. 107 (B11), 2285.

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

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., Jesset, S., Thébault, E., Bouillon, J., Lefèvre, A., Le Goff, M., 2016. New archaeointensity data from French early medieval pottery production (6th–10th century AD). Tracing 1500 years of Geomagnetic field intensity variations in Western Europe. Phys. Earth Planet. Int. 257, 205–219.

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 archaeomagnetic 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), 2858–2872.

Genevey, A., Kondopoulou, D., Pétridis, P., Aidona, E., Muller, A., Blondé, F., Gros, J.S., 2017; al, in press. New constraints on geomagnetic field intensity variations in the Balkans during the Early Byzantine period from ceramics unearthed at Thasos and Delphi. Greece. J. Archaeol Sci.

Gillet, N., Jault, D., Finlay, C.C., 2015. Planetary gyre and time-dependent midlatitude
eddies at the Earth's core surface. J. Geophys. Res. 120, 39914013.

Gillet, N., Jault, D., Finlay, C.C., Olsen, N., 2013. Stochastic modeling of the Earth's magnetic field: inversion for covariances over the observatory era. Geochem. Geophys. Geosyst. 14, 766–786.

Goguitchaichvili, A., Greco, C., Morales, J., 2011. Geomagnetic field intensity behavior in South America between 400 AD and 1800 AD: first archeointensity results from Argentina. Phys. Earth Planet. Int. 186, 191–197.

Goguitchaichvili, A., Morales, J., Schavelzon, D., Vásquez, C., Gogorza, C.S., Loponte, D., Rapalini, A., 2015. Variation of the Earth's magnetic field strength in South America during the last two millennia: new results from historical buildings of Buenos Aires and re-evaluation of regional data. Phys. Earth Planet. Int. 245, 15–25.

Gómez-Paccard, M., Beamud, E., Mc Intosh, G., Larrasoana, J.C., 2013. New archaeomagnetic data recovered from the study of three Roman kilns from north-east Spain: A contribution to the Iberian palaeosecular variation curve. Archaeometry 55, 159–177.

Gómez-Paccard, M., Catanzariti, G., Ruiz-Martínez, V.C., McIntosh, G., Núñez, J.I., Osete, M.L., Chauvin, A., Lanos, Ph., Tarling, D.H., Bernal-Casasola, D., Thiriot, J., 2006a. A catalogue of Spanish archaeomagnetic data. Geophys. J. Int. 166, 1125–1143.

Gómez-Paccard, M., Chauvin, A., Lanos, P., Thiriot, J., 2008. New archeointensity data from Spain and the geomagnetic dipole moment in western Europe over the past 2000 years. J. Geophys. Res. 113 (B9).

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

Gómez-Paccard, M., Chauvin, A., Lanos, P., Thiriot, J., Jiménez-Castillo, P., 2006b. Archeomagnetic study of seven contemporaneous kilns from Murcia (Spain). Phys. Earth Planet. Int. 157, 16–32.

Gómez-Paccard, M., McIntosh, G., Chauvin, A., Beamud, E., Pavón-Carrasco, F.J., Thiriot, J., 2012b. Archaeomagnetic and rock magnetic study of six kilns from North Africa (Tunisia and Morocco). Geophys. J. Int. 189, 169–186.

Gómez-Paccard, M., Osete, M.L., Chauvin, A., Pavón-Carrasco, F.J., Pérez-Asensio, M., Jiménez, P., Lanos, P., 2016. New constraints on the most significant palaeointensity change in Western Europe over the last two millennia. A non-dipolar origin? Earth Planet. Sci. Lett. 454, 55–64.

Gratton, M.N., Goguitchaichvili, A., Conte, G., Shaw, J., Urrutia-Fucugauchi, J., 2005a. Microwave palaeointensity study of the Jorullo volcano (Central Mexico). Geophys. J. Int. 161, 627–634.

Gratton, M.N., Shaw, J., Herrero-Bervera, E., 2005b. An absolute palaeointensity record from SOH1 lava core, Hawaii using the microwave technique. Phys. Earth Planet. Int. 148, 193–214.

Gubbins, D., Jones, A.L., Finlay, C.C., 2006. Fall in Earth's magnetic field is erratic. Science 312, 900–902.

Halgedahl, S.L., Day, R., Fuller, M., 1980. The effect of cooling rate on the intensity of weak-field TRM in single-domain magnetite. J. Geophys. Res. 85 (B7), 3690–3698.

Hammer, R., Harper, D.A.T., Ryan, P.D., 2001. PAST: paleontological statistics software package for education and data analysis. Palaeontol. Electron. 4 9pp.

Hartmann, G.A., Genevey, A., Gallet, Y., Trindade, R.I., Etchevarne, C., Le Goff, M., Afonso, M.C., 2010. Archaeointensity in Northeast Brazil over the past five centuries. Earth Planet. Sci. Lett. 296, 340–352.

Hartmann, G.A., Genevey, A., Gallet, Y., Trindade, R.I., Le Goff, M., Najjar, R., Etchevarne, C., Afonso, M.C., 2011. New historical archeointensity data from Brazil: Evidence for a large regional non-dipole field contribution over the past few centuries. Earth Planet. Sci. Lett. 306, 66–76.

Hartmann, G.A., Trindade, R.I., Goguitchaichvili, A., Etchevarne, C., Morales, J., Afonso, M.C., 2009. First archeointensity results from Portuguese potteries (1550–1750 AD). Earth Planets Space 61, 93–100.

Hill, M.J., Gratton, M.N., Shaw, J., 2002a. A comparison of thermal and microwave palaeomagnetic techniques using lava containing laboratory induced remanence. Geophys. J. Int. 151, 157–163.

Hill, M.J., Gratton, M.N., Shaw, J., 2002b. Palaeomagnetic investigation of Tertiary lava from Barrington Tops, NSW, Australia, using thermal and microwave techniques. Earth Planet. Sci. Lett. 198, 245–256.

Hill, M.J., Shaw, J., 1999. Palaeointensity results for historic lavas from Mt Etna using microwave remagnetization/remagnetization in a modified Thellier-type experiment. Geophys. J. Int. 139, 583–590.

Hulot, G., Le Mouël, J.L., 1994. A statistical approach to the Earth's main dipolar field. Phys. Earth planet. Inter. 82, 162–183.

Hulot, G., Khokhlov, A., Le Mouël, J.L., 1997. Uniqueness of mainly dipolar magnetic fields recovered from directional data. Geophys. J. Int. 129, 347–354.

Jackson, A., 2003. Intense equatorial flux spots on the surface of Earth's core. Nature 424, 760–763.

Jackson, A., Jonkers, A.R., Walker, M.R., 2000. Four centuries of geomagnetic secular variation from historical records. Philos. Trans. R. Soc. Lond. A 358, 957–990.

Kissel, C., Laj, C., Rodriguez-Gonzalez, A., Perez-Torrado, F., Carracedo, J.C., Wandres, C., 2015. Holocene geomagnetic field intensity variations: Contribution from the low latitude Canary Islands site. Earth Planet. Sci. Lett. 430, 178–190.

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.

**Kono, M., 2007. In: Kono, M. (Ed.), Treatise on Geophysics, vol. 5, Elsevier Science, pp. 1–32.**

Korte, M., Constable, C.G., 2011. Improving geomagnetic field reconstructions for 0–3 ka. Phys. Earth Planet. Inter. 188, 247–259.

Korte, M., Donadini, F., Constable, C.G., 2009. Geomagnetic field for 0–3 ka: 2. A new series of time-varying global models. Geochem. Geophys. Geosyst. 10, Q06008.

Krása, D., Heunemann, C., Leonhardt, R., Petersen, N., 2003. Experimental procedure to detect multidomain remanence during Thellier-Thellier experiments. Phys. Chem. Earth. 28, 681–687.

Le Goff, M., Gallet, Y., 2004. A new three-axis vibrating sample magnetometer for continuous high-temperature magnetization measurements: applications to paleo- and archeointensity determinations. Earth Planet. Sci. Lett. 229, 31–43.

Le Goff, M., Gallet, Y., 2017. A reappraisal of instrumental magnetic measurements made in Western Europe before AD 1750: confronting historical geomagnetism and archeomagnetism. Earth Planets Space 69, 32.

Lhuillier, F., Fournier, A., Hulot, G., Aubert, J., 2011. The geomagnetic secular variation timescale in observations and numerical dynamo models. Geophys. Res. Lett. 38, L055306.

Licht, A., Hulot, G., Gallet, Y., Thebault, E., 2013. Ensembles of low degree archaeomagnetic field models for the past three millennia. Phys. Earth planet. Inter. 224, 38–67.

Livermore, P.W., Fournier, A., Gallet, Y., 2014. Core-flow constraints on extreme archeomagnetic intensity changes. Earth Planet. Sci. Lett. 387, 145–156.

Mankinen, E.A., Champion, D.E., 1993. Broad trends in geomagnetic paleointensity on Hawaii during Holocene time. J. Geophys. Res. 98 (B5), 7959–7976.

Michalk, D.M., Muxworthy, A.R., Böhnel, H.N., MacIennan, P., Nowaczyk, N., 2008. Evaluation of the multispecimen parallel differential pTRM method: a test on historical lavas from Iceland and Mexico. Geophys. J. Int. 173, 409–420.

Mitra, R., Tauxe, L., McIntosh, S.K., 2013. Two thousand years of archeointensity from West Africa. Earth Planet. Sci. Lett. 364, 123–133.

Nilsson, A., Holme, R., Korte, M., Suttie, N., Hill, M., 2014. Reconstructing Holocene geomagnetic field variation: new methods, models and implications. Geophys. J. Int. 198, 229–248.

Olson, P., Amit, H., 2006. Changes in Earth's dipole. Naturwissenschaften 93, 519–542.

Osete, M.L., Catanzariti, G., Chauvin, A., Pavón-Carrasco, F.J., Roperch, P., Fernández, V.M., 2015. First archaeomagnetic field intensity data from Ethiopia, Africa (1615 ± 12 AD). Phys. Earth Planet. Inter. 242, 24–35.

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

Paterson, G.A., Tauxe, L., Biggin, Shaar, R., Jonestrask, L.C., 2014. On improving the selection of Thellier-type palaeointensity data. Geochem. Geophys. Geosyst. 15, 1180–1192.

Pavón-Carrasco, F.J., Osete, M.L., Torta, J.M., De Santis, A., 2014. A geomagnetic field model for the Holocene based on archaeomagnetic and lava flow data. Earth Planet. Sci. Lett. 388, 98–109.

Pick, T., Tauxe, L., 1993. Holocene paleointensities: Thellier experiments on submarine basaltic glass from the East Pacific Rise. J. Geophys. Res. 98 (B10), 17949–17964.

Poletti, W., Hartmann, G.A., Hill, M.J., Biggin, A.J., Trindade, R.I.F., 2013. The cooling rate effect on microwave archeointensity estimates. Geophys. Res. Lett. 40, 3847–3852.

Poletti, W., Trindade, R.I.F., Hartmann, G.A., Damiani, N., Rech, R.M., 2016. Archaeomagnetism of Jesuit Missions in South Brazil (1657–1706 AD) and assessment of the South American database. Earth Planet. Sci. Lett. 445, 36–47.

Pressling, N., Brown, M.C., Gratton, M.N., Shaw, J., Gubbins, D., 2007. Microwave palaeointensities from Holocene age Hawaiian lavas: investigation of magnetic properties and comparison with thermal palaeointensities. Phys. Earth planet. Inter. 162, 99–118.

Pressling, N., Laj, C., Kissel, C., Champion, D., Gubbins, D., 2006. Palaeomagnetic intensities from 14 C-dated lava flows on the Big Island, Hawaii: 0–21 kyr. Earth Planet. Sci. Lett. 247, 26–40.

Riisager, P., Riisager, J., 2001. Detecting multidomain magnetic grains in Thellier palaeointensity experiments. Phys. Earth planet. Inter. 125, 111–117.

Rogers, J., Fox, J., Aitken, M.J., 1979. Magnetic anisotropy in ancient pottery. Nature 277, 644–646.

Roperch, P., Chauvin, A., Lara, L.E., Moreno, H., 2015. Secular variation of the Earth's magnetic field and application to paleomagnetic dating of historical lava flows in Chile. Phys. Earth planet. Inter. 242, 65–78.

Salnaia, N., Gallet, Y., Genevey, A., Antipov, I., 2017. New archeointensity data from Novgorod (North-Western Russia) between c. 1100 and 1700 AD. Implications for the European intensity secular variation. Phys. Earth planet. Inter. 269, 18–28.

Sanchez, S., Fournier, A., Aubert, J., Cosme, E., Gallet, Y., 2016. Modelling the archaeomagnetic field under spatial constraints from dynamo simulations: a resolution analysis. Geophys. J. Int. 207, 983–1002.

Schnepf, E., Lanos, P., Chauvin, A., 2009. Geomagnetic palaeointensity between 1300 and 1750 AD derived from a bread oven floor sequence in Lübeck, Germany. Geochem. Geophys. Geosyst. 10 (8).

Shaar, 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.

Shaar, R., Tauxe, L., Ben-Yosef, E., Kassianidou, V., Lorentzen, B., Feinberg, J.M., Levy, T.E., 2015. Decadal-scale variations in geomagnetic field intensity from ancient Cypriot slag mounds. Geochem. Geophys. Geosyst. 16, 195–214.

Shaar, R., Tauxe, L., Goguitchaichvili, A., Devidze, M., Licheli, V., 2017. Further evidence of the Levantine Iron Age geomagnetic anomaly from Georgian pottery. Geophys. Res. Lett. 44, 2229–2236.

Shaar, R., Tauxe, L., Ron, H., Ebert, Y., Finkelstein, I., Agnon, A., 2016. Large geomagnetic field anomalies revealed in Bronze to Iron Age archaeomagnetic data from Tel Megiddo and Tel Hazor. Israel. Earth Planet. Sci. Lett. 442, 173–185.

Shaw, J., Walton, D., Yang, S., Rolph, T.C., Share, J.A., 1996. Microwave archaeointensities from Peruvian ceramics. Geophys. J. Int. 124, 241–244.

Shaw, J., Yang, S., Rolph, T.A., Sun, F.Y., 1999. A comparison of archeointensity results from Chinese ceramics using microwave and conventional Thellier's and Shaw's
methods. Geophys. J. Int. 136, 714–718.

Spassov, S., Valet, J.P., Kondopoulou, D., Zananiri, I., Casas, L., Le Goff, M., 2010. Rock magnetic property and paleointensity determination on historical Santorini lava flows. Geochem. Geophys. Geosyst. 11 (7).

Stanton, T., Riisager, P., Knudsen, M.F., Thordarson, T., 2011. New palaeointensity data from Holocene Icelandic lavas. Phys. Earth planet. Inter. 186, 1–10.

Stark, F., Cassidy, J., Hill, M.J., Shaw, J., Sheppard, P., 2010. Establishing a first ar-
chaeointensity record for the SW Pacific. Earth Planet. Sci. Lett. 298, 113–124.

Stephenson, A., Sadikun, S.T., Potter, D.K., 1986. A theoretical and experimental com-
parison of the anisotropies of magnetic susceptibility and remanence in rocks and
minerals. Geophys. J. Int. 84, 185–200.

Suttie, N., Holme, R., Hill, M.J., Shaw, J., 2011. Consistent treatment of errors in ar-
chaeointensity implies rapid decay of the dipole prior to 1840. Earth Planet. Sci. Lett.
304, 13–21.

Tanaka, H., Kono, M., 1991. Preliminary Results and Reliability of Palaeointensity Studies on Historical and 14C Dated Hawaiian Lvas. J. Geomag. Geoelectr. 43, 375–388.

Tanaka, H., Komuro, N., Turner, G.M., 2009. Palaeosecular variation for 0.1-21 Ka from the Okataina Volcanic Centre. New Zealand. Earth Planets Space 61, 213–225.

Tarduno, J.A., Watkeys, M.K., Huffman, T.N., Cottrell, R.D., Blackman, E.G., Wendt, A., Scribner, C.A., Wagner, C.L., 2015. Antiquity of the South Atlantic Anomaly and evidence for top-down control on the geodynamo. Nat. Commun. 6, 7865.

Tema, E., Gómez-Paccard, M., Kondopoulou, D., Almar, Y., 2012. Intensity of the Earth's magnetic field in Greece during the last five millennia: New data from Greek pottery. Phys. Earth Planet. Inter. 202, 14–26.

Tema, E., Morales, J., Goguitchaichvili, A., Camps, P., 2013. New archaeointensity data from Italy and geomagnetic field intensity variation in the Italian Peninsula. Geophys. J. Int. 193, 603–614.

Terra-Nova, F., Amit, H., Hartmann, G.A., Trindade, R.I.F., 2016. Using archaeomagnetic

field models to constrain the physics of the core: robustness and preferred locations of reversed flux patches. Geophys. J. Int. 206, 1890–1913.

Terra-Nova, F., Amit, H., Hartmann, G.A., Trindade, R.I.F., 2015. The time-dependence of reversed archaeomagnetic flux patches. J. Geophys. Res. 120 (2), 691–704.

Thellier, E., Thellier, O., 1959. Sur l'intensité du champ magnétique terrestre dans le passé historique et géologique. Ann. Geophys. 15, 285–376.

Usoskin, I.G., Gallet, Y., Lopes, F., Kovaltsov, G.A., Hulot, G., 2016. Solar activity during the Holocene: the Hallstatt cycle and its consequence for grand minima and maxima. Astron. Astrophys. 587, A150.

Valet, J.P., Herrero-Bervera, E., LeMouel, J.L., Plenier, G., 2008. Secular variation of the geomagnetic dipole during the past 2000 years. Geochem. Geophys. Geosyst. 9, Q01008.

Veitch, R., Hedley, I., Wagner, J.J., 1984. An investigation of the intensity of the geo-
magnetic field during Roman times sing magnetically anisotropic bricks and tiles.
Arch. Sci. 37, 359–373.

Whaler, K.A., Holme, R., 2011. The axial dipole strength and flow in the outer core. Phys. Earth Planet. Inter. 188, 235–246.

Yoshihara, A., Kondo, A., Ohno, M., Hamano, Y., 2003. Secular variation of the geo-
magnetic field intensity during the past 2000 years in Japan. Earth Planet. Sci. Lett.
210, 219–231.

Yu, Y., Dunlop, D.J., 2003. On partial thermoremanent magnetization tail checks in Thellier paleointensity determination. J. Geophys. Res. 108 (B11).

Yu, Y., 2012. High-fidelity paleointensity determination from historic volcanoes in Japan. J. Geophys. Res. 117 (B8).

Yu, Y., Dunlop, D.J., Özdemir, Ö., 2003. Are ARM and TRM analogs? Thellier analysis of ARM and pseudo-Thellier analysis of TRM. Earth Planet. Sci. Lett. 205, 325–336.

Yu, Y., Tauxe, L., Genevey, A., 2004. Toward an optimal geomagnetic field intensity determination technique. Geochem. Geophys. Geosyst. 5 Q02H07.