# Using archaeomagnetic field models to constrain the physics of the core: robustness and preferred locations of reversed flux patches

Filipe Terra-Nova, Hagay Amit, Gelvam Hartmann, Ricardo I.F. Trindade

► To cite this version:

Filipe Terra-Nova, Hagay Amit, Gelvam Hartmann, Ricardo I.F. Trindade. Using archaeomagnetic field models to constrain the physics of the core: robustness and preferred locations of reversed flux patches. Geophysical Journal International, 2016, 206 (3), pp.1890-1913. 10.1093/gji/ggw248. hal-02344349

HAL Id: hal-02344349

https://hal.science/hal-02344349

Submitted on 10 Sep 2021

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
<figure><img src="image_2.png" /></figure>

# Using archaeomagnetic field models to constrain the physics of the core: robustness and preferred locations of reversed flux patches

Filipe Terra-Nova,<fnref n="1" /> Hagay Amit,<fnref n="1" /> Gelvam A. Hartmann<fnref n="2" />,<fnref n="3" /> and Ricardo I.F. Trindade<fnref n="4" />

Accepted 2016 June 27. Received 2016 June 27; in original form 2016 January 21

## SUMMARY

Archaeomagnetic field models cover longer timescales than historical models and may therefore resolve the motion of geomagnetic features on the core–mantle boundary (CMB) in a more meaningful statistical sense. Here we perform a detailed appraisal of archaeomagnetic field models to infer some aspects of the physics of the outer core. We characterize and compare the identification and tracking of reversed flux patches (RFPs) in order to assess the RFPs robustness. We find similar behaviour within a family of models but differences among different families, suggesting that modelling strategy is more influential than data set. Similarities involve recurrent positions of RFPs, but no preferred direction of motion is found. The tracking of normal flux patches shows similar qualitative behaviour confirming that RFPs identification and tracking is not strongly biased by their relative weakness. We also compare the tracking of RFPs with that of the historical field model gufm1 and with seismic anomalies of the lowermost mantle to explore the possibility that RFPs have preferred locations prescribed by lower mantle lateral heterogeneity. The archaeomagnetic field model that most resembles the historical field is interpreted in terms of core dynamics and core–mantle thermal interactions. This model exhibits correlation between RFPs and low seismic shear velocity in co-latitude and a shift in longitude. These results shed light on core processes, in particular we infer toroidal field lines with azimuthal orientation below the CMB and large fluid upwelling structures with a width of about 80° (Africa) and 110° (Pacific) at the top of the core. Finally, similar preferred locations of RFPs in the past 9 and 3 kyr of the same archaeomagnetic field model suggest that a 3 kyr period is sufficiently long to reliably detect mantle control on core dynamics. This allows estimating an upper bound of 220–310 km for the magnetic boundary layer thickness below the CMB.

**Key words:** Archaeomagnetism; Dynamo: theories and simulations; Palaeomagnetic secular variation.

## 1 INTRODUCTION

The movement of an electrically conductive fluid in the Earth's outer core generates the geomagnetic field and controls its secular variation (SV). The SV represents the changes of the Earth's magnetic field over various timescales. Observation-based timescales range from yearly as geomagnetic jerks (Courtillot et al. 1978; Mandea et al. 2010), decadal to centennial as drift of geomagnetic flux patches (Holme 2015) and dipole moment decrease (Olson & Amit 2006), to millennial as the duration of a reversal or an excursion of the geomagnetic dipole (Merrill et al. 1998). Longer timescales are also observed, such as the average duration of a chron, the average

duration of a superchron and reversal frequency irregular variability, though these very long timescales are probably controlled by the mantle (rather than the core) via core–mantle coupling (e.g. Olson et al. 2013).

Timescales of core dynamics can also be assessed theoretically. The magnetic induction equation relates the flow in the Earth's core to changes in the magnetic field. The SV ($\partial \vec{B}/\partial t$) is given by:

$$ \frac{\partial \vec{B}}{\partial t} = \nabla \times (\vec{u} \times \vec{B}) + \eta \nabla^2 \vec{B} \quad (1) $$

where $\vec{u}$ is the velocity, $\vec{B}$ is the magnetic field and $\eta$ is the magnetic diffusivity. The SV is the outcome of electromagnetic induction

<footnote n="1">¹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. E-mail: filipe.terranova@univ-nantes.fr</footnote>

<footnote n="2">²Instituto de Geociências, Universidade Estadual de Campinas, Rua João Pandiá Calógeras, 51, 13083-870 Campinas, Brasil</footnote>

<footnote n="3">³Observatório Nacional, Rua General José Cristino, 77, 20921-400 Rio de Janeiro, Brasil</footnote>

<footnote n="4">⁴Departamento de Geofísica, Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo, Rua do Matão, 1226, Cidade Universitária, 05508-090 São Paulo, Brasil</footnote>
$\nabla \times (\vec{u} \times \vec{B})$ (field generation) and ohmic diffusion $\eta \nabla^2 \vec{B}$ (field destruction). The relative importance of these effects is quantified by the magnetic Reynolds number, which is the ratio of the magnetic diffusion time $\tau_{\eta} = L^2/\eta$ and the fluid advection time $\tau_{u} = L/U$,

$$R_m = \frac{\tau_{\eta}}{\tau_{u}} = \frac{UL}{\eta}, \qquad (2)$$

where $L$ and $U$ are typical length and velocity scales respectively. Using $L \sim 1000$ km, $U \sim 5 \times 10^{-4}$ ms$^{-1}$ and $\eta \sim 1$ m$^2$s$^{-1}$ the timescales are $\tau_{\eta} \sim 30$ kyr and $\tau_{u} \sim 60$ yr, and $R_m \sim 500$ (Bloxham & Jackson 1991; Holme 2015). Accordingly, over short timescales (decadal to centennial), the SV is expected to be strongly dominated by advection.

To properly investigate field changes using observational data two ingredients are needed. First, a sufficiently high resolution model is required to resolve the field morphology. Second, a long period is essential to cover as many timescales as possible. Earth's magnetic field models constructed with satellites data (e.g. Olsen et al. 2014; Finlay et al. 2015) resolve well the smaller lengthscales, but only cover a very short period. The historical geomagnetic field model gufm1 (Jackson et al. 2000) covers some core timescales but not all, in particular the gufm1 period is of the same order as the advection time. Archaeomagnetic field models typically extend to longer time periods and cover more core timescales, therefore resolving the motion of geomagnetic features on the core-mantle boundary (CMB) in a more meaningful statistical sense. However, these models have lower resolution due to lack of global coverage and high uncertainties of their data sets (e.g. Donadini et al. 2009; Brown et al. 2015a,b), therefore their robustness must be assessed before drawing any geophysical implications.

Intense flux patches on the CMB are considered main features of the Earth's magnetic field (Christensen et al. 2010). Normal flux patches (NFPs) are observed near the tangent cylinder (Jackson et al. 2000; Olsen et al. 2014) possibly due to fluid downwelling concentrating field lines there (Olson et al. 1999). NFPs also appear at low- and mid-latitudes (Jackson 2003). Their westward drift is one of the most prominent features of the geomagnetic field (Finlay & Jackson 2003; Aubert et al. 2013). Regions bounded by a null-flux curve with a reversed radial field with respect to its hemisphere are the so-called reversed flux patches (RFPs). The most intense RFPs over the past decades are found below the South Atlantic (e.g. Jackson et al. 2000; Olsen et al. 2014).

It has been proposed that the dynamical origin of RFPs is toroidal field expulsion by fluid upwelling. Gubbins & Bloxham (1986) considered a simple 2-D kinematic model of convection to illustrate that flux expulsion, through the action of radial diffusion on a strong toroidal field advected towards the CMB by upwelling, is a viable mechanism to explain RFPs. This process introduces a much smaller radial lengthscale due to the upwelling, hence diffusion effects could locally become comparable to advection (Amit & Christensen 2008). Numerical dynamo studies confirm that the RFPs emerge by the expulsion of toroidal field due to fluid upwelling, a process which may eventually trigger reversals (Wicht & Olson 2004; Takahashi et al. 2007; Aubert et al. 2008; Olson et al. 2009). It is thought that the proliferation of RFPs on the CMB and their poleward motion as well as equatorward advection of NFPs cause the rapid decrease in the dipole intensity over the historical period (Gubbins 1987; Gubbins et al. 2006; Olson & Amit 2006; Finlay et al. 2012, 2016). In contrast, frozen-flux core field models may also be constructed (O'Brien et al. 1997; Jackson et al. 2007). Therefore, tracking RFPs and understanding the physical

mechanism responsible for their temporal evolution is fundamental to better understand the dynamics of the geomagnetic field.

Analyses of archaeomagnetic field models mostly focused on the kinematics of high-latitudes intense normal polarity flux patches. These patches were found to be mobile with alternating eastward-westward drifts (Dumberry & Finlay 2007; Wardinski & Korte 2008; Korte & Holme 2010; Amit et al. 2011). Recently, Terra-Nova et al. (2015) designed topological algorithms to identify and track RFPs in time and applied it to the archaeomagnetic field model of Korte & Constable (2011). Terra-Nova et al. (2015) used the magnetic equator to define reversed flux regions instead of the geographic equator thus eliminating ambiguity in RFP identification at the equator. They found that most RFPs exhibit westward drift and migrate towards higher latitudes. They also suggested that the power in spherical harmonic degrees 4 and above strongly affect the existence of RFPs. Based on low-pass filter sensitivity tests Terra-Nova et al. (2015) suggested that RFPs are robust features of the archaeomagnetic field model.

The main goal of this paper is to characterize and compare time-dependent RFPs in various archaeomagnetic field models in order to assess their robustness. We examine the similarity within a family of archaeomagnetic field models as well as among different families of models. We compare the tracking results with that obtained from the historical field model gufm1 (Jackson et al. 2000) and with seismic anomalies of the lowermost mantle (Masters et al. 2000) to explore the possibility that RFPs have preferred locations prescribed by lower mantle heterogeneity.

We analyse the archaeomagnetic field models constructed by Korte & Constable (2011), Licht et al. (2013) and Nilsson et al. (2014) for the past three millennia to explore the dependence of RFPs temporal evolution on data source, data treatment and modelling strategy including the way these models deal with outliers. From hereafter these papers are denoted as KC11, LHGT13 and NHKSH14, respectively. Section 2 summarizes differences and similarities among these archaeomagnetic field models. In section 3 we briefly summarize the topological algorithms introduced by Terra-Nova et al. (2015) to define, identify and track RFPs. In order to assess the reliability of the RFPs identification and tracking results we also apply these algorithms to NFPs. In section 4 we present the results of the tracking of RFPs and NFPs. We further explore the time-dependence of RFPs for a wider time window with the pfm9k models (NHKSH14) that cover 9 kyr. Finally, in Section 5 we discuss the robustness of RFPs and some possible geophysical implications, in particular the mantle control on the locations of RFPs.

## 2 ARCHAOMAGNETIC FIELD MODELS

Here we compare the archaeomagnetic field models of KC11, LHGT13 and NHKSH14. We focus on data source and treatment as well as modelling strategy including dealing with outliers. KC11 introduced two models: CALS3k.4 and CALS3k.4b. Here we only consider CALS3k.4b, since it is a more conservative model which maintains only the most robust spatial and temporal features, thus it is more suitable for the representation of the field at the CMB. LHGT13 built three families of models based on three types of data sets: A_FM, ASD_FM and ASDI_FM. NHKSH14 presented another family of field models: pfm9k.1, pfm9k.1a, pfm9k.1b. All these models or ensembles of models are discussed below. The models' characteristics are summarized in Tables S1 and S2 (Supporting Information).
## 2.1 Data set and treatment

There are differences in the data sources among all models. The CALS3k.4b field model was constructed from two data sources. First, from the GEOMAGIA V.2 data base (http://geomagia.ucsd.edu/) of archaeomagnetic and volcanic data sets updated until August 2009. Second, from the sedimentary data set consisting of SED3K_dat0 (Donadini et al. 2009) plus another compilation of data from 13 new locations (KC11). LHGT13 used the same sedimentary data set as in CALS3k.4b, preserving preliminary data corrections and adjustments done by KC11. For archaeomagnetic and volcanic data set they also used the GEOMAGIA V.2 but updated until July 2011. NHKSH14 obtained archaeomagnetic data from the online GEOMAGIA50 data base (Donadini et al. 2006; Korhonen et al. 2008) updated until August 2013 and sedimentary data from the SED12k data compilation (Donadini et al. 2009; Korte et al. 2011).

The types of data were slightly varied among all models analysed. CALS3k.4b uses archaeomagnetic, volcanic and sedimentary data, while each family of models from LHGT13 uses different types of data source. A_FM uses only archaeomagnetic and volcanic data, ASD_FM uses the latter data set plus sedimentary direction data, and ASDI_FM uses the entire data set. All pfm9k models use archaeomagnetic and sedimentary data (NHKSH14). A summary of the data sources used by the various archaeomagnetic field models is given in Supporting Information Table S1.

Sedimentary intensity data are always relative, while declination may be oriented, although usually that is not the case. KC11 nevertheless considered the sedimentary data as being oriented and calibrated relative sedimentary intensity data by comparing their values with predictions from the CALS3k.3 model (Korte et al. 2009) to construct a first guess model. LHGT13 calibrated the relative intensity sedimentary data using the archaeomagnetic based A_FM-M mean model. NHKSH14 resampled the sedimentary data to avoid inappropriate weighting of the data. They binned all sedimentary data records in 50 yr intervals, giving equal weight to each site at any given time. This led to more than 70 per cent reduction of the sedimentary data, adding weight to the archaeomagnetic data (NHKSH14). Uncertainties in intensity were obtained by calibrating the model against the dipole model. The model used information from ¹⁰Be flux data from grip ice core (Muscheler et al. 2004; Vonmoos et al. 2006) and from the historical geomagnetic field model gufm1 (Jackson et al. 2000). Missing information from 1350 AD until 1590 AD (although in practice gufm1 is constrained by intensity measurements only from 1840 AD) was filled by linear interpolation. Then, NHKSH14 calibrated each sedimentary data declination by a constant number of degrees, akin to the preliminary calibration of relative palaeointensity, where each entire record was multiplied by a constant factor to obtain absolute intensities. These constants were based on comparisons with the prior dipole field model or with archaeomagnetic data, whichever better fitted the data.

Error assessment may play an important role in the construction of field models (Suttie et al. 2011). To assign error estimates, KC11 used the same scheme applied for the CALS3k.3 model. Based on the average deviation between archaeomagnetic data and the gufm1 model during the historical era, the confidence circle of direction (Fisher 1953) for archaeomagnetic data ($\alpha_{95a}$) was assigned a minimum value (Supporting Information Table S1). Since the statistics are less reliable for the sedimentary data because of the lack of such recent records, the minimum value of $\alpha_{95a}$ that was set by KC11 for sediments is higher (Supporting Information Table S1). It is worth

noting that very few sediment records come with magnetic uncertainties and thus nearly all these records are weighted equally. For intensity data, the minimum uncertainty $\sigma_f$ (the standard deviation of the intensity measurements) was the same for both data types (Supporting Information Table S1).

Special attention was given by LHGT13 to data uncertainties conserving previous published errors. They argued that even high quality data should not be fitted too closely by the models due to their intrinsically limited resolution. The approach of LHGT13 is therefore different from that of Donadini et al. (2009), in order to better identify high quality data. Their approach concerning data with unknown measurements errors is similar to that of Donadini et al. (2009), although LHGT13 used RMS values of the published errors of each data type multiplied by 1.5 to slightly penalize data with unknown errors (Supporting Information Table S1). Models of LHGT13 were truncated at spherical harmonic degree and order 5, which led them to introduce a modelling error due to the unknown contribution from higher degrees, in addition to the independent measurement error. Their assigned error estimates are a combination of both types of errors. The gufm1 historical model (Jackson et al. 2000) up to degree 5 was used to derive the modelling error, with the results of $\sigma_f^{m} = 2 \mu\text{T}$ and $\alpha_{95a}^{m} \sim 3.5^{\circ}$ being quite conservative (LHGT13). Nevertheless, the modelling errors chosen by LHGT13 leave the possibility for high quality data to be assigned with a total error (measurement error combined with modelling error) less than the minimum error assigned to all data by Donadini et al. (2009). Since measurements and modelling errors are independent, the total errors were added quadratically.

The error estimate of NHKSH14 is slightly different from the scheme of Donadini et al. (2009). To penalize archaeomagnetic data with less known uncertainties, different minimum errors were assigned depending on the number of samples/specimens (N/n) used to calculate the direction or intensity. NHKSH14 assigned errors as follows. The $\alpha_{95a}$ parameter depends on N, and if data have unknown uncertainties, the higher the N the lower the error (Supporting Information Table S1). For archaeomagnetic intensity, the standard error of the mean $s_f$ depends on n; $s_f$ is given in per cent of true field, which is calculated with the prior dipole model. The higher the n the lower is $s_f$ (Supporting Information Table S1). For two sedimentary records of directions with good error estimates $\alpha_{95a}$ is low. For all other records, a fixed $s_f$ was attributed (Supporting Information Table S1). Lastly, for sedimentary intensity $s_f$ is fixed for all data (Supporting Information Table S1). The assigned uncertainties provided errors to sedimentary data that are on average larger than in Donadini et al. (2009), with some data with slightly smaller error estimations (NHKSH14).

No minimum value for age uncertainty was used in CALS3k.4b since most archaeomagnetic artefacts and lavas from historical time can be dated relatively accurately. For archaeomagnetic data with no age uncertainties KC11 fixed $\sigma_a$ (Supporting Information Table S1), where $\sigma_a$ is the standard deviation of age uncertainties in archaeomagnetic data. Sediment age uncertainties were all fixed to $\sigma_s$ (Supporting Information Table S1), where $\sigma_s$ is the standard deviation of sediment age uncertainties, which was applied randomly by the bootstrap method described below. In the families of models A_FM, ASD_FM and ASDI_FM age uncertainties of archaeomagnetic and volcanic data were assigned with an approach similar to that used for direction and intensity. LHGT13 used bins of 500 yr, multiplying the mean average uncertainty from all archaeomagnetic and lava data by 1.5 (Supporting Information Table S1). NHKSH14 addressed age uncertainties of archaeomagnetic data by a minimum value of $\sigma_a$ according to the age of the sample; the older the record
the higher is σₐ (Supporting Information Table S1). This avoids overestimating age uncertainties for recent data (historical data). NHKSH14 made a special effort to treat age uncertainties of sedimentary data, as each of their models used a different approach. The pfm9k.1 model is free of age errors in sedimentary records. The pfm9k.1a model uses an alternative way of adjusting the sediment data chronologies by randomly stretching and compressing the individual timescales of each record in order to assess the likely range of age uncertainties and adjust the timescale of each record based on comparisons with predictions from pfm9k.1. The pfm9k.1b addresses all data uncertainties using the MAST bootstrap methodology (detailed in Korte et al. 2009) using the archaeomagnetic data normal distribution centred in the age estimate and their respective standard errors. Yet for the temporal sampling the timescale of sedimentary data was stretched and compressed using the same routine used in pfm9k.1a (Supporting Information Table S1).

## 2.2 Modelling strategy

For spatial regularization all models minimized the Ohmic dissipation at the CMB (Gubbins 1975). The time-dependence of the Gauss coefficients is expanded by cubic B-splines. The models vectors contain the Gauss coefficients to be recovered. The temporal regularization minimizes the second time derivative of the radial field at the CMB (Bloxham & Jackson 1992). Models are obtained by minimizing a cost function that combines data misfit and spatial and temporal complexities. All models follow Korte et al. (2009) and exclude the dipole field from the spatial regularization.

Differences among the models are related to some parameter choices, such as spherical harmonic truncation, knot spacing and damping coefficients of spatial and temporal complexities. The models of KC11 and NHKSH14 were expanded until the same spherical harmonic degree and order of 10 but with different knot spacings, while the models of LHGT13 were expanded until a lower degree of 5 with a different knot spacing (Supporting Information Table S2). CALS3k.4b damping parameters were chosen by visual comparison of power spectra of the main field and SV to those of the gufm1 (Jackson et al. 2000) and the International Geomagnetic Reference Field (IGRF) 10th generation for 2005 (Maus et al. 2005) as in CALS3k.3 (Korte et al. 2009; Supporting Information Table S2). LHGT13 based their damping parameters on the functional form of the core field spectrum proposed by Hongre et al. (1998) and applied to degrees 2–7 of the historical field spectrum of Bloxham & Jackson (1992). They took advantage of the different data sets to construct reference models (A_FM-0000, ASD_FM-0000 and ASDI_FM-0000) based on pairs of damping parameters that best comply with the spectral expectations. The values obtained were different for each family of models (Supporting Information Table S2). The family of models A_FM had the highest value of λ but the lowest of τ. The families of models ASD_FM and ASDI_FM had the same value of λ with the former having the highest value of τ (Supporting Information Table S2). All LHGT13 models had both λ and τ higher than the damping coefficients of CALS3k.4b (Supporting Information Table S2). NHKSH14 chose the damping parameters by visual comparison (Lodge & Holme 2009) of the time-averaged geomagnetic power spectra of the main field and SV with those of gufm1 (Jackson et al. 2000) and CALS3k.4 (KC11).

Next we compare modelling strategies including the use of bootstraps and the way each study dealt with outliers (Supporting Information Table S2). To provide an empirical measure of the errors affecting archaeomagnetic field models and their predictions, the

standard approach consists in relying on the concept of an empirical covariance matrix. Normally the covariance matrix only addresses measurements errors (Backus 1988a,b). To avoid this pitfall, Korte et al. (2009) introduced the MAST (magnetic, age, spatial and temporal) bootstrap method. CALS3k.4b is produced by averaging the 2000 individual MAST bootstrap models using the final data set of CALS3k.4 (initial data set after three iterations of rejection of outliers and recalibration of intensity sedimentary data). The errors used in the bootstrap method of CALS3k.4b are classified into four distributions. A normal distribution is centred on each data type: archaeomagnetic data, sedimentary data, and archaeological age data, using the assigned uncertainties described above. For the sedimentary data a uniform distribution within ±300 yr of the original ages is used (KC11).

The approach above provides only estimates of variances of the Gauss coefficients in time, whereas the ST bootstrap procedure is useful to assess how sensitive models are to possible removal of information provided by some sites. LHGT13 used bootstrap methods in almost all models. The reference models were built with the final data set and outlier reweighting (further explained in this section), but did not use bootstrap methods. The ensembles of models were built using the MA (magnetic age) bootstrap method (Korte et al. 2009; Thébault & Gallet 2010) to account for magnetic data and timing errors (LHGT13). To produce an ensemble of 1000 plausible data sets for each of the three original data sets, LHGT13 added Gaussian random error to the observed magnetic value independently of the type of data, using the standard deviations detailed above. For archaeomagnetic and volcanic age data they used a random value in a uniform distribution with the age bracket defined by the age uncertainties described above. For sedimentary age data with known uncertainties a Gaussian random error is added to that provided with the data with standard deviation typically ranging from 40 to 100 yr (LHGT13). The full ensemble of models provides relevant statistics for quantities that are linear or non-linear to the Gauss coefficients as well as possible correlations in both space and time of Gauss coefficients. The models of LHGT13 do not discard outliers, but re-weight them. Data classified as outliers are given less weight by increasing their uncertainties. The models A_FM-M, ASD_FM-M and ASDI_FM-M are the mean models of the 1000 MA bootstrap models for each ensemble, while A_FM-XXXX, ASD_FM-XXXX and ASDI_FM-XXXX (XXXX from 0001 to 1000) indicate a certain sample of the MA bootstrap method of each ensemble.

Bootstrap is not used in pfm9k.1 and pfm9k.1a models (NHKSH14). These models are built by three iterations of rejection of outliers and recalibration of sedimentary intensity and declination data. The models differ in addressing sedimentary age uncertainties (see above). The pfm9k.1b model investigates uncertainties using the MAST bootstrap method. Initially, it also uses the above methodology of rejection/calibration but with a relaxed temporal damping (NHKSH14), chosen by visual comparison to the SV power spectrum of gufm1 (Jackson et al. 2000). NHKSH14 drew the individual data from the bootstrap method as KC11, using normal distributions centred on each data type (archaeomagnetic data, sedimentary data, and archaeological age data), but using their uncertainties (see above). However, for the temporal sampling of sedimentary data they applied an alternative method, by randomly stretching and compressing each data timescale (described above and detailed in NHKSH14). On this data set, 2000 MAST bootstrap models were averaged to obtain the pfm9k.1b model.

Hellio et al. (2014) showed that in order to properly account for dating uncertainties only the most probable set of dates should be
considered. Indeed, randomly sampling all possible dates of mea-
surements results in a large amount of unprobable models. Conse-
quently, the obtained models might be severely smoothed.

It is worth noting some modelling particularities. The
CALS3k.4b model is strongly constrained by gufm1 during histori-
cal times. This constraint causes a discontinuity that is smoothed out
around 16th and 17th century AD (KC11). The ensembles of MA
bootstrap models (A_FM, ASD_FM and ASDI_FM) allow for sig-
nificant time-varying non-diagonal terms in the covariance matrix
of covariance between pairs of time-dependent Gauss coefficients
(LHGT13). In an attempt to minimize the underestimation of in-
tensity by regularized models, NHKSH14 resampled the data and
increased the weight given to all intensity data by 50 per cent by
reducing the uncertainty estimates of the data.

Finally, we note that the data record probably contains unre-
alistic variations in the Australian/South-East Asian region due
to some declination records that remained unoriented and hence
systematically offset (Panovska et al. 2015). Thus we also con-
sider the CALS3k.3 field model, which does not contain these
unoriented data, to test the RFPs robustness. Perhaps more im-
portantly, the gufm1 constraint is much weaker in CALS3k.3
than in CALS3k.4b (Korte et al. 2009, KC11), hence RFPs in
CALS3k.3 may be equally distributed in time as opposed to the
concentration of RFPs during the historical period in CALS3k.4b
(Terra-Nova et al. 2015).

3 METHOD

**3.1 Identification and tracking**

To identify RFPs, Terra-Nova et al. (2015) introduced topological
algorithms that consist of four steps at each snapshot: (i) Identifi-
cation of the magnetic equator; (ii) attribution of each grid point to
one magnetic hemisphere; (iii) identification of peaks as centres of
RFPs; and lastly (iv) imposition of an intensity criterion to filter out
weak RFPs. Using the magnetic equator to define the local polarity
avoids interpreting undulation of the magnetic equator as RFPs and
equatorial patches as partly normal partly reversed. To track RFPs
in time, Terra-Nova et al. (2015) coded an algorithm that calculates
the spherical distance along a great circle of each RFP to all RFPs
in the next snapshot. Then they considered a critical velocity of
70 km yr⁻¹ (based on the maximum of time-average core flow plus
1 standard deviation from the geomagnetic SV inversions of Amit
& Olson (2006)). For each model this critical velocity is multiplied
by the time step between snapshots to obtain the critical distance
for tracking. Two RFPs at two successive snapshots that have a
spherical distance lower than the critical value are denoted as the
same RFP. In the new snapshot if an RFP is farther than the critical
distance from all RFPs of the previous snapshot, it is denoted as a
new RFP. For more details see Terra-Nova et al. (2015).

We apply the same algorithms of Terra-Nova et al. (2015) to iden-
tify and track RFPs, including the same magnetic intensity criterion
that filters all RFPs with peaks lower than half the most intense
RFP, as well as the same 70 km yr⁻¹ velocity threshold for tracking.
In addition, in order to assess the reliability of the identification
and tracking of RFPs, we adjust and apply the algorithms described
above to NFPs. Because NFPs are more intense than RFPs and thus
more abundant, a stronger intensity criteria is applied. Only NFPs
with peak values larger than three quarters the most intense NFP of
the same snapshot are considered. In Fig. 1, we show results of the
identification of RFPs (red X symbols) and NFPs (green diamonds)

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The radial geomagnetic field at the core–mantle boundary at (a) 0160 BC, (b) 0920 AD and (c) 1520 AD in the A_FM-M model. The dashed black line is the mapped magnetic equator, the red crosses are the identified reversed flux patches and the green diamonds are the identified normal flux patches.</figcaption></figure>

for three snapshots from the A_FM-M model. The algorithms per-
form well for both RFPs and NFPs.

**3.2 Choosing samples from ensembles of models**

Each mean model of LHGT13 (A_FM-M, ASD_FM-M and ASDI_FM-M) was built from 1000 samples of the bootstrap method, with the full ensemble of models also provided (A_FM-XXXX, ASD_FM-XXXX and ASDI_FM-XXXX, where XXXX indicates sample number). LHGT13 argued that each model of these ensembles is an adequate fit to the data within the error. From each ensemble of models we sample three models which contain distinctive power spectrum behaviour. This selection is relevant because strongly damped models are less likely to accommodate many patches (in particular RFPs) whereas weakly damped models are more likely to contain such smaller scale features. Indeed, Terra-Nova et al. (2015) demonstrated that the level of filtering applied to field models may affect the identified RFPs.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) The power spectrum $R_n$ in µT² versus spherical harmonic degree n for all samples of the A_FM-M ensemble of models in yellow curves. Three coloured curves represent different types of damping. (b) Slopes of linear regression *a* for the curves in (a) excluding n = 1 versus sample number. (c) As in (a), only for the chosen models based on (b), and including standard deviations in time. Inset shows zoom into spherical harmonic degrees 3 and 4.</figcaption></figure>

Nevertheless other selection processes of samples from an ensemble, for example Empirical Orthogonal Function/Principal Component Analysis (e.g. Pais *et al.* 2015), are worth exploring.

First we calculate for each model the time average power spectrum for the period 1000 BC until 2000 AD (yellow curves in Fig. 2(a) for the A_FM ensemble). Then, we fit the non-dipole part of the spectra with linear regressions, and plot the slopes for each sample (Fig. 2b). Next, to avoid too extreme behaviours, we factor the maximum value of all slopes by 0.99 and the minimum by 1.01, and we choose the samples that are the closest to these values. We also calculate the mean of all slopes and choose the sample that has the closest value to it.

In Fig. 2(a), we show the chosen samples, green for the closest to the maximum times 0.99 representing weak damping, blue for the

closest to the mean (moderate damping) and black for the closest to the minimum times 1.01 (strong damping). Lastly, we draw again the time averaged power spectra but only for the three chosen models. Here we also add the standard deviation representing the time-dependence of the spectra (see red error bars in Fig. 2c, with a zoom into spherical harmonics 3 and 4 for better visualization). Weak temporal variability is observed, that is, the time averaged spectra are representative of individual snapshots. Moreover, the chosen models showed different spectrum behaviour among them, therefore each chosen model may indeed represent a different level of spatial damping. The above methodology is applied for each ensemble of LHGT13 models. The chosen models are listed in Supporting Information Table S2.

# 4 RESULTS

## 4.1 Reversed flux patches tracking

Figs 3 and 4 show the tracking of RFPs in longitude ($\phi$) and co-latitude ($\theta$) for the ASDI_FM family of models of LHGT13. The similarity among ASDI_FM models is high. For example, two robust RFPs are resolved in all models, around 500 AD at 130°E and around 1600 AD at longitude 0°, with little longitudinal drift for both RFPs in all models. The number of RFPs that migrate equatorward and poleward are similar in all models. We note a high recurrence of RFPs along co-latitude 60° before 500 BC. Robust RFPs are found around 500 AD at co-latitude 110° and around 1600 AD at co-latitude 150°, exhibiting little net north-south motion. For results of the families of models A_FM and ASD_FM see Figs S1–S4 and their description in the Supporting Information section.

Overall within each family of LHGT13 models tracking results have similarities. An exception is the A_FM-0677 model that shows a weak similarity with the other models of its family. The three mean models (A_FM-M, ASD_FM-M, ASD_FM-M), as expected, usually show RFPs that already appear in the other models of the family, and just in rare cases the M models contain RFPs that are not resolved in any other model of its family.

Figs 5 and 6 show the tracking for pfm9k.1, pfm9k.1a and pfm9k.1b models of NHKSH14, for CALS3k.3 of Korte *et al.* (2009) and the results of Terra-Nova *et al.* (2015) for CALS3k.4b model of KC11. The tracking in longitude of pfm9k.1, pfm9k.1a and pfm9k.1b (Figs 5a–c) show most RFPs with westward drift before 1000 AD and with eastward drift after. For example, a very long-lived RFP is seen drifting westward in all three models around longitude 150°E before 0 AD. The westward drift of the RFPs in the pfm9k models is much weaker than that obtained for RFPs in the historical field model gufm1 (Terra-Nova *et al.* 2015), probably indicating that the reliability of the RFPs drift at ancient times is low. The similarity of tracking among the three models is very high, both in position and in drift direction. At least six RFPs are clearly resolved in all three models, five of them after 320 AD. Most of these RFPs appear near longitude 0°, and they do not have a preferred drift direction. In Figs 6(a)–(c), a similar number of RFPs move equatorward and poleward for all pfm9k models. The long-lived RFP between 1000 BC and 0 AD is located near co-latitude 60° in all pfm9k models. The similarity of the longitudinal tracking among the pfm9k models is therefore also high. It is worth noting the scarcity of RFPs resolved in the Southern Hemisphere. The longitudinal tracking of the CALS3k.4b model (Fig. 5d) shows two RFPs with westward drift and two quasi-stationary RFPs
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Time-dependent longitude in degrees of reversed flux patches for the family of models ASDI_FM from 1000 BC to 2000 AD. The same coloured diamonds are used for the position of a given RFP. The ASDI_FM-0962, ASDI_FM-0425 and ASDI_FM-0188 models represent minimum, mean and maximum damping, respectively.</figcaption></figure>

before 550 AD, while for recent times (after 1450 AD) there is a higher occurrence of RFPs with dominantly westward drift (Terra-Nova *et al.* 2015). In co-latitude most of the CALS3k.4b RFPs migrate poleward (see Fig. 6d; Terra-Nova *et al.* 2015). In lon-

gitude CALS3k.3 shows more westward drifting RFPs, especially for recent times (Fig. 5e). Its longest-lived RFP has been drifting eastward before about 700 BC. In co-latitude, CALS3k.3 shows significantly more Northern Hemisphere RFPs with little preference
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Time-dependent co-latitude in degrees of reversed flux patches for the family of models ASDI_FM from 1000 BC to 2000 AD. The same coloured diamonds are used for the position of a given RFP. The ASDI_FM-0962, ASDI_FM-0425 and ASDI_FM-0188 models represent minimum, mean and maximum damping, respectively.</figcaption></figure>

of drift direction (Fig. 6e). Note that CALS3k.3 has roughly equal distribution of RFPs over time (Figs 5e and 6e), as opposed to the concentration of RFPs in the historical period in CALS3k.4b (Terra-Nova *et al.* 2015).

Next we compare the tracking of all families of models. It is clear that the similarity between models is lower among different families of models than within the same family. Nevertheless, there are robust RFPs that are common to most models. For example, a
<figure><img src="image_7.png" /><figcaption>**Figure 5.** As in Fig. 3 for the families of models pfm9k and CALS3k. For the CALS3k models the temporal sampling is higher (points are denser).</figcaption></figure>

long-lived RFP is seen in all models, except in CALS3k.4b, between longitude 120°E and 150°E and near co-latitude 60°. It shows no preferred drift motion. In the ASDI_FM family of models it appears more as a region of proliferation of RFPs, which may be a result of the intensity criterion (Terra-Nova *et al.* 2015). In the A_FM family

this RFP has a short lifetime. We find a high similarity between the ASD_FM and ASDI_FM families of models with both models showing RFPs mainly in the Southern Hemisphere. For example, an RFP that is resolved after 1500 AD is seen for all ten analysed models of these two families. In addition, the pfm9k family of
<figure><img src="image_8.png" /><figcaption>**Figure 6.** As in Fig. 4 for the families of models pfm9k and CALS3k.</figcaption></figure>

models is the only one that shows similarity to CALS3k.4b, and only for the northern RFPs after 1500 AD since pfm9k does not show RFPs in the Southern Hemisphere for this time.

Table 1 quantifies the statistical behaviour of the tracked RFPs. To distinguish drifting RFPs from quasi-stationary RFPs we use one quarter of the weighted average rate of all RFPs at each model

as a critical value, with the weight corresponding to the lifetime of each RFP. The A_FM models (LHGT13) present slightly more RFPs with westward drift, although the eastward rate is either larger than or comparable to the westward rate (Table 1). Only the A_FM-0677 and A_FM-M models resolved RFPs with quasi-stationary behaviour. The horizontal vector velocity is highest for
<table><caption>**Table 1.** Statistics of azimuthal motion of RFPs.</caption><thead><tr><th rowspan="3">Family</th><th rowspan="3">Model</th><th colspan="6">Φ drift</th></tr><tr><th colspan="2">Quasi-stationary</th><th colspan="2">Westward</th><th colspan="2">Eastward</th></tr><tr><th>N<sub>RFPs</sub></th><th>Rate</th><th>N<sub>RFPs</sub></th><th>Rate</th><th>N<sub>RFPs</sub></th><th>Rate</th><th>Velocity</th></tr></thead><tbody><tr><td rowspan="5">A_FM</td><td>A_FM-0000</td><td>0</td><td>—</td><td>4</td><td>-0.05</td><td>3</td><td>0.14</td><td>6.17</td></tr><tr><td>A_FM-0318</td><td>0</td><td>—</td><td>3</td><td>-0.06</td><td>4</td><td>0.07</td><td>5.40</td></tr><tr><td>A_FM-0677</td><td>1</td><td>*</td><td>2</td><td>-0.08</td><td>3</td><td>0.13</td><td>5.88</td></tr><tr><td>A_FM-0918</td><td>0</td><td>—</td><td>5</td><td>-0.07</td><td>2</td><td>0.07</td><td>4.86</td></tr><tr><td>A_FM-M</td><td>2</td><td>*</td><td>4</td><td>-0.03</td><td>1</td><td>0.28</td><td>3.25</td></tr><tr><td rowspan="5">ASD_FM</td><td>ASD_FM-0000</td><td>0</td><td>—</td><td>2</td><td>-0.11</td><td>2</td><td>0.04</td><td>2.87</td></tr><tr><td>ASD_FM-0097</td><td>0</td><td>—</td><td>3</td><td>-0.01</td><td>2</td><td>0.05</td><td>1.67</td></tr><tr><td>ASD_FM-0377</td><td>1</td><td>*</td><td>3</td><td>-0.09</td><td>3</td><td>0.04</td><td>3.40</td></tr><tr><td>ASD_FM-0683</td><td>1</td><td>*</td><td>3</td><td>-0.04</td><td>2</td><td>0.03</td><td>2.12</td></tr><tr><td>ASD_FM-M</td><td>0</td><td>—</td><td>2</td><td>-0.05</td><td>4</td><td>0.05</td><td>3.69</td></tr><tr><td rowspan="4">ASDI_FM</td><td>ASDI_FM-0000</td><td>2</td><td>*</td><td>3</td><td>-0.09</td><td>7</td><td>0.08</td><td>4.29</td></tr><tr><td>ASDI_FM-0188</td><td>2</td><td>-0.01</td><td>0</td><td>—</td><td>6</td><td>0.08</td><td>4.14</td></tr><tr><td>ASDI_FM-0425</td><td>1</td><td>*</td><td>5</td><td>-0.06</td><td>4</td><td>0.05</td><td>3.39</td></tr><tr><td>ASDI_FM-0962</td><td>1</td><td>*</td><td>4</td><td>-0.04</td><td>3</td><td>0.05</td><td>3.15</td></tr><tr><td rowspan="3">pfm9k</td><td>ASDI_FM-M</td><td>3</td><td>*</td><td>2</td><td>-0.12</td><td>5</td><td>0.05</td><td>3.34</td></tr><tr><td>pfm9k.1</td><td>1</td><td>*</td><td>6</td><td>-0.02</td><td>4</td><td>0.06</td><td>2.49</td></tr><tr><td>pfm9k.1a</td><td>1</td><td>*</td><td>6</td><td>-0.03</td><td>3</td><td>0.06</td><td>1.86</td></tr><tr><td rowspan="2">CALSXk</td><td>pfm9k.1b</td><td>1</td><td>*</td><td>4</td><td>-0.02</td><td>7</td><td>0.07</td><td>2.49</td></tr><tr><td>CALS3k.4b</td><td>2</td><td>*</td><td>9</td><td>-0.10</td><td>2</td><td>0.07</td><td>4.58</td></tr></tbody></table>

Φ drift is the displacement in the east-west direction. N<sub>RFPs</sub> is the number of RFPs. Rate is the angular velocity in °yr<sup>-1</sup>. Quarter of the weighted average rate is the critical value to distinguish quasi-stationary and drifting features. Velocity is the horizontal vector speed in km yr<sup>-1</sup>, * denotes absolute rate values smaller than 0.005°yr<sup>-1</sup>, NHKSH14 models apply only for the past 3 kyr and CALS3k.4b results are from Terra-Nova et al. (2015). The specific samples from the ensembles of models were chosen based on the same methodology used to construct Fig. 2. 0000 is the reference models, while M models are the means of ensembles of 1000 models.

the reference model (A_FM-0000) and lowest in the mean model (A_FM-M), with the former being almost twice larger (Table 1). The ASD_FM models (LHGT13) present similar number of RFPs with westward or eastward drift, with the rate of westward drift either larger or comparable (Table 1). Only the ASD_FM-0377 and ASD_FM-0683 models resolved RFPs with quasi-stationary behaviour. The velocity is highest for the mean model (ASD_FM-M) and lowest in the ASD_FM-0097 model, with the former being more than twice larger (Table 1). The ASD_FM family of models (LHGT13) presents either more RFPs with eastward drift or comparable drifts, with the ASD_FM-0188 model not resolving any RFPs with westward drift. The rate of westward drift is usually larger than the rate of eastward drift (Table 1). All models of this family resolved RFPs with quasi-stationary behaviour, with ASD_FM-M resolving more quasi-stationary RFPs than RFPs with westward drift. The velocity is highest for the reference model (ASD_FM-0000) and lowest in the ASD_FM-0962 model (Table 1).

The pfm9k.1 and pfm9k.1a models (NHKSH14) present more RFPs with westward drift, whereas the pfm9k.1b model resolves more RFPs with eastward drift. The eastward rate is at least twice larger than the westward rate for all models (Table 1). All models resolved one RFP with quasi-stationary behaviour. The velocity is practically the same for the pfm9k.1 and pfm9k.1b models and lower for the pfm9k.1b model (Table 1). Lastly, The CALS3k.4b model (KC11) presents significantly more RFPs with westward drift, as well as two RFPs with quasi-stationary behaviour (Table 1). The rate of westward drift is larger than the eastward drift (Terra-Nova et al. 2015).

Fig. 7 shows an example of westward drift of an RFP in the ASD_FM-M model (LHGT13). Over a period of 200 yr, the RFP

moved from 23°E to 7°E in an average rate of 0.08°yr<sup>-1</sup>. Two interesting aspects are worth noting. First, a cut of the magnetic equator is observed from 1480 AD to 1520 AD (Figs 7a and b). This cut shows a good example of how the algorithm performs. In Fig. 7(a), the magnetic equator reaches high latitudes forming a neck, which is cut in the next snapshot (Fig. 7b) forming an RFP in the Northern Hemisphere, which later disappears (Fig. 7c). Second, the westward drifting RFP resolved in the Southern Hemisphere is part of a pair. North of this RFP, a patch with intense normal flux appears (Figs 7a and b) and later it fades (Fig. 7c). This bipolar configuration is reminiscent of the present-day field below the South Atlantic (Jackson et al. 2000; Olsen et al. 2014).

In summary, there are significant similarities of the tracking results within each family of models, but the results differ substantially from one family to another. The reason for these differences may be that the field models are significantly different. Alternatively, these differences may arise because RFPs are weak features that are strongly sensitive to small changes in the field. In order to test these hypotheses we will compare the tracking of RFPs with tracking of NFPs which are strong features of the field.

## 4.2 Normal flux patches tracking

Here we compare the identification and tracking of NFPs and RFPs to assess the robustness of the results of RFPs among different families. Figs 8 and 9 show the tracking of NFPs for all mean models of LHGT13: A_FM-M, ASD_FM-M and ASD_FM-M, as well as pfm9k.1b (NHKSH14) and CALS3k.4b (KC11). Two persistent NFPs are clearly observed. One is located near longitude 0° between 0 AD and 1000 AD. This NFP shows both westward and eastward
<figure><img src="image_9.png" /><figcaption>(b) ASDI_FM-M 1520AD</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c) ASDI_FM-M 1680AD</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** As in Fig. 1 at (a) 1480 AD, (b) 1520 AD and (c) 1680 AD for the ASDI_FM-M model.</figcaption></figure>

drift in all models. The second appears around 1720 AD in almost all models near longitude 120°. This NFP shows little azimuthal drift. We observed in all models a proliferation of NFPs between longitudes 60°W and 120°W after 750 AD. We find significant similarity between the ASD_FM-M and ASD_FM-M models. For example, in both models one NFP with strong westward drift moves from about longitude 60°E to about longitude 120°W between 1000 BC and 20 AD. Some similarities are also observed between the CALS3k.4b and pfm9k.1b models, but to a lesser extent. No preferred longitudinal drift direction is seen for NFPs. In all models NFPs cluster around the tangent cylinder. Northern NFPs are more frequent than southern. The A_FM-M model has the lowest number of NFPs in the Southern Hemisphere. The pfm9k.1b model exhibits the oldest NFPs in the Southern Hemisphere, one of them persisting for 1500 yr near 60°W. The number of RFPs that migrate poleward and equatorward is similar. The fact that the strong NFPs are reproduced within model families but are distinct among different families, in qualitative agreement with the weaker RFPs, suggests that differences between the archaeomagnetic field models are relevant in tracking RFPs. For NFPs results within a given family see Figs S5 and S6 and their description in the Supporting Information section.

## 4.3 Persistent locations of reversed flux patches

We now compare the persistent locations of RFPs in the archaeomagnetic field models with those of the gufm1 model (Jackson *et al.* 2000). Fig. 10 shows the results from Terra-Nova *et al.* (2015) for the tracking of RFPs for the period between 1840 AD and 1990 AD of gufm1 (Jackson *et al.* 2000). The majority of RFPs in gufm1 drift westward (Terra-Nova *et al.* 2015). To analyse the persistent locations of RFPs we plotted their histograms in 20° bins (see Figs 11 and 12). It is worth noting that here we count RFPs positions at each snapshot, that is, all RFPs are considered regardless of their longevity.

In Fig. 11, we show histograms of the ASDI_FM family of models. We find two regions of persistent locations of RFPs, between longitudes 80°W and 40°E and between 100°E and 160°E, with the latter containing the largest amount of RFPs. In co-latitude we note three peaks of RFPs, between 60° and 80°, between 100° and 120°, and between 140° and 160°. The first and the second regions concentrate the largest amounts of RFPs. No RFPs are found close to polar regions, and only a few are found in equatorial regions. For results of the families of models A_FM and ASD_FM see Figs S7 and S8 and their description in the Supporting Information section.

Fig. 12 shows three regions of persistent locations of RFPs in longitude in all pfm9k models, between 120°W and 100°W, between 40°W and 40°E, and between 120°E and 160°E. The second region contains the largest amount of RFPs in all models except in model pfm9k.1a where the third region has a number of RFPs comparable to the second region. The RFPs are observed more frequently in the Northern Hemisphere. Almost no RFP is resolved near the equator. All pfm9k models present the highest quantity of RFPs between co-latitudes 40° and 60°. Only pfm9k.1a finds RFPs in the northern polar region. Overall the results show strong similarity among the pfm9k models.

The results of Fig. 12 for CALS3k.4b and gufm1 are very similar simply because CALS3k.4b is strongly constrained by gufm1 (KC11) and most of its RFPs appear in the historical period (Terra-Nova *et al.* 2015). Two regions contain all RFPs from gufm1, in the Atlantic Hemisphere between longitudes 80°W and 60°E and in the west Pacific between longitudes 140°E and 180°E. The gufm1 RFPs exhibit larger concentration in the Southern Hemisphere. The co-latitude highest peak in both gufm1 and CALS3k.4b is at mid-latitudes of the Southern Hemisphere between 120° and 140°. Both models exhibit RFPs in the northern polar region. The CALS3k.3 model, which is less constrained by gufm1 (KC11), shows in longitude an extended region of recurrence of RFPs between 120°W and 60°E, with a peak between 60°W and 40°W. RFPs are resolved along all co-latitudes except close to the polar regions.

The density of RFPs is defined as the number of RFPs within one longitude or co-latitude bin divided by the total number of RFPs of each model. Fig. 13 shows the density of RFPs in longitude and co-latitude for models A_FM-M, ASD_FM-M, ASD_FM-M, pfm9k.1b, CALS3k.4b, CALS3k.3 and gufm1. We note two longitudinal regions that exhibit RFPs concentration (Fig. 13a). The first is between longitudes 80°W and 50°E. The second is between longitudes 90°E and 170°E, which is the region where the density of RFPs in the ASD_FM-M and ASD_FM-M models is higher. In Fig. 13(b), only models ASD_FM-M, pfm9k.1b, and CALS3k.3 have more RFPs in the Northern Hemisphere while the other models find more RFPs in the Southern Hemisphere. Low density of RFPs are observed in the equatorial region in all models.

If there is strong enough mantle control on the geodynamo (e.g. Amit *et al.* 2015, and references therein), possibly its signature may
<figure><img src="image_12.png" /><figcaption>**Figure 8.** As in Fig. 3 for normal flux patches for the A_FM-M, ASD_FM-M, ASDI_FM-M, pfm9k.1b and CALS3k.4b models.</figcaption></figure>

be visible already in the historical period (Constable *et al.*, in preparation). The gufm1 histograms peaks in longitude and co-latitude (Figs 12i and j, and solid black line in Fig. 13) may represent persistent sites of RFPs. We thus compare gufm1 histograms with other archeomagnetic models in Fig. 13. The agreement with CALS3k.4b

is not considered for the reasons discussed above. In longitude, the A_FM-M, ASD_FM-M, pfm9k.1b and CALS3k.3 models overlap the range of gufm1 in two regions, but these regions are almost all distinctive among the models. For A_FM between longitudes 80°W and 60°W and between 140°E and 160°E, for ASD_FM-M
<figure><img src="image_13.png" /><figcaption>**Figure 9.** As in Fig. 4 for normal flux patches for the A_FM-M, ASD_FM-M, ASDI_FM-M, pfm9k.1b and CALS3k.4b models. The dotted lines denote the tangent cylinder.</figcaption></figure>

and pfm9k.1b between longitudes 0° and 20°E and between 140°E and 160°E, and for CALS3k.3 between longitude 80°W and 40°W and between 0° and 60°E. The ASDI_FM-M model correlates with gufm1 in particular in three regions, between longitudes

0° and 20°E, between 140°E and 160°E, and between 80°W and 40°W. In co-latitude, the ASDI_FM-M model shows some agreement with the gufm1 model, yet it is not strong. In summary, the closest longitudinal agreement with gufm1 among these models is
<figure><img src="image_14.png" /><figcaption>Figure 10. (a) As in Fig. 3, and (b) as in Fig. 4, both for the historical field model gufm1 from 1840 AD to 1990 AD.</figcaption></figure>

thus obtained by ASDI_FM-M model. The correlation coefficients between the longitude density of RFPs in each archaeomagnetic field model and gufm1 (Table 2) indeed show that ASDI_FM-M is the best correlated. Perhaps surprisingly the ASD_FM-M model mimics less good the results of gufm1 than A_FM-M. In addition, CALS3k.3 better correlates with gufm1 than CALS3k.4b.

Next we compare the locations of RFPs in our preferred ASD_FM-M model with a model of lower mantle heterogeneity. Seismic shear waves velocity ($v_s$) anomalies may approximate the temperature heterogeneity in the lowermost mantle (e.g. Amit et al. 2015). Gubbins (2003) compared integrated magnetic flux with integrated seismic anomalies. Similarly, we compare the longitudinal distribution of RFPs with the longitudinal distribution of $\delta v_s/v_s$ which is given by

$$V_s(\phi) = \int_{0}^{\pi} \delta \ln v_s(\phi, \theta) \sin(\theta) d\theta. \qquad (3)$$

We also examine the distributions in co-latitude by

$$V_s(\theta) = \int_{-\pi}^{\pi} \delta \ln v_s(\phi, \theta) d\phi. \qquad (4)$$

We consider both the full $\delta v_s/v_s$ as well as only its negative part (Fig. 13) which is given in longitude by

$$V_s^{-}(\phi) = \begin{cases} \int_{0}^{\pi} \delta \ln v_s(\phi, \theta) \sin(\theta) d\theta & \text{if } \delta v_s/v_s < 0 \\ 0 & \text{if } \delta v_s/v_s > 0 \end{cases} \qquad (5)$$

and co-latitude by

$$V_s^{-}(\theta) = \begin{cases} \int_{-\pi}^{\pi} \delta \ln v_s(\phi, \theta) d\phi & \text{if } \delta v_s/v_s < 0 \\ 0 & \text{if } \delta v_s/v_s > 0 \end{cases} \qquad (6)$$

In addition, for $V_s$ we set the plotted range to show only negative integrated $\delta v_s/v_s$. For visualization purposes we plot $-V_s$ (and likewise for the negative part) since we want to examine the correlation between hot mantle structures (associated with negative $v_s$ anomalies) and the density of RFPs. The seismic anomalies are based on the tomography model of Masters et al. (2000) at the lowermost mantle. Since ASDI_FM-M is the archaeomagnetic field that best resembles gufm1 we compare it with the seismic distributions (Fig. 13). In longitude, the seismic anomalies are concentrated in two regions (Fig. 13a). The negative part of $V_s$ (pink curve) is around Africa between ~30°W and ~60°E, and in the middle of the Pacific between ~150°E and ~120°W. The $V_s^-$ (dark yellow curve) is in similar regions (though more spread) between 40°W and 110°E and between 130°E and 90°W. Our preferred ASDI_FM-M model also has two regions of high density of RFPs, which are shifted to the west with respect to the seismic regions of $V_s$ by about 40° (Africa) and 55° (Pacific). In co-latitude, seismic distributions have remarkable correlation with the peaks in density of RFPs from the ASDI_FM-M model (Fig. 13b), with the three curves indicating negative $\delta v_s$ and more RFPs at low-latitudes of the Southern Hemisphere with a peak near 110°, that is, no shift between seismic distributions and RFPs curves is observed in co-latitude. Indeed, the correlation coefficients between $V_s$ and ASDI_FM-M as well as between $V_s^-$ and
Q (number of patches)

<figure><img src="image_15.png" /><figcaption>(c) ASDI_FM_0188</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(e) ASDI_FM_0425</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(g) ASDI_FM_0962</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(i) ASDI_FM_M</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>$\phi(°)$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(d) ASDI_FM_0188</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(f) ASDI_FM_0425</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(h) ASDI_FM_0962</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(j) ASDI_FM_M</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>$\theta(°)$</figcaption></figure>

**Figure 11.** Histograms of reversed flux patches occurrences for the last 3 kyr in longitude (left) and co-latitude (right) for the ASDI_FM family of models. Q is the number of reversed flux patches at a position range integrated over all snapshots. The histograms are constructed with bins of 20°.

ASDI_FM-M are significantly lower in longitude (0.28 and 0.27, respectively) than in co-latitude (0.51 and 0.56, respectively).

## 4.4 Archaeomagnetic field models for the last 9 kyr

In order to investigate the longer term behaviour of RFPs we used the pfm9k models, which enable to expand the timescale of our analysis to the past 9 kyr. Figs 14 and 15 show the results of RFPs tracking of the pfm9k family of models (NHKSH4) for the past 9 kyr. The three models share common features, in particular the models pfm9k.1 and pfm9k.1a. These two models resolved more RFPs than pfm9k.1b, with westward being the more common drift direction (Fig. 14). In contrast, pfm9k.1b shows more balance between westward and eastward drifts. Several RFPs are resolved in all models, for example, near longitude 30°E between 3150 BC and 2800 BC, and near longitude 90°W between 1500 BC and 750 BC. However, the three models exhibit different drift directions for the most robust RFPs. All models show similar number of RFPs that move equatorward and poleward (Fig. 15). Some RFPs alternate between the two motions, such as the long-lived RFP in the pfm9k.1b model at around co-latitude 100° between 5000 BC and 4000 BC. RFPs are prominent near co-latitude 60° in all models. Two long-

lived RFPs near co-latitude 60° move northward between 1500 BC and 750 BC. In all models RFPs are dominantly in the Northern Hemisphere. It is worth noting that the two RFPs that are resolved in all models do not necessarily show the same motion behaviour (Figs 14 and 15).

In Fig. 16, we search persistent locations of tracked RFPs in the pfm9k family of models for the past 9 kyr. In longitude, RFPs appear in all models in three regions: between 120°W and 60°W, between 40°W and 60°W and between 120°E and 160°E. The second region has the largest amount of RFPs in all models. Some RFPs appear at other longitudes, for example in the pfm9k.1a model between longitudes 100°E and 160°E. More RFPs are found in the Northern Hemisphere with the largest peak between co-latitudes 40° and 60° in all models. The pfm9k.1 and pfm9k.1a models show two regions of RFPs in the Southern Hemisphere, between 100° and 120° and between 160° and 180°. The pfm9k.1b model only shows RFPs in the former region. RFPs are found in the equatorial region in the pfm9k.1b model more than in the other two models. Both longitude and co-latitude results of Fig. 16 show a high similarity among all models.

The agreement between the spatial distribution in Fig. 16 for the past 9 kyr and in Figs 12(a)–(f) for the past 3 kyr is striking. Unlike
Q (number of patches)

<figure><img src="image_25.png" /><figcaption>**Figure 12.** As in Fig. 11 for the pfm9k and CALS3K families of models, and the gufm1 model. Note that the pfm9k histograms apply for the past 3 kyr only, and that gufm1 histograms apply for 1840–1990 AD.</figcaption></figure>

the agreement between CALS3k.4b and gufm1 that arises from the emergence of RFPs in the former mostly during the historical period (Figs 5d, 6d and 10), in the pfm9k models RFPs appear throughout the entire 9 kyr period in a rather uniform way (Figs 14 and 15). This agreement is therefore a strong evidence that a period of 3 kyr already captures the most important aspects of the core kinematics. We elaborate on this point in the discussion section below.

## 5 DISCUSSION

Uncertainties in archaeomagnetic field models due to low temporal and spatial resolution may cast doubt on the tracking and even on the existence of RFPs. Terra-Nova *et al.* (2015) performed sensitivity tests using low-pass filters applied to CALS3k.4b (KC11). Their results are not sensitive to model damping, indicating that RFPs are robust in CALS3k.4b (Terra-Nova *et al.* 2015). Nevertheless assessing RFPs robustness over several archaeomagnetic field models

is a crucial sensitivity test because it examines more rigorously the effects of differences in data sources, data selection and treatment as well as in modelling strategies and the consequences of these choices on RFPs identification and tracking in the models.

LHGT13 used three types of data sets (see Supporting Information Table S1) with similar data treatment and modelling strategies to build three different families of models. The A_FM family resolved less RFPs than the ASD_FM and ASDI_FM families, possibly due to the lower geographic coverage of data in A_FM. It also had the highest velocity of RFPs (Table 1). Conversely, ASDI_FM had the largest number of RFPs probably due to its highest spatiotemporal resolution obtained by the inclusion of sedimentary intensity data. These results corroborate those of Panovska *et al.* (2015) who showed that without sediment data the Holocene field variations in the Southern Hemisphere are not well constrained by archaeomagnetic field models.

Data uncertainties prescribed by LHGT13 are lower than those of NHKSH14, although data with unpublished uncertainties were
<figure><img src="image_26.png" /><figcaption>**Figure 13.** Density of RFPs in longitude (a) and co-latitude (b) for models A_FM-M (dashed grey line), ASD_FM-M (dashed green line), ASDI_FM-M (solid blue line), pfm9k.1b (dashed cyan line), CALS3k.4b (dashed red line), CALS3k.3 (dashed black line) and gufm1 (solid black line). The curves are constructed based on the histograms in Figs S7, S8, 11 and 12 normalized by the total number of RFPs at each model. Also included are the density of $-V_s$ (solid thick pink line) and $-V_s^-$ (solid thick dark yellow line) calculated based on the model of Masters *et al.* (2000) truncated at spherical harmonic degree 6 (see eqs 3–6).</figcaption></figure>

<table><caption>**Table 2.** Correlation coefficients *C* between the longitude density of RFPs in each archaeomagnetic field model and gufm1.</caption><thead><tr><th>Models</th><th><i>C</i></th></tr></thead><tbody><tr><td>A_FM-M</td><td>0.21</td></tr><tr><td>ASD_FM-M</td><td>0.05</td></tr><tr><td>ASDI_FM-M</td><td>0.35</td></tr><tr><td>pmf9k.1b</td><td>0.02</td></tr><tr><td>CALS3k.4b</td><td>0.25</td></tr><tr><td>CALS3k.3</td><td>0.30</td></tr></tbody></table>

strongly penalized by LGHT13. Because LHGT13 provided ensembles of models, it was possible to investigate variations in RFPs behaviour only due to magnetic and age uncertainties on the data. Panovska *et al.* (2015) also found that determination of robust field

structures and variations in the models are strongly affected by uncertainty assignments which give weight to different data in the modelling procedure. NHKSH14 built models with larger $\Delta t$, hence the choice of NHKSH14 for larger uncertainties, and with larger spherical harmonic truncation $n_{max} = 10$ versus $n_{max} = 5$ of LGHT13 (Supporting Information Table S2). NHKSH14 gave more weight to archaeomagnetic data which is scarce in the Southern Hemisphere, thus recent RFPs resolved with the pfm9k family of models were usually at the Northern Hemisphere as in the A_FM family. Sedimentary ages uncertainties in NHKSH14 and LGHT13 models do not seem to play a major role in resolving RFPs position, but may influence the drift direction, as suggested by the similarity between models pfm9k.1 and pfm9k.1a in Figs 5, 6, 14 and 15. The reference model influences more than the sediment data uncertainties assignments, as also shown by Panovska *et al.* (2015). The pfm9k.1b model resembles the other NHKSH14 models only for the past 5 kyr
<figure><img src="image_27.png" /><figcaption>**Figure 14.** As in Fig. 3 for the pfm9k family of models for the past 9 kyr.</figcaption></figure>

(Figs 14 and 15), indicating less robustness in results of NHKSH14 before 3000 BC, when significantly less archaeomagnetic data are available.

Supporting Information Table S2 compares modelling strategies. KC11 used the lowest temporal knot spacing and relied on regularization parameters to fill the gap of missing data, whereas LHGT13 used a larger temporal knot spacing to have always data between knots. NHKSH14 applied the largest knot spacing but also the largest temporal window. LHGT13 used a knot of 40 yr to account for the dominant timescales until $n_{max} = 5$, which turns out to be sufficient for us to identify RFPs. It is important to stress that the effort made by LHGT13 to use outliers as source of data by reducing their weight has led to the conservation of the geographic coverage of data, in particular in the Southern Hemisphere, thus allowing to better resolve RFPs in this hemisphere.

The similarities within each family is tested by the recurrence of RFPs position and direction of displacement in longitude and co-latitude. The similarity is larger for the ASD_FM, ASDI_FM and pfm9k families of models than among the A_FM family of models. It indicates that inclusion of sedimentary data gives more stability and robustness to the models due to the better geographic coverage. Our analysis also shows that the direction of displacement is much less robust than the position of RFPs. Long-lived RFPs that appear in a given model will always appear in the other models of the same

family, at least as a region of proliferation of RFPs at the position of the long-lived RFP.

Comparing tracking among families shows that only a few RFPs are commonly resolved in the ASD_FM, ASDI_FM and pfm9k models, and even those exhibit different drift directions. The similarities among families of models are therefore limited to persistent locations of RFPs rather than similar mobility. None of the three families of models ASD_FM, ASDI_FM and pfm9k showed clear similarity with results from CALS3k.4b or A_FM-M. Only gufm1 and CALS3k.4b showed clearly preferred westward drift.

Obviously, caution is required in interpreting archaeomagnetic field models. In particular, RFPs are weak magnetic features, and as such they are sensitive to the choices made by field modellers. Fig. 7 shows an example of resolution issues. An intrusion of the magnetic equator is cut, forming a transient RFP. RFPs may be formed by toroidal field expulsion by fluid upwelling (Bloxham 1986), in which case an RFP is associated with an adjacent NFP. But no intense NFP pairs with the northern RFP in Fig. 7. It might be that this RFP is only resolved due to a resolution bias of the archaeomagnetic field model. We therefore complement our RFPs analysis with that of NFPs which are more intense and hence expected to be more robust. All models identify NFPs around the tangent cylinder. The Northern Hemisphere NFPs are more recurrent and clearly move more in longitude than in co-latitude. As with
<figure><img src="image_28.png" /><figcaption>Figure 15. As in Fig. 4 for the pfm9k family of models for the past 9 kyr.</figcaption></figure>

RFPs, the tracking of NFPs is much more similar within one family than among different families of models. Considerable differences in NFPs tracking among families of models are observed. Although there is somewhat more similarity among different families in NFPs tracking than in the RFPs tracking, results are qualitatively similar, and differences between NFPs even within one family lead us to argue that differences in RFPs tracking are mostly associated to different data sets and modelling strategies rather than to the bias of the tracking due to the small scale of RFPs.

To evaluate RFPs persistent locations we used histograms of their longitude and co-latitude. Among models of the same family, the histograms showed high similarity. We note similarities of position of RFPs between CALS3k.4b, A_FM-M and pfm9k.1b indicating that the strong weighting of archaeomagnetic data in these models provides a strong constraint. As in the tracking results, histograms within a family show higher similarity than between different families.

In summary, we show that incorporating more data to the models (bearing in mind that these data have higher uncertainties) increases the number of identified NFPs and RFPs. Therefore, modelling strategies that opt for the inclusion of the highest temporal and spatial coverage of data are more successful in producing flux patches at centennial and millennial timescales.

A high-resolution model is required to resolve the field morphology, and a long period is essential to cover as many timescales

as possible. We hypothesize that the resolution problem is more severe, thus gufm1 is arguably more reliable than archaeomagnetic field models. Using a palaeomagnetic SV index Constable *et al.* (in preparation) found that both the historical and the Holocene fields are weaker and more time variable in the Southern Hemisphere. If we assume that mantle control is significant at timescales already covered by gufm1, we may apply the comparison of all archaeomagnetic field models with gufm1 to select a candidate for the most reliable archaeomagnetic field model.

We searched for the archaeomagnetic field model that is most similar to gufm1. We excluded comparison with CALS3k.4b because it is strongly constrained by gufm1 and has most of its RFPs during the historical period. In contrast, in CALS3k.3 the RFPs are more distributed over the entire period. Here we did not compare preferred drift directions which are less robust, but instead we compared locations of persistent RFPs which are robust within the same family and in some cases among models of different families (see Section 4). We found that ASDI_FM-M is the model that most resembles gufm1, in particular in terms of longitude of RFPs. This model has strong damping, uses bootstraps, preserves outliers and uses all types of data.

Next, we examined whether the ASDI_FM-M model can be interpreted in terms of core–mantle thermal coupling. Heterogeneous CMB conditions may provide statistical preference to the positions
<figure><img src="image_29.png" /><figcaption>**Figure 16.** As in Fig. 11 for the pfm9k family of models for the past 9 kyr.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 17.** Schematic illustration of reversed flux patch emergence on the CMB due to flux expulsion of a toroidal field line by fluid upwelling. Red represents normal flux and blue reversed flux in the Southern Hemisphere. $\delta$ is the shift between the centre of fluid upwelling and the centre of an RFP (or an NFP).</figcaption></figure>

of RFPs (Tarduno *et al.* 2015). Usually, seismic lowermost mantle anomalies are interpreted as thermal heterogeneities. Negative heat flux anomalies localize core fluid upwelling (Gubbins 2003). We compared the seismic shear wave velocity anomalies at the lowermost mantle (Masters *et al.* 2000) with RFPs density in our preferred archaeomagnetic field model ASDI_FM-M. We found similar trends with a shift between maxima values in longitude and a strong correlation in co-latitude. These results may indicate mantle control on RFPs locations. Fig. 17 shows a schematic diagram of toroidal field line expulsion by fluid upwelling in the Southern Hemisphere. The shift $\delta$ between a centre of fluid upwelling and a centre of an RFP is roughly half of the fluid upwelling structure size. Possible shifts in co-latitude may be explained by the same scheme, with the axis in the north-south direction.

Since the shift is found in longitude with no shift in co-latitude, it may imply that toroidal field lines near RFPs are oriented azimuthally (Fig. 17). In addition, the amount of shift suggests that RFPs are formed by upwelling structures with a width of about 80° (Africa) and 110° (Pacific). Numerical dynamo simulations also show that the toroidal field on the outer boundary has a significant east-west component (Glatzmaier & Roberts 1995; Olson *et al.* 1999). Aubert *et al.* (2008) termed the intense localized azimuthal toroidal field structures as ‘magnetic cyclones’. These features are formed in their dynamo models at the equatorial part of columnar flow cyclones. Magnetic cyclones expel toroidal field through the CMB, resulting in pairs of NFPs/RFPs. These patches are relatively short-lived due to the dominant diffusion in these regions (Aubert *et al.* 2008).
It is worth considering the longer term archaeomagnetic field models pfm9k in the early epoch. More RFPs are found in the Southern Hemisphere compared to the last 3 kyr, since these models are more strongly controlled by sedimentary data due to the scarcity of archaeomagnetic field data at ancient times. Histograms for the last 3 kyr and for the last 9 kyr are very similar for each NHKSH14 model (Figs 12 and 16). We suggest that this is evidence that 3 kyr is already enough time to detect with some confidence the mantle control on the Earth's magnetic field, in particular to identify statistically significant recurrence of RFPs at specific regions. These results may therefore provide an upper bound for the magnetic boundary layer thickness in the Earth's outer core (Gubbins 1996). If an effective magnetic diffusion time obeys $\tau_{\eta} < 3000$ yr, and since radial diffusion is expected to be much more dominant than tangential diffusion in the SV (Amit & Christensen 2008) then:

$$ \frac{\mathcal{L}_r^2}{\eta} < 3000. \qquad (7) $$

For $\eta \sim 0.5 - 1 \text{ m}^2\text{s}^{-1}$ (Pozzo et al. 2012; Zhang et al. 2015) we get $\mathcal{L}_r < 220 - 310$ km for the typical radial lengthscale of the field just below the CMB, which represents the magnetic boundary layer thickness.

Amit & Christensen (2008) used numerical dynamos to evaluate the ratio of radial to horizontal diffusion $D_f$. For the observed large-scale field they extrapolated $D_f \sim 260$. Because this ratio scales like $D_f \sim (\frac{\mathcal{L}_r}{\eta})^2$, and assuming a horizontal lengthscale of $\mathcal{L} \sim 1000$ km, the radial lengthscale is $\mathcal{L}_r \sim 60$ km. Chulliat & Olsen (2010) used two snapshots of satellites data to estimate the magnetic boundary layer thickness based on the intensification of reversed flux surrounded by a null-curve as:

$$ \mathcal{L}_r = \sqrt{\eta \Delta t \frac{B_r}{\Delta B_r}} \qquad (8) $$

where $\Delta t$ is the time difference between the snapshots, $\Delta B_r$ is the variation of magnetic flux within the null-curve during the period, and $B_r$ is the mean magnetic flux. Using $\Delta t = 25$ yr, $B_r \sim 200$ µT, $\Delta B_r \sim 100$ µT (Chulliat & Olsen 2010) and $\eta \sim 0.5 - 1 \text{ m}^2\text{s}^{-1}$ (Pozzo et al. 2012; Zhang et al. 2015), eq. (8) gives $\mathcal{L}_r \sim 30-40$ km. In summary, our estimations of the upper bound of the magnetic boundary layer thickness are $\sim 4-8$ times larger than the estimations of Amit & Christensen (2008) and Chulliat & Olsen (2010) respectively, further supporting the likelihood that a period of 3 kyr is sufficiently long to capture meaningful kinematics at the top of the core.

## 6 CONCLUSIONS

We identified and tracked RFPs in 20 different archaeomagnetic field models. We applied different analyses to extract robust results. Based on these results we inferred some aspects of outer core physics. Our main findings are:

(i) RFPs positions and motions are sensitive to the type of data source because it affects the geographic data coverage.

(ii) Similarities in RFPs behaviour appear within a family of archaeomagnetic field models, but models from different families yield distinctive RFPs behaviour, suggesting that modelling strategy is more influential than data set.

(iii) Similarities among archaeomagnetic field models are better observed by recurrence of positions of RFPs. Motion direction seems to be strongly sensitive to the assignment of uncertainties to

the data, thus we cannot confirm the persistent westward drift in the outer core on millennial timescales found in the previous analysis of Terra-Nova et al. (2015).

(iv) Preferred RFPs positions in the ASDI_FM-M model are shifted in longitude and correlated in co-latitude with respect to negative lower mantle seismic anomalies, suggesting that toroidal field lines below the CMB are oriented azimuthally. We also infer large fluid upwelling structures of $80^{\circ}-110^{\circ}$ at the top of the core.

(v) The similarity between 3 kyr and 9 kyr position histograms at each of the pfm9k models suggests that a 3 kyr period is already long enough to detect with confidence the mantle control on the Earth's field, allowing to estimate an upper bound of 220-310 km for the magnetic boundary layer thickness below the CMB.

Clearly new Holocene geomagnetic field data are always needed to better constrain archaeomagnetic field models, especially in the Southern Hemisphere (Poletti et al. 2016), Central Asia, oceans and high-latitude regions, as well as for epochs that precede 500 BC (e.g. Donadini et al. 2009; Brown et al. 2015a,b). Incorporation of recent local time-series (e.g. Brown et al. 2015a,b; Tarduno et al. 2015, and references therein) and regional models (e.g. Hellio et al. 2014) may improve archaeomagnetic field models. These improved models may lead to more reliable identification and tracking of archaeomagnetic field structures, and consequently better understanding of core-mantle coupling and core dynamics on millennial timescales.

## ACKNOWLEDGEMENTS

FT-N acknowledges The National Council for Scientific and Technological Development (CNPq/Brazil) for grant 206997/2014-0. FT-N and HA were partly supported by the Centre National des Études Spatiales (CNES). GAH thanks CAPES (grant AUXPE 2043/2014) and CNPq/Brazil (grant 454609/2014-0). RIFT is supported by CNPq/Brazil PQ grant 304934/2014-3. We are grateful to Richard Holme, Monika Korte and an anonymous reviewer for their constructive comments that improved the manuscript.

## REFERENCES

Amit, H. & Christensen, U.R., 2008. Accounting for magnetic diffusion in core flow inversions from geomagnetic secular variation, *Geophys. J. Int.*, **175**, 913-924.

Amit, H. & Olson, P., 2006. Time-average and time-dependent parts of core flow, *Phys. Earth planet. Inter.*, **155**, 120-139.

Amit, H., Korte, M., Aubert, J., Constable, C.G. & Hulot, G., 2011. The time-dependence of intense archeomagnetic flux patches, *J. geophys. Res.*, **116**, B12106, doi:10.1029/2011JB008538.

Amit, H., Choblet, G., Olson, P., Monteux, J., Deschamps, F., Langlais, B. & Tobie, G., 2015. Towards more realistic core-mantle boundary heat flux patterns: a source of diversity in planetary dynamos, *Prog. Earth Planet. Sci.*, **2**(26), doi:10.1186/s40645-015-0056-3.

Aubert, J., Aurnou, J. & Wicht, J., 2008. The magnetic structure of convection-driven numerical dynamos, *Geophys. J. Int.*, **172**, 945-956.

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

Backus, G.E., 1988a. Comparing hard and soft prior bounds in geophysical inverse problems, *Geophys. J. R. astr. Soc.*, **92**, 125-142.

Backus, G.E., 1988b. Bayesian inference in geomagnetism, *Geophys. J. R. astr. Soc.*, **94**, 249-261.

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

Bloxham, J. & Jackson, A., 1991. Fluid flow near the surface of the Earth's outer core, *Rev. Geophys.*, **29**, 97-120.
Bloxham, J. & Jackson, A., 1992. Time-dependent mapping of the magnetic field at the core-mantle boundary, *J. geophys. Res.*, **97**, 19 537–19 563.

Brown, M.C., Donadini, F., Korte, M., Nilsson, A., Korhonen, K., Lodge, A., Lengyel, S.N. & Constable, C.G., 2015a. GEOMAGIA50.v3: 1. General structure and modifications to the archeological and volcanic database, *Earth Planets Space*, doi:10.1186/s40623-015-0232-0.

Brown, M.C. *et al.*, 2015b. GEOMAGIA50.v3: 2. A new paleomagnetic database for lake and marine sediments, *Earth Planets Space*, doi:10.1186/s40623-015-0233-z.

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

Chulliat, A. & Olsen, N., 2010. Observation of magnetic diffusion in the Earth's outer core from Magsat, Orsted and CHAMP data, *J. geophys. Res.*, **115**, doi:10.1029/2009JB006994.

Courtillot, V., Ducruix, J. & Le Mouel, J.-L., 1978. Sur une accélération récente de la variation séculaire du champ magnétique terrestre, *C. R. Acad. Sci.*, *Paris D*, **287**, 1095–1098.

Donadini, F., Korhonen, K., Riisager, P. & Pesonen, L.J., 2006. Database for Holocene geomagnetic intensity information, *EOS, Trans. Am. geophys. Un.*, **87**, 137–138.

Donadini, F., Korte, M. & Constable, C.G., 2009. Geomagnetic field for 0-3ka: 1. new data sets for global modeling, *Geochem. Geophys. Geosyst.*, **10**, doi:10.1029/2008GC002295.

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.

Finlay, C.C. & Jackson, A., 2003. Equatorially dominated magnetic field change at the surface of Earth's core, *Science*, **300**, 2084–2086.

Finlay, C.C., Jackson, A., Gillet, N. & Olsen, N., 2012. Core surface magnetic field evolution 2000–2010, *Geophys. J. Int.*, **189**, 761–781.

Finlay, C.C., Olsen, N. & Tøffner-Clausen, L., 2015. DTU candidate field models for IGRF-12 and the CHAOS-5 geomagnetic field model, *Earth Planets Space*, **67**(114), doi:10.1186/s40623-015-0274-3.

Finlay, C.C., Aubert, J. & Gillet, N., 2016. Gyre-driven decay of the earth's magnetic dipole. *Nat. Commum.*, **7**, 10422, doi:10.1038/ncomms10422.

Fisher, R., 1953. Dispersion on a sphere, *Proc. R. Soc. A*, **217**, 295–035.

Glatzmaier, G.A. & Roberts, P.H., 1995. A three-dimensional self-consistent computer simulation of a geomagnetic field reversal, *Nature*, **377**, 203–209.

Gubbins, D., 1975. Can the Earth's magnetic field be sustained by core oscillations?, *Geophys. Res. Lett.*, **2**, 409–412.

Gubbins, D., 1987. Mechanism for geomagnetic polarity reversals, *Nature*, **326**, 167–169.

Gubbins, D., 1996. A formalism for the inversion of geomagnetic data for core motions with diffusion, *Phys. Earth planet. Inter.*, **98**, 193–206.

Gubbins, D., 2003. Thermal core-mantle interactions: theory and observations, in *Earth's Core: Dynamics, Structure and Rotation*, eds Dehant, V., Creager, K., Karato, S. & Zatman, S., American Geophysical Union, doi:10.1029/GD031p0163.

Gubbins, D. & Bloxham, J., 1986. Geomagnetic field analysis IV - Testing the frozen-flux hypothesis, *Geophys. J. R. astr. Soc.*, **84**, 139–152.

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

Hellio, G., Gillet, N., Bouligand, C. & Jault, D., 2014. Stochastic modelling of regional archaeomagnetic series, *Geophys. J. Int.*, **199**, 931–943.

Holme, R., 2015. Large-scale flow in the core, in *Treatise on Geophysics*, 2nd edn, vol. 8, pp. 107–130, ed. Olson, P., Elsevier.

Hongre, L., Hulot, G. & Khokhlov, A., 1998. An analysis of the geomagnetic field over the past 2000 years, *Phys. Earth planet. Inter.*, **106**, 311–335.

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

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

Jackson, A., Constable, C.G., Walker, M.R. & Parker, R.L., 2007. Models of Earth's main magnetic field incorporating flux and radial vorticity constraints, *Geophys. J. Int.*, **171**, 133–147.

Korhonen, K., Donadini, F., Riisager, P. & Pesonen, L.J., 2008. GEOMAGIA50: an archeointensity database with PHP and MySQL, *Geochem. Geophys. Geosyst.*, **9**, Q04029, doi:10.1029/2007GC001893.

Korte, M. & Constable, C.G., 2011. Improving geomagnetic field reconstructions for 0–3 ka, *Phys. Earth planet. Inter.*, **188**(3–4), 247–259.

Korte, M. & Holme, R., 2010. On the persistence of geomagnetic flux lobes in global field models, *Phys. Earth planet. Inter.*, **182**, 179–186.

Korte, M., Donadini, F. & Constable, C.G., 2009. The geomagnetic field for 0–3 ka: 2. A new series of time-varying global models, *J. geophys. Res.*, **10**, Q06008, doi:10.1029/2008GC002297.

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

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

Lodge, A. & Holme, R., 2009. Towards a new approach to archaeomagnetic dating in Europe using geomagnetic field modelling, *Archaeometry*, **51**, 309–322.

Mandea, M., Holme, R., Pais, A., Pintheiro, K., Jackson, A. & Verbanac, G., 2010. Geomagnetic jerks: rapid core field variations and core dynamics, *Space Sci. Rev.*, **155**, 147–175.

Masters, G., Laske, G., Bolton, H. & Dziewonski, A., 2000. The relative behavior of shear velocity, bulk sound velocity, and compressional velocity in the mantle: implications for chemical and thermal structure, in *Earths Deep Interior*, vol. 117, eds Karato, S., Forte, A., Liebermann, R., Masters, G. & Stixrude, L., American Geophysical Union, doi:10.1029/GM117p0063.

Maus, S. *et al.*, 2005. The 10th generation international geomagnetic reference field, *Geophys. J. Int.*, **161**, 561–565.

Merrill, R., McElhinny, M. & McFadden, P., 1998. *The Magnetic Field of the Earth: Paleomagnetism, the Core, and the Deep Mantle*, Academic Press.

Muscheler, R., Beer, J., Wagner, G., Laj, C., Kissel, C., Raisbeck, G.M., Yiou, F. & Kubik, P.W., 2004. Changes in the carbon cycle during the last deglaciation as indicated by the comparison of $^{10}$Be and $^{14}$C records, *Earth planet. Sci. Lett.*, **219**, 325–340.

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.

O'Brien, M.S., Constable, C.G. & Parker, R.L., 1997. Frozen-flux modelling for ephocs 1915 and 1980, *Geophys. J. Int.*, **171**, 434–450.

Olsen, N., Luehr, H., Sabaka, T.J., Michaelis, I., Rauberg, J. & Tøffner-Clausen, L., 2014. The CHAOS-4 geomagnetic field model, *Geophys. J. Int.*, **197**, 815–827.

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

Olson, P., Christensen, U.R. & Glatzmaier, G.A., 1999. Numerical modeling of the geodynamo: mechanisms of field generation and equilibration, *J. geophys. Res.*, **104**, 10 383–10 404.

Olson, P., Driscoll, P. & Amit, H., 2009. Dipole collapse and reversal precursors in a numerical dynamo, *Phys. Earth planet. Inter.*, **173**, 121–140.

Olson, P., Deguen, R., Hinnov, L. & Zhong, S., 2013. Controls on geomagnetic reversals and core evolution by mantle convection in the Phanerozoic, *Phys. Earth planet. Inter.*, **214**, 87–103.

Pais, M.A., Morozova, A.L. & Schaeffer, N., 2015. Variability modes in core flows inverted from geomagnetic field models, *Geophys. J. Int.*, **200**(1), 402–420.

Panovska, S., Korte, M., Finlay, C.C. & Constable, C.G., 2015. Limitations in paleomagnetism data and modelling techniques and their impact on Holocene geomagnetic field models, *Geophys. J. Int.*, **202**, 402–418.

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

Pozzo, M., Davies, C., Gubbins, D. & Alfè, D., 2012. Thermal and electrical conductivity of iron at Earth's core conditions, *Nature*, **485**, 355–358.
Suttie, N., Holme, R., Hill, M.J. & Shaw, J., 2011. Consistent treatment of errors in archaeointensity implies rapid decay of the dipole prior to 1840, *Earth planet. Sci. Lett.*, **304**(1, 13–21.

Takahashi, F., Matsushima, M. & Honkura, Y., 2007. A numerical study on magnetic polarity transition in an MHD dynamo model, *Earth Planets Space*, **59**, 665–673.

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

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

Thébault, E. & Gallet, Y., 2010. A bootstrap algorithm for deriving the archeomagnetic field intensity variation curve in the Middle East over the past 4 millennia BC, *Geophys. Res. Lett.*, **37**, L22303, doi:10.1029/2010GL044788.

Vonmoos, M., Beer, J. & Muscheler, R., 2006. Large variations in Holocene solar activity: constraints from $^{10}$Be in the Greenland Ice Core Project ice core, *J. geophys. Res.*, **111**, A10105, doi:10.1029/2005JA011500.

Wardinski, I. & Korte, M., 2008. The evolution of the core-surface flow over the last seven thousands years, *J. geophys. Res.*, **133**, doi:10.1029/2007JB005024.

Wicht, J. & Olson, P., 2004. A detailed study of the polarity reversal mechanism in a numerical dynamo model, *Geophys. Geochem. Geosyst.*, **5**, doi:10.1029/2003GC000602.

Zhang, P., Cohen, R. & Haule, K., 2015. Effects of electron correlations on transport properties of iron at Earth's core conditions, *Nature*, **517**, 605–607.

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

**Figure S1.** Time dependent longitude in degrees of reversed flux patches in the family of models A_FM from 1000 BC to

2000 AD. The same coloured diamonds are used for the position of a given RFP. The A_FM-0918, A_FM-0677 and A_FM-0318 models represent minimum, mean and maximum damping, respectively.

**Figure S2.** Time dependent co-latitude in degrees of the reversed flux patches in the family of models A_FM from 1000 BC to 2000 AD. The same coloured diamonds are used for the position of a given RFP. The A_FM-0918, A_FM-0677 and A_FM-0318 models represent minimum, mean and maximum damping, respectively.

**Figure S3.** As in Fig. S1 for the family of models ASD_FM. The ASD_FM-0683, ASD_FM-0377 and ASD_FM-0097 models represent minimum, mean and maximum damping, respectively.

**Figure S4.** As in Fig. S2 for the family of models ASD_FM. The ASD_FM-0683, ASD_FM-0377 and ASD_FM-0097 models represent minimum, mean and maximum damping, respectively.

**Figure S5.** As in Fig. S1 for normal flux patches for the ASDI_FM family of models.

**Figure S6.** As in Fig. S2 for normal flux patches for the ASDI_FM family of models. The dotted lines denote the tangent cylinder.

**Figure S7.** Histograms of reversed flux patches occurrences for the past 3 kyr in longitude (left) and co-latitude (right) for the A_FM family of models. Q is the number of reversed flux patches at a position range integrated over all snapshots. The histograms are constructed with bins of 20°.

**Figure S8.** As in Fig. S7 for the ASD_FM family of models.

**Table S1.** Data sources and treatment.

**Table S2.** Modelling strategies.

(http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggw248/-/DC1)

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.