# THE EVOLVING GEOMAGNETIC FIELD

EDITED BY : Greig A. Paterson, Christopher J. Davies and Ron Shaar PUBLISHED IN : Frontiers in Earth Science

#### Frontiers Copyright Statement

© Copyright 2007-2019 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88945-895-0 DOI 10.3389/978-2-88945-895-0

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

## THE EVOLVING GEOMAGNETIC FIELD

Topic Editors:

Greig A. Paterson, University of Liverpool, United Kingdom Christopher J. Davies, University of Leeds, United Kingdom Ron Shaar, The Hebrew University of Jerusalem, Israel

Earth's magnetic field has protected our planet for billions of years and provides key insights into the internal workings of our home planet. The geomagnetic field varies in distinctive fashions across a broad spectrum of timescales from milliseconds to millions of years. To understand these variations, Earth scientists utilize a diverse arsenal of tools from hi-tech satellites, such as the Swarm array, to archeological pottery and geological materials, through to advanced numerical simulations that harness the power of supercomputers. Armed with these tools we tackle problems related to the ancient magnetic field, how the geodynamo works and what this means for modern life. Despite being studied for more than 400 years, there are many unanswered questions about the geomagnetic field. This Research Topic on "The Evolving Geomagnetic Field" brings together these varied approaches to present our latest understanding of the workings of the geodynamo and the geomagnetic field across all timescales.

Citation: Paterson, G. A., Davies, C. J., Shaar, R., eds. (2019). The Evolving Geomagnetic Field. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-895-0

## Table of Contents


Peter E. Driscoll and Cian Wilson


Ron Shaar, Erez Hassul, Kate Raphael, Yael Ebert, Yael Segal, Ittai Eden, Yoav Vaknin, Shmuel Marco, Norbert R. Nowaczyk, Annick Chauvin and Amotz Agnon

*54 Archeomagnetic Intensity Spikes: Global or Regional Geomagnetic Field Features?*

Monika Korte and Catherine G. Constable


Andreas Nilsson, Neil Suttie and Mimi J. Hill

## Editorial: The Evolving Geomagnetic Field

#### Greig A. Paterson<sup>1</sup> \*, Christopher J. Davies <sup>2</sup> and Ron Shaar <sup>3</sup>

*<sup>1</sup> Department of Earth, Ocean and Ecological Sciences, University of Liverpool, Liverpool, United Kingdom, <sup>2</sup> School of Earth & Environment, University of Leeds, Leeds, United Kingdom, <sup>3</sup> The Institute of Earth Sciences, The Hebrew University of Jerusalem, Jerusalem, Israel*

Keywords: paleomagnetism, archeomagnetism, geodynamo modeling, geomagnetism, core dynamics, secular variation, global magnetic field models, time-average fields

**Editorial on the Research Topic**

#### **The Evolving Geomagnetic Field**

Earth's magnetic field, generated over 2,800 km below the surface in the liquid core, is a long-lived and complicated feature of our planet that requires a multi-stranded approach to quantify and comprehend. From its early beginning in the Archean to the modern day, the geomagnetic field exhibits striking spatial and temporal behavior across a vast range of scales, including sustained periods of atypical stability (superchrons), polarity reversals, rapid and spatially localized intensity "spikes," and impulsive changes in field strength and direction such as geomagnetic jerks.

Characterizing and understanding these complex variations requires a diverse approach that synthesizes experimental, observational, theoretical, and computational approaches. In this collection of papers that contribute to the topic "The Evolving Geomagnetic Field," we bring together studies on secular variation, geodynamo simulations, paleomagnetism, archeological magnetism, magnetic field modeling, and sedimentary magnetism to shed light on the latest advances in exploring our magnetic field.

Our collection opens with Muxworthy who explores the latitudinal dependence of the magnetic field strength (paleointensity) over the last 5 Myrs. Building on existing models (Tauxe and Kent, 2004), Muxworthy demonstrates that models where non-axial dipole terms time-average to zero are incapable of capturing the variability of the observed paleointensity data. He concludes that the current paleointensity data requires a constant axial quadrupole term of up to −10% of the axial dipole combined with an octupole term of −15% of the dipole. Interestingly, this octupole term has the opposite sign to the axial dipole, which contrasts with estimates from directional studies and suggests that time-average field models should incorporate full vector analysis. Driscoll and Wilson also consider long-timescale behavior by utilizing geodynamo simulations to explore time-averaging of the field to isolate a dominant axial dipole. Their analysis suggests that averaging over 20–120 kyr is required to obtain stable paleomagnetic poles, but for a frequently reversing field, longer periods may be needed. Obtaining reliable averages of the axial dipole moment, however, requires longer time periods and intensity estimates from less stable fields may introduce a latitudinal bias away from the true axial dipole moment, which may be corrected for if the reversal rate is known.

On shorter millennial timescales, archeomagnetic studies in China flourished in the 1980's and 90's, but went through a period of quiescence for more than a decade. Cai et al. review the rejuvenation of Chinese archeomagnetism over the past 5 years and present the state-of-the-art knowledge for how Earth's magnetic field has varied across China for the last ∼7000 years. From the most extensive collection to-date, Cai et al. construct a new series of archeomagnetic master curves for China, which will serve as an invaluable dating tool for Chinese archeologists.

Edited and reviewed by: *Kenneth Philip Kodama, Lehigh University, United States*

\*Correspondence: *Greig A. Paterson greig.paterson@liverpool.ac.uk*

#### Specialty section:

*This article was submitted to Geomagnetism and Paleomagnetism, a section of the journal Frontiers in Earth Science*

> Received: *06 March 2019* Accepted: *25 March 2019* Published: *10 April 2019*

#### Citation:

*Paterson GA, Davies CJ and Shaar R (2019) Editorial: The Evolving Geomagnetic Field. Front. Earth Sci. 7:77. doi: 10.3389/feart.2019.00077*

On similar timescales, but the other side of the world, archeomagnetism in the Middle East continues to produce valuable results as exemplified by Shaar et al., who compile the first catalog of archeomagnetic directional data from Israel, consisting of 76 directions. When pooled together with directional data from nearby regions, strong deviations from an axial dipole field are observed at 800–1000 BCE, lending further support to anomalous field behavior at this time, known as the Levantine Iron Age Anomaly (LIAA). The LIAA is also the focus of Korte and Constable, who investigate the two rapid spikes in field intensity during this time period using a new series of global magnetic field models. By preferentially up-weighting or downweighting data to construct their models, they conclude that these rapid field changes are the result of the growth and decay of quasistationary flux patches superposed on a stronger, but variable dipole. Requiring axial dipole variations no more that 60% higher than the present day, their analyses suggest that the intensity changes associated with the LIAA may be less anomalous than previously thought.

Yamamoto and Yamaoka turned to lava flows from the island of Hawaii to obtain absolute paleointensity over the last ∼24 kyr, extending the time range of archeomagnetic studies. The results of Yamamoto and Yamaoka were obtained from the Tsunakawa-Shaw paleointensity method, which is based on a specimen's coercivity spectrum. These new results are consistent with those from the IZZI method (Cromwell et al., 2015, 2018), which are based on a specimen's blocking temperature spectrum. This agreement between two distinct methods boosts our confidence in their overall reliability and is a great example of different approaches complimenting each other to yield robust records.

The final two contributions to this Research Topic use sedimentary sequences as detailed recorders of geomagnetic field

#### REFERENCES


behavior. Firstly, Lund brings together marine and lake sediment records from six regions of the world, spanning the Americas and the Philippines/Indonesia and dating back to ∼70 ka. Lund not only shows that the records are regionally consistent, but that the variability of the records is generally globally coherent with relative paleointensity having the highest consistency. Finally, Nilsson et al. focus on varved lake sediments from Sweden to develop and test a new Bayesian method to account for sedimentary lock-in depth (the depth to which sediments must be buried to "freeze-in" their magnetic recording). This is a potentially valuable new tool, not only to understand the lockin process of sedimentary records, but also for incorporating sedimentary records in long-term models of field behavior, such as those of Korte and Constable, which are part of the family of related models that have been used here by Nilsson et al. themselves, but also by Cai et al., Shaar et al., and Yamamoto and Yamaoka.

The breadth and inter-connectivity of the science presented here is the tip of the large iceberg that is the study of the geomagnetic field. We hope that the articles compiled in this Research Topic provide an enjoyable read that stimulates valuable discussion. These contributions have pushed forward our understanding of field behavior across a range of timescales, but also new questions and challenges. All of which are the seeds for new collaborations and renewed understanding of our magnetic field.

#### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Paleomagnetic Field, eds J. E. T.Channell, D. V. Kent, W. Lowrie, and J. G. Meert (Washington, DC: AGU), 101–115.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Paterson, Davies and Shaar. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Considerations for Latitudinal Time-Averaged-Field Palaeointensity Analysis of the Last Five Million Years

#### Adrian R. Muxworthy\*

*Natural Magnetism Group, Department of Earth Science and Engineering, Imperial College London, London, United Kingdom*

Central to palaeomagnetism and geophysics is the assumption that the time-averaged geomagnetic field is approximated by a geocentric-axial-dipole (GAD). In this paper, it is demonstrated through the use of a simple cap model that due to secular variation the time-averaged palaeointensity record will always have a smaller latitudinal dependency than a true GAD field. However, the simple cap model does not fully explain the behavior of the palaeointensity database (averaged over 0–5 Ma) especially at high-latitudes. To investigate this dependency I use a Giant Gaussian Processes (GGP) model to estimate the contribution of permanent non-dipole features and determine their statistical significance. It was found that an axial quadrupole term between −5 and −10% of the GAD field combined with octupole term ∼ −15% of the GAD field, best explained palaeointensity latitudinal behavior. In particular, the octupole term with a sign opposite to that of the GAD, is required to describe the palaeointensity behavior at high latitudes, i.e., >60◦ .

#### Edited by:

*Ron Shaar, Hebrew University of Jerusalem, Israel*

#### Reviewed by:

*Toni Veikkolainen, University of Helsinki, Finland Peter Aaron Selkin, University of Washington Tacoma, United States*

#### \*Correspondence:

*Adrian R. Muxworthy adrian.muxworthy@imperial.ac.uk*

#### Specialty section:

*This article was submitted to Geomagnetism and Paleomagnetism, a section of the journal Frontiers in Earth Science*

> Received: *16 August 2017* Accepted: *28 September 2017* Published: *12 October 2017*

#### Citation:

*Muxworthy AR (2017) Considerations for Latitudinal Time-Averaged-Field Palaeointensity Analysis of the Last Five Million Years. Front. Earth Sci. 5:79. doi: 10.3389/feart.2017.00079* Keywords: geocentric axial dipole hypothesis, time-averaged field, palaeointensity, palaeosecular variation, non-dipole field

#### INTRODUCTION

From geomagnetic observations and palaeomagnetic records, we know that the Earth's magnetic field is dominated by a dipole component, which is also dynamic and changing in terms of both its direction and intensity (Valet, 2003). Central to much palaeomagnetic research, e.g., palaeogeographic reconstructions, is the assumption that the time-averaged field (TAF) is a dipole field aligned parallel with the Earth's spin axis, the so-called geocentric-axial-dipole (GAD) hypothesis. Over very long periods of time, i.e., 200 million years, from palaeomagnetic directional data this hypothesis has been shown to hold (Merrill and McFadden, 2003); however, it has been known for some time (Wilson, 1970) that there are systematic departures from this simple model over shorter timescales. Direct field observations only encompass the last 400 years or so (Jackson et al., 2000), and reveal a strongly non-GAD field. If undetected, these non-GAD fields can have a major impact on palaeogeographical reconstructions creating spurious mismatches, for example, palaeomagnetic results from the timespan 0–5 Ma (Kono et al., 2000; Johnson et al., 2008) suggest a significant non-GAD contribution to the TAF of the order ∼5% (mostly quadrupole and octupole) that can correspond to ∼500 km reconstruction mismatches. Deviations from GAD are greatest at high-latitudes, but the exact nature and duration of these deviations is still debated (Constable, 2007). These high-latitudes deviations are greater in the ancient geomagnetic field intensity (palaeointensity) record than the directional record (Constable, 2007).

In the analysis of palaeointensity data as a function of latitude, it is common to plot palaeointensity data as a function of latitude for latitudinally binned "reliable" absolute palaeointensity data from the last 5 Myr using the PINT database (Biggin et al., 2009) available at http://earth.liv.ac.uk/pint; when binning the data sometimes the average of each bin is taken (e.g., Lawrence et al., 2009; Wang et al., 2015) and sometimes the median (e.g., Cromwell et al., 2015; Døssing et al., 2016). The definition of "reliable" varies between studies, but the same trend is observed: the measured intensity data deviates most from the GAD intensity model B <sup>D</sup> at high latitudes, where B <sup>D</sup> varies with the (magnetic) co-latitude θ by.

$$B^D = B\_{eq}^D (1 + 3 \cos^2 \theta)^{\frac{1}{2}} \tag{1}$$

where B D eq is the field intensity at the equator for an axial dipole field.

If the geomagnetic field was a steady-state GAD field, then a plot of the binned palaeointensity data from the last 5 Myr against latitude should lead to a dipole field, however, we know from the current-day geomagnetic field that the field it is not a GAD and it displays secular variation. However, by comparing the GAD intensity model and the binned palaeointensity data against latitude, it is intrinsically assumed that a TAF has the possibility of resembling a GAD field. In this paper, by considering the contribution of secular variation to the geomagnetic field, I show that it is simply incorrect; comparing binned palaeointensity data against a dipole field aligned along the rotation axis leads to an overestimation which is greatest at high latitudes, i.e., it is not possible for a TAF to be a GAD field due to secular variation. However, including dipole secular variation, does not explain the data trend found in the PINT database (Biggin et al., 2009); I further then investigate how permanent non-dipole structures might explain the latitudinal discrepancy.

#### A SIMPLE CAP MODEL FOR DIPOLAR SECULAR VARIATION

In studies that compare intensity vs. latitude (e.g., Lawrence et al., 2009; Cromwell et al., 2015; Wang et al., 2015; Døssing et al., 2016), it is the geographical latitude, not the palaeomagnetic latitude, that is used; this removes the circular logic of using GAD theory to calculate the palaeomagnetic latitude in a study of the robustness of the GAD hypothesis. Only palaeointensity data from the last 5 Myr are usually considered, because, as a first-order approximation, this removes the need to include any tectonic corrections to the sample locations' geographic latitudes. Once selected and latitudinally binned, the palaeointensity data are often directly compared to the intensity of the GAD field B D Equation (1).

However, the field experienced at a sample location S is not a GAD field, but as a first-order approximation a dipole field with a pole position P (**Figure 1**). The position of the current-day dipole field pole drifts, i.e., westward drift, 360◦ around the rotation axis and geographic pole G, with a tilt angle of ∼11◦ ; the tilt angle

also varies (Jackson et al., 2000; Constable, 2007). In determining the TAF intensity behavior of the palaeointensity database, we need to include secular variation of pole position P; it is simply incorrect to average all the data points from a given latitude as the field experienced at a given latitude depends on the magnetic colatitude, which varies with time and depends on a cosine function of the co-latitude Equation (1).

corresponding integration surface for Equation (3) in gray.

To demonstrate the effect of latitude on dipolar secular variation, consider a simple cap model of dipolar secular variation: The average great-circle length between P and S, i.e., g (**Figure 1**), which equals the magnetic co-latitude on a unitsphere. For a fixed sample location S, the arc-length g depends on both γ (westward drift) and tilt angle (arc-length p). γ varies between 0 and 360◦ , and the tilt-angle/arc-length p between pmin and pmax; if pmin = 0 the area is simply a cap. The average g for a point S can be determined using the Law of Cosines, which on a sphere with unit radius gives:

$$\frac{\frac{\mathcal{P}\_{\text{max}}}{\int \cos \overline{g}} - \frac{2\pi}{\int \sin \overline{\rho} \cos s + \sin \overline{\rho} \sin s \cos \overline{\rho} \, \sin \overline{\rho} \, d\overline{\rho}}{\int \frac{\mathcal{P}\_{\text{max}}}{\int \sin \overline{\rho} \sin \overline{\rho} \, d\overline{\rho}} \, d\overline{\rho}}$$

$$= \frac{\cos s}{4} \frac{(\cos 2\overline{\rho}\_{\text{min}} - \cos 2\overline{\rho}\_{\text{max}})}{(\cos \overline{\rho}\_{\text{min}} - \cos \overline{\rho}\_{\text{max}})}\tag{2}$$

where cos g is the average cosg for the area under consideration.

To study the TAF, the palaeointensity data are binned from different time periods over intervals of typically 10–15◦ latitude (Wang et al., 2015; Døssing et al., 2016). The average of cos g smax smin over a sample latitude window smin to smax is simply given by:

$$\frac{\varepsilon\_{\text{max}}}{\cos \mathcal{G}\_{s\_{\text{min}}}^{s\_{\text{max}}}} = \frac{\int\_{s=s\_{\text{min}}}^{s\_{\text{max}}} \left( \frac{\cos \mathcal{D}p\_{\text{min}} - \cos \mathcal{D}p\_{\text{max}}}{(\cos p\_{\text{min}} - \cos p\_{\text{max}})} \right) \sin s ds}{\int\_{s\_{\text{max}}}^{s\_{\text{max}}} \sin s ds}$$

$$= \frac{1}{16} \frac{(\cos 2p\_{\text{min}} - \cos 2p\_{\text{max}})}{(\cos p\_{\text{min}} - \cos p\_{\text{max}})} \frac{(\cos 2s\_{\text{min}} - \cos 2s\_{\text{max}})}{(\cos s\_{\text{min}} - \cos s\_{\text{max}})} \tag{3}$$

To calculate the time-averaged intensity B TAF as a function of latitude for a dipole field that displays secular variation as recorded by latitudinally binned palaeointensity data, simply put Equation (3) into Equation (1) to give:

$$B^{TAF} = B\_{eq}^{TAF} \left( 1 + 3 \cos^2 \overline{g}\_{s\_{\rm min}}^{s\_{\rm max}} \right)^{\frac{1}{2}} \tag{4}$$

where g smax smin is the average of g determined from Equation (3) for a site-latitude window between smin to smax, and B TAF eq is the field intensity at the equator for the TAF.

Using a latitudinal bin size of 10◦ , I have calculated Equation (4) for three model scenarios, and compared these to the simple dipole model Equation (1) in **Figure 2** for B TAF eq <sup>=</sup> <sup>B</sup> D eq <sup>=</sup> <sup>25</sup> <sup>µ</sup>T. The three models are: (1) <sup>p</sup>max <sup>=</sup> <sup>45</sup>◦ and pmin = 0 ◦ , (2) pmax = 20◦ and pmin = 0 ◦ , (3) <sup>p</sup>max <sup>=</sup> <sup>45</sup>◦ and <sup>p</sup>min <sup>=</sup> <sup>20</sup>◦ and (3) <sup>p</sup>max <sup>=</sup> <sup>15</sup>◦ and pmin = 5 ◦ . A maximum tilt angle of 45◦ was chosen as this is commonly argued as the latitudinal limit of a non-transitional pole position (Tauxe, 2010); however, the exact cut-off is not significant.

All three models have a weaker latitudinal dependency than the dipole model of Equation (1) (**Figure 2**), that is, the effect of including dipolar secular variation into the TAF intensity model is to reduce the expected latitudinal dependence of the observed field, i.e., a TAF palaeointensity record cannot have the same latitudinal dependency as a GAD field. For a constant equatorial field the deviation between the dipole model and the TAF model is enhanced at high latitudes. This is expected: consider the case where the sampling locality is approximately at the geographic North Pole, i.e., smax is slightly larger than smin, but both are close to 0◦ , and pmin >> 0 ◦ , then the value of g smax smin used in Equation (4) is >0 ◦ reducing B TAF at high latitudes; this effect is less pronounced for equatorial sampling localities. Generally the value of smax was found to be more important than smin, and the effect of bin size is relatively unimportant for bin sizes between 5 ◦ and 15◦ .

#### A GIANT GAUSSIAN PROCESSES (GGP) MODEL APPROACH TO SECULAR VARIATION

While the simple cap dipolar model demonstrates the effect of secular variation on binned palaeointensity data, it does not

FIGURE 2 | *B* latitude using a latitudinal bin size of 10◦ : 1) *<sup>p</sup>*max <sup>=</sup> <sup>45</sup>◦ , *p*min = 0 ◦ , 2) *p*max = <sup>20</sup>◦ , *p*min = 0 ◦ and 3) *<sup>p</sup>*max <sup>=</sup> <sup>45</sup>◦ , *<sup>p</sup>*min <sup>=</sup> <sup>20</sup>◦ . Additionally the standard dipole model (Equation1) is also calculated. In all four models the equatorial field was set to 25 µT, i.e., *B TAF eq* <sup>=</sup> *<sup>B</sup> D eq* <sup>=</sup> <sup>25</sup> <sup>µ</sup>T.

attempt to describe secular variation, as it is flat distribution. There is a class of models based on Giant Gaussian Processes (GGP) (Constable and Parker, 1988; Kono and Hiroi, 1996; Quidelleur and Courtillot, 1996; Constable and Johnson, 1999; Tauxe and Kent, 2004; Shcherbakov et al., 2014), that express secular variation by varying the Gauss coefficients g m l , h m l ; the gauss coefficients have a mean and standard deviations (σ<sup>l</sup> ) that are a function of degree l. Generally the mean g 0 1 (axial dipole) and g 0 2 (axial quadrupole) are non-zero, and the higher order g 0 m are zero, though occasionally non-zero mean g 0 3 (axial octupole) terms are considered (Tauxe and Kent, 2004) (**Table 1**).

The variation is partially accommodated by the α parameter, and is given by

$$
\sigma\_l^2 = \frac{\langle \mathfrak{c}/a \rangle^{2l} \alpha^2}{(l+1)(2l+1)} \tag{5}
$$

where c/a is the ratio of the core radius to that of the Earth's, and α a fitted parameter. Through time the models have developed, with increasing complexity used to describe the variation in the individual Gauss coefficients (**Table 1**). The early model of Constable and Parker (1988), did not accommodate the observed latitudinal variation in the scatter of the virtual geomagnetic pole (VGP) data. Later models did this by weighting some of the variance parameters (Quidelleur and Courtillot, 1996; Constable and Johnson, 1999; Tauxe and Kent, 2004) (**Table 1**).

Constable and Johnson (1999) compared their GGP model to the then palaeointensity record, and found a reasonable correlation. However, in the last 20 years the palaeointensity database has changed significantly, with the inclusion of far more high-latitude data points. Additionally, the selection criteria for drawing samples from the database are far more stringent now than it was in the late 1990s (e.g., Lawrence et al., 2009; Cromwell et al., 2015; Wang et al., 2015; Døssing et al., 2016). As we will



*The "best-fit" parameters from the model in this study were determined using the results shown in* Figure 5 *for g*<sup>0</sup> <sup>2</sup> <sup>=</sup> <sup>−</sup>*0.1 g*<sup>0</sup> 1 *and g*<sup>0</sup> <sup>3</sup> = −*0.15 g*<sup>0</sup> 1 *.*

see in the next section the GGP model of Constable and Johnson (1999) no longer explains the current PINT database.

#### COMPARISON OF CURRENT GGP MODELS WITH THE LATITUDINAL PALAEOINTENSITY DATA RECORD

To compare GGP models with the current palaeointensity data record, I selected Thellier data with pTRM checks from the PINT database (Biggin et al., 2009) from non-transitional (or unknown polarity) samples ≤5 Ma in age with at least three intensity estimates per unit, i.e., N ≥ 3, with a unit palaeointensity estimate standard deviation of <50%, as in previous studies (Døssing et al., 2016). This yields 496 field estimates. I also combined the northern and southern hemisphere data, because the southern hemisphere data is very limited both spatially and temporally for the last 5 Myr. For example, for a 10◦ bin size, several of the southern hemisphere bins have no data, and others have only one dataset covering narrow time periods, e.g., data from Tristan da Cunha (Shah et al., 2016) is the single data set in the −30 to −40◦ bin and covers a time period of <50 kyr. After folding and binning at 10◦ steps, the 80–90◦ bin has no data, but the next least populated bin, the 0–10◦ bin, has eight points.

To determine a value for each bin, Døssing et al. (2016) calculated the median of the data in each bin. Whilst this approach makes less assumptions, if we compare the bin with the most data points (bin 10–20◦ , 226 points,), we see that at 95% confidence it can be described as Gaussian (**Figure 3**). I therefore calculate the mean for each bin, rather than the median (**Figures 4**, **5**). I also determine a 95% confidence limit for each bin. The limits are effected by the number of points in each bin; the 10–20◦ bin has the narrowest confidence interval (**Figure 4**). The mean data are seen to increase with latitude until the 50– 60◦ bin, then decrease; however, the site data are very scattered (**Figure 3**). The simple dipole model Equation (1) increases more rapidly than the data (**Figure 4**), and does not describe the data.

FIGURE 3 | Histogram of bin count vs. intensity for the palaeointensity data from the PINT palaeointensity database (Biggin et al., 2009, http://earth.liv.ac. uk/pint), for the 10–20◦ bin. The palaeointensity data were selected from the PINT database from sites ≤5 Ma, for which there were at least three intensity estimates. The northern and southern hemisphere data are combined. The data is not distinguishable from a Gaussian distribution at 95% confidence.

Nor do the simple cap models (**Figure 2**) describe the data trend (**Figure 4**).

I next compare the CJ98 and TK03 GGP models (**Table 1**) with the selected PINT data, and determine 95% confidence intervals using the number of PINT data in each bin to calculate the models' confidence limits for comparison with the data (**Figure 4**). It is visually clear that neither of these two models follows the binned palaeointensity latitudinal trend. CJ98 and TK03 both increase too rapidly with latitude, plus TK03 is too low at the equator. Adjusting the g 0 1 term for TK03 model does not resolve the latitudinal variation inconsistency. To statistically estimate if the model and PINT data distributions are distinct, I calculated two parameters: (1) the Kolmogorov-Smirnov (KS) distribution probability (Press et al., 2007), and (2) the Bayes error (Braga Neto and Dougherty, 2015). The KS distribution probability estimates the likelihood that the data distributions for each bin are drawn from the same distribution, i.e., low probability values suggest that the distributions are significantly different (Press et al., 2007). The KS distribution probability examines only the shape of the distributions at each bin, but is invariant to absolute values. To quantify the absolute overlap between two normal distributions, I calculate the upper bound on the Bayes error (Braga Neto and Dougherty, 2015), to estimate the maximum misclassification error between two normal distributions with means µ<sup>0</sup> and µ<sup>1</sup> and standard deviations σ<sup>0</sup> and σ1. This is done by finding the value of γ, where 0 < γ < 1, which minimizes cγ.

$$c\_{\mathcal{V}} = \min\_{ax(0,1)} \left[ \frac{\sigma\_0^{1-\mathcal{V}} \sigma\_1^{\mathcal{V}}}{\sqrt{(1-\mathcal{V})\sigma\_0^2 + \mathcal{V}\sigma\_1^2}} \exp\left(\frac{(\mathcal{V}-1)\mathcal{Y}(\mu\_0 - \mu\_1)^2}{2(\mathcal{V}\sigma\_1^2 - (\mathcal{V}-1)\sigma\_0^2)}\right) \right] \tag{6}$$

FIGURE 4 | GGP model intensity predictions vs. latitude for: (A) CJ98 model, and (B) TK03 model. Also depicted are palaeointensity data (green dots) from the PINT palaeointensity database (Biggin et al., 2009, http://earth.liv.ac.uk/pint), which has been binned at 10◦ degree intervals and the mean determined. The palaeointensity data were selected from the PINT database from sites ≤5 Ma, for which there were at least three intensity estimates. The northern and southern hemisphere data are combined; there is no data <sup>&</sup>gt;80◦ . 95% confidence intervals are plotted for the PINT data (shaded green) and the model data (black lines). Additionally the standard dipole model (Equation1) is also calculated, with *B D eq* <sup>=</sup> <sup>25</sup> <sup>µ</sup>T (gray line), in (A) a cap model as shown in Figure 3 (*p*max <sup>=</sup> <sup>45</sup>◦ , *<sup>p</sup>*min <sup>=</sup> <sup>20</sup>◦ ). In the lower section of the figures the Kolmogorov-Smirnov probability and the Bayes error *E* are plotted as a function of latitude.

If we assume a priori that the probabilities that the observation came from each of the distributions is equal, then the upper bound on the probability of a misclassification error is E ≤ cγ/2. If E = 0, there is no probability of a misclassification, i.e., the two distributions are completely distinct from each other, and if E = 0.5 there is a 50:50 chance that misclassification will take place, therefore the two distributions must be identical. However, the probabilistic nature of the analysis can just be taken as an intuitive measure of the similarity between two distributions.

For CJ98 and TK03, both the KS distribution probability and the Bayes E are generally low, indicating that the models do not accurately describe the data. For CJ98, Bayes E is >0.45 for the 0– 10◦ bin, this is partially an artifact of the two mean values being similar.

#### PERMANENT NON-DIPOLE FEATURES AND THE PALAEOINTENSITY DATA RECORD

It is clear both by inspection and the statistical tests that the CJ98 and TK03 models to not accommodate the paleointensity data (**Figure 4**). To explore this mismatch, I examined the GGP models, starting with the TK03 model of Tauxe and Kent (2004). There are many parameters (**Table 1**) that can potentially change the shape of the geomagnetic field (Constable and Johnson, 1999; Merrill and McFadden, 2003), however, whilst changes to the standard deviation parameters can drastically increase the secular variation, when averaged over long periods of time, variations in these parameters contribute more to the scatter rather than the mean trend, which still resembles a dipole field, albeit reduced, as demonstrated in the simple cap model (**Figure 2**).

To change the model trends to resemble something closer to the PINT database, it is necessary to add permanent non-dipole terms, in particular g 0 2 (axial quadrupole) and g 0 3 (axial octupole). In TK03 (**Table 1**) g 0 2 is set to zero, however, in most other GGP models (Constable and Parker, 1988; Quidelleur and Courtillot, 1996; Constable and Johnson, 1999; Tauxe, 2005; Shcherbakov et al., 2014), this term is non-zero and of the same sign as g 0 1 . The literature provides some indications of what values of g 0 2 and g 0 3 are considered "reasonable": Using palaeomagnetic data for the last 5 Myrs, various studies have inverted the palaeomagnetic data to determine the time-averaged Gauss coefficients (e.g., Johnson and Constable, 1997; Carlut and Courtillot, 1998; Kono et al., 2000; Johnson et al., 2008). Of these studies, only that of Kono et al. (2000) attempted to fully invert the palaeointensity data record. They found that g 0 <sup>2</sup> ≈ −0.06 <sup>g</sup> 0 1 and g 0 <sup>3</sup> <sup>≈</sup> 0.06 <sup>g</sup> 0 1 , which is in contrast to the palaeodirectional data only inversions that generally predict that g 0 2 and g 0 3 have the same sign as g 0 1 ; the exception being some of the inversions for single polarities (Johnson and Constable, 1997; Johnson et al., 2008), where g 0 <sup>3</sup> was found to be of opposite sign to g 0 1 . Using geodynamo simulations, Veikkolainen et al. (2017) examined the entire PINT database , in an attempt to quantify the contribution of g 0 2 and g 0 3 to the TAF through examination of palaeointensity frequency distribution. They found that the paleointensity record is best described by periods of different g 0 1 intensities, however, each period had a ±10% octupole field contribution, with a minimal quadrupole contribution.

Given the limited number of binned data (8 points), the aim of the exploration was to reproduce the same trends as observed in the data and not to determine an exact fit. Therefore, in the following calculations, I varied only the g 0 1 , g 0 2 , and g 0 3 terms within the TK03 model. I normalized the TK03 model with

FIGURE 5 | GGP model intensity predictions vs. latitude for: (A) *g* 0 2 ≈ −0.05 *g* 0 1 (−5%) and *g* 0 3 ≈ 0 (0%), (B) *g* 0 2 ≈ −0.05 *g* 0 1 and *g* 0 3 ≈ −0.05 *g* 0 1 , (C) *g* 0 2 ≈ 0 and *g* 0 3 ≈ −0.1 *g* 0 1 , (D) *g* 0 2 ≈ 0.05 *g* 0 1 and *g* 0 3 ≈ −0.1 *g* 0 1 , (E) *g* 0 2 ≈ −0.05 *g* 0 1 and *g* 0 3 ≈ −0.15, and (F) *g* 0 2 ≈ −0.1 *g* 0 1 and *g* 0 3 ≈ −0.15 *g* 0 1 . Included are the binned palaeointensity data (green dots) from the PINT palaeointensity database. The northern and southern hemisphere data are combined; there is no data <sup>&</sup>gt;80◦ , 95% confidence intervals are plotted for the PINT data (shaded green) and the model data (black lines). Additionally, the standard dipole model (Equation1) is also calculated, with *B D eq* <sup>=</sup> <sup>25</sup> <sup>µ</sup>T (gray line). In the lower section of the figures the Kolmogorov-Smirnov probability and the Bayes error *<sup>E</sup>* are plotted as a function of latitude.

respect to the most accurate PINT database bin (10–20◦ bin), i.e., g 0 <sup>1</sup>was calculated for each model. This normalizing clearly affects the Bayes error, however, it is consistent across the models for comparison.

Adding g 0 2 and g 0 3 terms of the same sign as g 0 1 simply increases the latitudinal variation, and makes the fit between the models and the data worse than when g 0 2 and g 0 3 are both zero. Therefore, to match the PINT database trends for the last 5 Myrs, permanent g 0 2 and g 0 3 terms of the opposite sign to g 0 1 are required. Various combinations of g 0 2 and g 0 3 are plotted in **Figure 5**; g 0 2 and g 0 3 are expressed as percentages of g 0 1 . The effect of introducing g 0 2 and g 0 3 of the opposite sign to g 0 1 is to decrease the rate of latitudinal increase, producing a 'flatter' trend. It was found that to represent the data trend at high latitudes under the constraint g 0 2 , g 0 <sup>3</sup> < 0.2 g 0 1 , both g 0 2 and g 0 3 have to be of the opposite sign to g 0 1 .

In **Figure 5A**, I consider a 5% quadrupole field with no octupole field, as suggested by several palaeodirectional studies (e.g., Carlut and Courtillot, 1998). I plot only the case where g 0 <sup>2</sup> <sup>=</sup> <sup>−</sup>0.05 <sup>g</sup> 0 1 , as this trend matches the data better than for g 0 <sup>2</sup> = +0.05 <sup>g</sup> 0 1 , however, its latitudinal dependency is too high even for the negative contribution. Adding a g 0 <sup>2</sup> = +0.05 <sup>g</sup> 0 1 term to this figure, i.e., approximately the same configuration as suggested by Kono et al. (2000), makes the mismatch between model and data worse than in **Figure 5A**. The effect of adding a g 0 <sup>3</sup> = −0.05 <sup>g</sup> 0 1 term to a g 0 <sup>2</sup> = −0.05 <sup>g</sup> 0 1 , sharply decreases the field intensity at high-latitudes (**Figure 5B**). The behavior for mixed quadrupole/octupole −5% contribution (**Figure 5B**) is similar to that of a pure octupole term of −0.1 g 0 1 (**Figure 5C**). This latter configuration is the same as that suggested Veikkolainen et al. (2017) who examined the entire PINT database. Attempting to combine quadrupole and octupole terms of opposite signs does not improve the data fit (**Figure 5D**). The models that best describe the data are shown in **Figures 5E,F**; that is, g 0 <sup>3</sup> = −0.15 <sup>g</sup> 0 1 and g 0 <sup>2</sup> = −0.05 and <sup>−</sup>0.1 <sup>g</sup> 0 1 . Even within these models there is still a large degree of mismatch at mid-latitude. The g 0 1 term determined by fitting the models to the 10–20◦ bin, was consistently between 33.0 and 33.2 µT, a little larger than the values commonly quoted, e.g., 30 µT (Constable and Johnson, 1999; Shcherbakov et al., 2014).

The g 0 2 and g 0 3 contributions are higher than that found from the directional palaeomagnetic data (e.g., Carlut and Courtillot, 1998), however, this is not inconsistent. The large g 0 <sup>3</sup> was needed to accommodate the high-latitude palaeointensity data behavior; it is well-documented that palaeointensity data is more sensitive to high-latitude geomagnetic field behavior than declination and inclination data (Constable, 2007). Therefore, it appears that the paleointensity database is highlighting high-latitudinal behavior not accessible to the declination and inclination

#### REFERENCES

Biggin, A. J., Strik, G. H. M. A., and Langereis, C. G. (2009). The intensity of the geomagnetic field in the late-Archaean: new measurements and an analysis of the updated IAGA palaeointensity database. Earth Planets Space 61, 9–22. doi: 10.1186/BF03352881

data. Veikkolainen et al. (2017) who also studied only the palaeointensity record, also concluded that large permanent octupole features are required to explain database trends. In contrast to this study, Kono et al. (2000) found g 0 <sup>2</sup> ≈ −0.06<sup>g</sup> 0 1 and g 0 <sup>3</sup> <sup>≈</sup> 0.06<sup>g</sup> 0 1 for the TAF (0–5 Ma), i.e., g 0 3 is of the same sign as g 0 1 ; however, as stated above the palaeointensity database has improved significantly since the late 1990s.

The KS probability was consistently low for bin 30–40◦ ; with 65 points it was the second most populated bin, and is normal at 95% confidence (**Figure 5**). That is, the GGP distribution prediction for this bin does not follow the observed behavior well-regardless of g 0 2 and g 0 3 contributions. Given the problematic nature of palaeointensity collection, it is suggested more data at this mid-latitude are required.

#### CONCLUSIONS

This paper has shown through a simple cap model of secular variation that a TAF intensity field cannot theoretically display the same behavior as a pure GAD, and will always have a trend that is less dependent on latitude (**Figure 2**). However, the simple cap model does not explain the latitudinal behavior seen in the PINT database (**Figure 4**). Through the use of GGP models it was found that the TAF (0–5 Ma) has likely permanent non-dipole features, in particular an axial quadrupole term ≈ −0.1g 0 1 and octupole term ≈ −0.15g 0 1 . The g 0 1 term was determined to be ∼33 µT. These permanent non-dipole features are larger than reported in other studies, however, these other studies either examined only the directional palaeomagnetic data (e.g., Carlut and Courtillot, 1998) or considered much early versions of the palaeointensity database where the high-latitude palaeointensity data is less constrained than today (Constable and Johnson, 1999). However, while the PINT database has far more data now than 20 years ago, the data at high-latitudes, even after folding, are very limited and subject to incomplete temporal sampling; more high-latitude data are required. Data within the 30–40◦ bin also displays unusual behavior, which would benefit from more data.

#### AUTHOR CONTRIBUTIONS

The author confirms being the sole contributor of this work and approved it for publication.

#### ACKNOWLEDGMENTS

AM would like to thank Dr. David Heslop for fruitful discussions regarding Bayes error analysis.

Braga-Neto, U. M., and Dougherty, E. R. (eds.). (2015). "Performance analysis," in Error Estimation for Pattern Recognition (Hoboken, NJ: John Wiley & Sons, Inc.). doi: 10.1002/9781119079507.ch3

Carlut, J., and Courtillot, V. (1998). How complex is the time-averaged geomagnetic field over the past 5 Myr? Geophys. J. Int. 134, 527–544.


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Muxworthy. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Paleomagnetic Biases Inferred From Numerical Dynamos and the Search for Geodynamo Evolution

Peter E. Driscoll\* and Cian Wilson

Department of Terrestrial Magnetism, Carnegie Institution for Science, Washington, DC, United States

The paleomagnetic record is central to our understanding of the history of the Earth. The orientation and intensity of magnetic minerals preserved in ancient rocks indicate the geodynamo has been alive since at least the Archean and possibly the Hadean. A paleomagnetic signature of the initial solidification of the inner core, arguably the singular most important event in core history, however, has remained elusive. In pursuit of this signature we investigate the assumption that the field is a geocentric axial dipole (GAD) over long time scales. We study a suite of numerical dynamo simulations from a paleomagnetic perspective to explore how long the field should be time-averaged to obtain stable paleomagnetic pole directions and intensities. We find that running averages over 20 − 40 kyr are needed to obtain stable paleomagnetic poles with <sup>α</sup><sup>95</sup> <sup>&</sup>lt; <sup>10</sup>◦ , and over 40 − 120 kyr for α<sup>95</sup> < 5 ◦ , depending on the variability of the field. We find that models with higher heat flux and more frequent polarity reversals require longer time averages, and that obtaining stable intensities requires longer time averaging than obtaining stable directions. Running averages of local field intensity and inclination produce reliable estimates of the underlying dipole moment when reversal frequency is low. However, when heat flux and reversal frequency are increased we find that local observations tend to underestimate virtual dipole moment (VDM) by up to 50% and overestimate virtual axial dipole moment (VADM) by up to 150%. A latitudinal dependence is found where VDM underestimates the true dipole moment more at low latitudes, while VADM overestimates the true axial dipole moment more at high latitudes. The cause for these observed intensity biases appears to be a contamination of the time averaged field by non-GAD terms, which grows with reversal frequency. We derive a scaling law connecting reversal frequency and site paleolatitude to paleointensity bias (ratio of observed to the true value). Finally we apply this adjustment to the PINT paleointensity record. These biases produce little change to the overall trend of a relatively flat but scattered intensity over the last 3.5 Ga. A more careful intensity adjustment applied during periods when the reversal frequency is known could reveal previously obscured features in the paleointensity record.

Keywords: geodynamo, paleointensity, earth evolution, paleomagnetism, core dynamics

Edited by:

Christopher Davies, University of Leeds, United Kingdom

#### Reviewed by:

Andrew John Biggin, University of Liverpool, United Kingdom Florian Lhuillier, Ludwig-Maximilians-Universität München, Germany

> \*Correspondence: Peter E. Driscoll pdriscoll@ciw.edu

#### Specialty section:

This article was submitted to Geomagnetism and Paleomagnetism, a section of the journal Frontiers in Earth Science

> Received: 12 March 2018 Accepted: 25 July 2018 Published: 21 August 2018

#### Citation:

Driscoll PE and Wilson C (2018) Paleomagnetic Biases Inferred From Numerical Dynamos and the Search for Geodynamo Evolution. Front. Earth Sci. 6:113. doi: 10.3389/feart.2018.00113

### 1. INTRODUCTION

Our knowledge of the history of Earth's magnetic field derives from paleomagnetic signals preserved in rocks. Many applications of paleomagnetism rely on an assumption that only the geocentric axial dipole (GAD) component remains after averaging the complex time-variable magnetic field over a sufficient amount of time, typically assumed to be around 10–20 kyr (Merrill and McFadden, 2003). This GAD field assumption has been extremely rewarding, for example in obtaining paleointensities (e.g., Biggin et al., 2009; Tauxe and Yamazaki, 2015), paleodirections and tectonic reconstructions (e.g., Torsvik et al., 2012; Raub et al., 2015), and even paleoclimate studies that rely on paleomagnetically derived paleolatitudes (Evans et al., 2000; Williams and Schmidt, 2004). Tests of the GAD field assumption have generally found support for its validity (e.g., Johnson et al., 1995; Acton et al., 1996; Meert et al., 2003; McElhinny, 2004; Evans, 2006; Swanson-Hysell et al., 2009; Panzik and Evans, 2014; Veikkolainen et al., 2014, 2017; Johnson and McFadden, 2015), although some have proposed long-term deviations from GAD in the Precambrian (Kent and Smethurst, 1998; Abrajevitch and Van der Voo, 2010).

Theoretically a GAD field is predicted over long time averaging because the equations governing dynamo action are symmetric about the transformation <sup>B</sup><sup>E</sup> → −BE, symmetric about the equator, and symmetric in rotation about the polar axis (Gubbins and Zhang, 1993). These symmetries imply that if random samples are time averaged long enough only the axial dipole term should retain a non-zero amplitude. However, the length of time required to average out all non-GAD terms is typically assumed rather than measured, and could itself vary significantly, especially during periods when the field is highly variable or frequently reversing (Merrill and McFadden, 2003).

The remarkable progress in reconstructing the motions of the continents has been fueled in large part by the success of the GAD assumption. However, some anomalous data are not explained by a simple GAD field. In particular, a growing number of anomalous directions in the Neoproterozoic (e.g., Maloof et al., 2006; Abrajevitch and Van der Voo, 2010; McCausland et al., 2011; Swanson-Hysell et al., 2012; Halls et al., 2015; Klein et al., 2015; Levashova et al., 2015; Bazhenov et al., 2016) have been variously interpreted as caused by extremely rapid plate motions, significant true polar wander, long-term non-GAD magnetic field components, or some mixture of these.

An additional puzzle that has garnered recent attention is the surprising lack of an obvious paleomagnetic signature of inner core nucleation (ICN) in the paleointensity record. Biggin et al. (2015) found a paleointensity peak between 1.0 and 1.5 Ga that could be a signature of ICN. Recently the primary signature of some of the underlying data has been questioned and revised (e.g., Smirnov et al., 2016; Sprain et al., 2018). Simultaneously, recent upward revisions of the thermal conductivity of iron have decreased estimates of the age of the inner core to Neoproterozoic time (e.g., Driscoll and Bercovici, 2014; Davies, 2015; Nimmo, 2015), although significant uncertainty remains.

There are at least three possible reasons for the lack of a clear paleomagnetic signature of ICN: (1) the paleomagnetic signature of ICN is too small or old to be preserved, (2) the paleointensity record is too sparse, or (3) the signature is obscured by non-GAD fields. Regarding the latter possibility, it has recently been proposed that prior to inner core nucleation around 600 Ma a non-GAD field may have been persistent in the Neoproterozoic as a consequence of the geodynamo being powered only by weak thermal convection at the time (Driscoll, 2016; Landeau et al., 2017). Unfortunately the paleointensity record around this time is sparse, possibly due to a lack of wide spread magmatism, a lack of preservation, inability to recover primary remanence, or low quality criteria. In any case, there is an impetus to investigate how paleomagnetic recordings are affected by a range of dynamo regimes and field morphologies from both empirical and theoretical grounds. Obtaining new high quality data, developing new analysis techniques of old data, and investigating synthetic data from numerical dynamo models all provide a way forward. In this paper we focus on the latter approach.

Several previous studies have generated synthetic observations from numerical dynamos for different purposes. Wicht (2005) found the that the observed length of reversal durations can change by an order of magnitude as a function of observed site latitude. A statistical analysis of several numerical dynamos by McMillan et al. (2001) found significant variation in field components when averaged over 100 kyr, and that a minimum of 10 dipole decays times were required to obtain stable estimates of the dipole field. Similarly Davies and Constable (2014) found that averaging over several hundred thousand years was required to obtain stable dipole field estimates, and longer averaging is needed for more turbulent (higher Rm) dynamo models. Lhuillier and Gilder (2013) found that roughly one million years was required to achieve stationary intensities and directions, and that these quantities correlate with stable chron duration.

In this paper we systematically explore how long a time average is required to obtain a GAD field from a range of dynamo regimes that span stable dipolar to reversing nondipolar. We generate local synthetic (or "virtual") geomagnetic observations from these models to investigate possible intrinsic biases generated by the core magnetic field itself; i.e., not caused by rock magnetic affects, alterations, or external forcings. In particular we aim to identify whether the dynamics of the core can produce predictable biases in the time averaged paleomagnetic field, in terms of both paleomagnetic directions and paleointensities, and whether such biases can be identified and removed from paleomagnetic data to reveal previously obscured features.

The paper is organized into two parts. The first part is an analysis of numerical dynamo models by synthetic observations. In section 2 we introduce the numerical dynamos and synthetic analysis methods, followed by results in section 3. The second part is an application to the paleointensity record. In section 4 we review the paleointensity record, accounting for several recently identified issues with certain paleointensity estimates, and adjust the record according to the results of the dynamo models. Finally, implications and conclusions are in section 5.

TABLE 1 | Dynamo model properties: basal heat flux dT/dri , time averaged volumetric kinetic energy (KE) and magnetic energy (ME), length of run 1t in kyr, number of dipole equator crossings Neq, number of polarity reversals Nrev, and critical running average length <sup>τ</sup>crit when <sup>α</sup><sup>95</sup> <sup>&</sup>lt; <sup>10</sup>◦ .


#### 2. NUMERICAL DYNAMOS AND ANALYSIS METHODS

We produce a suite of numerical dynamo models spanning a range of behavior from stable dipolar to reversing nondipolar. We analyze these models using synthetic "observations" analogous to paleomagnetic data (i.e., from locations on Earth's surface). Possible biases and correlations between "known" and "observed" quantities will be quantified as a function of the length of time averaging, reversal rate, and site co-latitude.

#### 2.1. Dynamo Model Setup

The dynamo models are computed using the Rayleigh dynamo code (Featherstone and Hindman, 2016; Matsui et al., 2016). All models share the following control parameters: E = 10−<sup>3</sup> , Ra = 10<sup>6</sup> , Pr = 1, Pm = 10, insulating magnetic boundary conditions, no-slip velocity conditions, inner-outer core radius ratio of 0.35, and fixed temperature gradient at both boundaries (see **Table 1**). These relatively high Ekman number simulations produce Earth-like large scale magnetic features (see below) and polarity reversals that resemble geomagnetic observations, and are numerically cheap so they can be run extremely long times to produce low frequency statistics. The inner boundary temperature gradient (i.e., basal heat flux) dT/dr<sup>i</sup> , fixed in time in each model, spans a range of 0.8 − 12 (in non-dimensional units). In the dynamo code the Rayleigh number is defined in terms of the temperature drop across the shell 1T for isothermal boundary conditions. For fixed flux Ra is defined by the specified lower boundary temperature gradient by 1T = D dT dri , where D = 1 is the dimensionless shell thickness. These basal heat fluxes produce a range of behavior from stable dipolar dynamos on the lower end to unstable, reversing non-dipolar dynamos at the high end (**Table 1**). The temperature gradient at the outer boundary is set to balance the heat flow at the base so that energy is conserved and there is no internal sink or source.

Time is scaled from thermal diffusion times (implemented in the code) to years by multiplying by a factor τdip/(Pm(ro/πD) 2 ), where r<sup>o</sup> = 1.5384 is dimensionless outer core radius, and τdip = 50 kyr is the assumed magnetic dipole decay time of the core. We adopt the same definition of a polarity reversal proposed by Driscoll and Olson (2009): that the dipole co-latitude θdip spend at least 20 kyr in a stable polarity before and after a reversal. We note that the time scaling is not unique and adopting an advective scaling could shift the time scale by a factor of ∼ 10 (Lhuillier and Gilder, 2013).

#### 2.2. Analysis Method

From each model we compute Gauss coefficients gl,m(t) and hl,m(t) over time, t, at Earth's surface from magnetic field spectra at the CMB (e.g., Merrill et al., 1996). Although the dynamo model spectra are resolved out to harmonic degree l ≤ lmax = 64 we only compute Gauss coefficients out to lmax = 8 because larger harmonics contribute very little to the surface magnetic field.

From the Gauss coefficients we compute the geocentric axial dipole (GAD) intensity,

$$\mathbf{g}\_{AD}(t) = \sqrt{\mathbf{g}\_{1,0}(t)^2} \,\mathrm{,}\tag{1}$$

the full dipole intensity,

$$\mathbf{g}\_{\rm D}(t) = \sqrt{\mathbf{g}\_{1,0}(t)^2 + \mathbf{g}\_{1,1}(t)^2 + h\_{1,1}(t)^2},\tag{2}$$

the RMS non-axial dipole (NAD) intensity,

$$\mathbf{g\_{NAD}}(t) = \left[ \sum\_{l=1}^{l\_{\text{max}}} \sum\_{m=0}^{l} \left\{ \mathbf{g\_{l,m}}(t)^2 + h\_{l,m}(t)^2 \right\} - \mathbf{g\_{1,0}}(t)^2 \right]^{1/2}, \tag{3}$$

and the total RMS Gauss field intensity,

$$\mathbf{g\_{RMS}}(t) = \left[ \sum\_{l=1}^{l\_{\text{max}}} \sum\_{m=0}^{l} \left\{ \mathbf{g\_{l,m}}(t)^2 + h\_{l,m}(t)^2 \right\} \right]^{1/2}.\tag{4}$$

We also compute local magnetic field quantities on Earth's surface that are interpreted as synthetic observations. From the Gauss coefficients we compute the surface vector magnetic field components Bx, By, and B<sup>z</sup> from the magnetic potential 9(r, θ, φ) at a point on the surface at radius r = a,

$$B\_{\chi} = \frac{1}{a} \frac{\partial \Psi}{\partial \theta} \; , \; B\_{\mathcal{V}} = -\frac{1}{a \sin \theta} \frac{\partial \Psi}{\partial \phi} \; , \; B\_z = \frac{\partial \Psi}{\partial r} \tag{5}$$

where

$$\Psi(r,\theta,\phi,t) = a \sum\_{l} \sum\_{m} \left(\frac{a}{r}\right)^{l+1} P\_l^m(\cos\theta) (\mathcal{g}\_{l,m}(t) \cos m\phi)$$

$$+ h\_{l,m}(t) \sin m\phi) \tag{6}$$

and a = 2.8157 is the dimensionless radius of Earth in these models, θ is co-latitude, φ is longitude, and P m l are Schmidt normalized Legendre polynomials (Merrill et al., 1996). From the local field components we compute the local magnetic field intensity F as

$$F(\theta,\phi,t) = \sqrt{B\_x(\theta,\phi,t)^2 + B\_\jmath(\theta,\phi,t)^2 + B\_z(\theta,\phi,t)^2}.\tag{7}$$

Synthetic paleomagnetic observations are generated at φ = 0 and 6 co-latitudes in 15◦ increments: θ = 15◦ , 30◦ , 45◦ , 60◦ , 75◦ , and 90◦ .

Synthetic observations of magnetic direction are also generated at these locations. These include local magnetic declination D

$$D = \tan^{-1} \left[ \frac{B\_{\mathcal{Y}}}{B\_{\mathcal{X}}} \right] \tag{8}$$

inclination I,

$$I = \tan^{-1} \left[ \frac{B\_z}{\sqrt{B\_x^2 + B\_y^2}} \right] \tag{9}$$

angular pole distribution α<sup>95</sup> with 95% probability,

$$a\_{95} = \cos^{-1}\left[1 - \frac{N-R}{R}\left\{\left(\frac{1}{0.05}\right)^{\frac{1}{N-1}} - 1\right\}\right] \tag{10}$$

where N is the number of observations, and Fisher's precision parameter k

$$k = \frac{N - 1}{N - R} \tag{11}$$

where the resultant is

$$R = \left[ \left( \sum l\_i \right)^2 + \left( \sum m\_i \right)^2 + \left( \sum n\_i \right)^2 \right]^{1/2} \tag{12}$$

and the directional cosines are l<sup>i</sup> = cos D<sup>i</sup> cosI<sup>i</sup> , m<sup>i</sup> = sin D<sup>i</sup> cosI<sup>i</sup> , and n<sup>i</sup> = sin I<sup>i</sup> (Merrill et al., 1996).

The "true" or "known" (dimensionless) dipole moment pDM is

$$
\mathfrak{p}\_{\rm DM} = 4\pi a^3 \mathfrak{g}\_{\rm D} \tag{13}
$$

where g<sup>D</sup> is from (2). The "true" geocentric axial dipole moment pADM is

$$
\mathfrak{p}\_{ADM} = 4\pi a^3 \mathfrak{g}\_{AD} \tag{14}
$$

where gAD is from (1). These "true" dipole moments pDM and pADM are obtained directly from the dynamo model and will be compared to synthetic or "virtual" observations that attempt to infer these quantities at the surface. The "virtual" dipole moment (VDM) is the amplitude of a geocentric dipole that gives rise to an observed magnetic field vector at the surface (Tauxe and Yamazaki, 2015) and is computed from the local magnetic intensity F and inclination I by (Merrill et al., 1996)

$$\rho\_{\rm VDM} = 2\pi a^3 F (1 + 3\cos^2 I)^{1/2}.\tag{15}$$

The virtual axial dipole moment pV ADM is the amplitude of a geocentric axial dipole that gives rise to an observed magnetic intensity at a known co-latitude θ<sup>k</sup>

$$p\_{\rm VADM} = 4\pi a^3 F (1 + 3\cos^2 \theta\_k)^{-1/2}.\tag{16}$$

For a GAD, the co-latitude and magnetic inclination are related by

$$\tan I = 2 \cot \theta. \tag{17}$$

Even with constant control parameters and boundary conditions the dynamo models produce a highly time variable magnetic field so that there is no single "true" moment for a single model, but rather a stationary time average with some characteristic fluctuations about that average. Therefore, it is important to consider the role of time averaging in producing a paleomagnetic direction or intensity. The length of the time average τ applied to the dynamo model time series could be interpreted as the time over which a series of paleomagnetic observations are averaged to get a single data point in time (e.g., segment of a sedimentary sequence), or the time over which the magnetic carrier obtains a remnant signal of the ambient field (e.g., cooling of a magmatic unit below its Curie temperature). We average each quantity of interest over a number of smoothing times τ from 5 to 500 kyr. For each τ the dynamo model time series is split into N = 1t/τ sub series, where 1t is the total length of the model in kyr. Within each sub-series a running mean is computed following the method of Davies and Constable (2014):

$$
\overline{\chi\_i} = \overline{\chi\_{i-1}} \frac{(i-1)}{i} + \frac{1}{i} \chi\_i \; ; \; i \ge 1 \tag{18}
$$

where x<sup>i</sup> is some output from the dynamo model at the ith sampling index within each sub-series. The final running average value x(τ ) is computed over a time window of length τ for each sub series and corresponds to the value in (18) at the final index in each sub series. Finally an average and standard deviation of these N running averages is computed. Both the "true" and "observed" quantities will be time averaged the same way. Dynamo model output quantities have an output sampling frequency of about once every 1 kyr for all models.

#### 3. RESULTS

We apply the analysis methods described above to the suite of dynamo simulations summaried in **Table 1**. In this section we investigate how long a time base-line of observations must be averaged to obtain a pure geocentric axial dipole (GAD) field, and how this time baseline depends on the dynamics of the model. Finally we investigate how local virtual observations of inclination, intensity, and inferred dipole moment compare to the true dipole values for the suite of dynamos.

#### 3.1. An Earth-Like Dynamo Model

To demonstrate that these models are in a relevant region of parameter space, we first focus the details of an "Earth-like" model with dT/dr<sup>i</sup> = 4, which we will refer to as "model 4." **Figure 1** shows the time series of dipole co-latitude and the axial dipole Gauss coefficient gl=1,m=<sup>0</sup> for model 4. This model reverses 20 times over 7.6 Myr (2.61 reversals per Myr) and the dipole spends about equal time in each hemisphere. It is nominally "Earth-like" according to the definition of Christensen et al. (2010) with an axial to non-axial dipole ratio "AD/NAD" of 0.29, an odd to even ratio "O/E" of 0.87, zonal to nonzonal ratio "Z/NZ" of 0.05, and flux concentration of 2.37, all within one standard deviation of the Earth-like values of 1.4, 1.0, 0.15, 1.5 respectively. Using the standard deviations expected by Christensen et al. (2010), these values produce a summary rating of χ <sup>2</sup> = 6.88, which is considered "marginally" Earthlike. Note that these statistics are time variable so that the model can be more or less Earth-like over time, and that we are reporting the time averaged statistics over the entire length of the run (7.6 Myr). These magnetic field statistics give us some confidence that our models produce magnetic fields that are "Earth-like" at the largest scales even though they are many orders of magnitude from the Earth in several non-dimensional parameters. A recent comparison of "Earth-like" dynamos that span a huge range in control parameters demonstrates that the large scale features and low frequency variability can be captured even when the small scale dynamics are not resolved (Aubert et al., 2017). Nevertheless, our results should be compared to higher resolution simulations in the future.

Also shown in **Figure 1** are time series of magnetic inclination and declination as observed at an arbitrary point on Earth's surface: at the equator θ = 90◦ and φ = 0, referred to as Ieq and Deq. Synthetic observations generated in this way will be analyzed from a paleomagnetic perspective and compared to the true solutions. Next we will test the ability to retrieve the true dipole directions and intensities from a time series of synthetic paleomagnetic observations.

#### 3.2. Effects of Time Averaging

**Figure 2** shows an example of the time series smoothing analysis in (18) applied to g1,0 from model 4 for four smoothing lengths τ . Clearly in this occasionally reversing model g1,0 has longterm variability on Myr time scales that is not captured even by smoothing over 500 kyr (**Figure 2D**). However, RMS quantities may converge to stationary values faster.

**Figure 3** shows the mean and standard deviation of the running mean of four output quantities (i.e., x(τ )) from model 4 for a range of τ . **Figure 3A** shows that the average of g1,0(τ ) (i.e., the average of the running averages of g1,0 from all sub series) converges to zero for all choices of τ , as expected for this reversing model that spends roughly equal time in normal (g1,0 < 0) and reversed (g1,0 > 0) polarity states (**Figure 1**). The relatively large standard deviation of g1,0(τ ) reflects the long-term variability that is not averaged out within each sub series. This is expected if τ is less than the longest time scale of intrinsic variability of the model. **Figure 3B** shows that the average gNAD(τ ) from (3) is also near zero for all τ , implying that all other field harmonics (everything other than g1,0) individually balance to zero. The standard deviation of gNAD(τ ) also approaches zero at large τ , implying that nonaxial dipole fields more consistently balance out over longer time averaging than g1,0. **Figure 3C** shows that the running mean of RMS axial dipole gAD(τ ) from (1) and RMS dipole gD(τ ) from (2) are stationary and non-zero over all τ . **Figure 3D** shows that the local magnetic amplitude at the equator Feq is similarly stationary and non-zero over all τ and similar in

amplitude to gAD. A more comprehensive comparison between observed and true intensities as a function of co-latitude is below.

**Figure 4** shows the average of running averages of four equatorial (θ = 90◦ , φ = 0) observations of magnetic pole-related quantities from model 4: declination Deq from (8),

inclination Ieq from (9), α<sup>95</sup> from (10), and Fisher's precision parameter k from (11). The average of Deq and Ieq hover around zero, which is the expected orientation of a geocentric axial dipole observed at the equator. The standard deviation of Deq and Ieq decrease steadily with τ as the non-axial magnetic terms balance out, similar to the trend in gNAD in **Figure 3B**. The angular spread α<sup>95</sup> in the paleomagnetic pole decreases rapidly with τ while the precision parameter k plateaus around 15. A vertical dashed line is drawn at the largest <sup>τ</sup> where <sup>α</sup><sup>95</sup> <sup>&</sup>lt; <sup>10</sup>◦ , which is a typical threshold value for computing a paleomagnetic pole (e.g., Van der Voo, 1990). For model 4 this occurs at τ = 30 kyr, implying that to obtain a stable magnetic pole orientation for this Earth-like model the local field must be averaged over about 30 kyr.

So far we have investigated synthetic observations at the equator of a single Earth-like dynamo model with a basal heat flux of dT/dr<sup>i</sup> = 4. Next we examine how synthetic observations of the time averaged field depend on the basal heat flux and site co-latitude θ.

#### 3.3. Dynamo Regimes

How does the dynamical regime of the core dynamo influence synthetic observations at the surface? How do they differ (if at all) from the known solutions? These questions can addressed by investigating the suite of dynamos that span dynamical regimes from stable non-reversing at low basal heat flux (dT/dr<sup>i</sup> = 0.8) to reversing non-dipolar at high basal heat flux (dT/dr<sup>i</sup> = 12). The major dynamo transition from dipolar non-reversing models to non-dipolar reversing models occurs around dT/dr<sup>i</sup> = 3. This transition is apparent in **Figure 5A** where the time average of the volume averaged magnetic energy (ME) drops into a minimum due to a weakening of the axial dipole, **Figure 5B** where timeaveraged gAD drops by a factor of ∼ 4, **Figure 5C** where reversals begin, and **Figure 5D** where the time-averaged axial dipolarity (i.e., time average of g1,0(t)/gRMS(t)) drops below ∼ 0.5.

Interestingly, **Figure 5A** shows that volume averaged ME drops to a minimum at the onset of reversals (dT/dr<sup>i</sup> = 3) and then increases with heat flux nearly in parallel with kinetic energy (KE) as heat flux increases. Because of the preference for low harmonic degree fields at the surface, the decrease in the dipole dominates the total surface magnetic field, leading to a sudden drop in gAD at the reversing onset and a floor of gAD = 0.05 for more energetic models (**Figure 5B**). The gAD floor may indicate saturation of the dipole field where generation of a stronger dipole by faster convective velocities is balanced by turbulent disruption of the large scale field. Dipole reversal frequency increases with basal heat flux (**Figure 5C**), implying that dT/dr<sup>i</sup> is a proxy for reversal frequency and dipole stability in these models. The plateau in reversal frequency around 8/Myr is an artifact of our requirement that a reversal be bracketed by stable periods longer than 20 kyr, which become less common as the heat flux increases.

#### 3.4. Obtaining Stable Poles and Intensities From Synthetic Observations

Next we apply the local paleomagnetic analysis methods from section 2.2 to all models with the goal of quantifying how the

time required to obtain a stable paleomagnetic pole and intensity depends on the dynamical state of the core. We define the critical smoothing time τcrit as the running mean length where α<sup>95</sup> falls below a threshold value of either 5◦ or 10◦ . This is the length of time averaging needed to obtain a stable virtual paleomagnetic pole from continuous observations at a single location.

coefficients for the full dipole gD/gRMS and the axial dipole g1,0/gRMS.

**Figure 6A** shows that the critical smoothing time τcrit increases for more energetic dynamos driven by larger basal heat fluxes. A threshold of <sup>α</sup><sup>95</sup> <sup>&</sup>lt; <sup>10</sup>◦ requires τcrit of 20–40 kyr, while a threshold of α<sup>95</sup> < 5 ◦ requires τcrit of 40 − 150 kyr. **Figures 6B–D** show the average running mean of several other dynamo statistics computed at τcrit. Surprisingly **Figure 6D** shows that gNAD(τcrit) does not converge to zero for stable dipolar models (dT/dr<sup>i</sup> < 4), which implies that directional scatter around the paleomagnetic pole converges faster than the intensity of the non-dipolar magnetic field. More generally, this implies that longer running averages are needed to converge to stationary intensities than directions.

#### 3.5. Synthetic Intensities and Inclinations

Next we analyze synthetic dipole moment observations at 6 points on the surface (φ = 0 and θ = 15◦ , 30◦ , 45◦ , 60◦ , 75◦ , 90◦ ). Running averages of the observed magnetic inclination I from (9) and magnetic intensity F from (7) are compared to their "true" values in **Figure 7** for a chosen averaging time of τ = 500 kyr. Each point in **Figure 7** is an average of the running average for each quantity. We take absolute value of the inclination so that the time average converges to a non-zero number for reversing models and limits the range of possible inclinations observed in the northern hemisphere to [0◦ , 90◦ ].

At low basal heat flux the dipole is dominant and stable, producing time-averaged synthetic observations of I and F close to their true values (**Figure 7**). As heat flux increases nondipolar field components increasingly contribute to the time averaged local magnetic vector (**Figure 7C**). Contamination by these non-GAD components may be equivalent to adding a randomly oriented vector to the GAD field, which skews I toward 45◦ , the average inclination of a randomly oriented vector in one quadrant. This contamination causes the inclination to be severely underestimated at low co-latitudes (θ = 15◦ ) and overestimated at high co-latitudes (θ = 90◦ ), with residuals that increase with heat flux (**Figures 7A,B**).

For all basal heat fluxes, synthetic observations of magnetic intensity F tend to systematically over-estimate the intensity of a pure GAD (**Figure 7C**). This over-estimate increases with heat flux as non-GAD components contribute more to the time averaged field. This overestimate is larger at low co-latitude (θ = 15◦ ) and weaker near the equator (θ = 90◦ ) (**Figures 7C,D**). Note that the residual time-averaged inclinations (**Figure 7B**) and intensities (**Figure 7D**) are both non-zero even for nonreversing dynamos (dT/dr<sup>i</sup> < 3) where the presence of nondipolar field components contaminate the time-averaged field.

Another way to test adherence to a GAD field is by plotting the average observed inclinations I vs. the known co-latitude θk (**Figure 8**), where the pure GAD relationship is known from (17). Similar to **Figure 7** this comparison also shows that the low heat flux models adhere closely to a pure GAD field but that

as heat flux (and reversal frequency) increase they stray away from the GAD curve (**Figure 8**). Also as seen in **Figure 7A**, I tends to underestimate GAD at low θ and overestimate at high θ (**Figure 8**). Mixtures of GAD with the axial quadrupole G2 = g2,0/g1,0 and axial octopole G3 = g3,0/g1,0 describe the high heat flux models slightly better than pure GAD, but only at high θ<sup>k</sup> .

#### 3.6. Synthetic Dipole Moments

Next we compute the virtual full (pV DM) and axial (pV ADM) dipole moments inferred from the time-averaged synthetic observations of F and I using (15,16), respectively (**Figure 9**). pV ADM only depends on one observed quantity, intensity F, and subsequently the inferred values of pV ADM (**Figures 9C,D**) follow F (**Figures 9C,D**). On the other hand, pV DM depends on both observed quantities, F and I, and the resulting residual from the known pDM has a more complicated dependence on both intensity and inclination (**Figure 9B**).

The ratio of virtual dipole moment pV DM in (15) to true dipole moment pDM in (13), or dipole moment bias, is

$$b\_{DM} = \frac{p\_{VDM}}{p\_{DM}} = \frac{1}{2} \frac{\overline{F} (1 + 3 \cos^2 \overline{I})^{1/2}}{\overline{g\_D}} \tag{19}$$

and similarly for the ratio of virtual axial dipole moment pV ADM in (16) to true axial dipole moment pADM,

$$b\_{ADM} = \frac{p\_{VADM}}{p\_{ADM}} = \frac{\overline{F}}{\overline{g\_{AD}}(1 + 3\cos^2\theta\_k)^{1/2}}.\tag{20}$$

**Figure 10** shows that bDM tends to scatter around the true value (i.e., bDM = 1) for non-reversing dipolar models, but then trends down to as low as ∼ 0.6 for reversing models. The scatter in bDM about this trend increases systematically with reversal frequency Rf . The decrease of bDM with R<sup>f</sup> in **Figure 10** is caused mainly by a drop in F/g<sup>D</sup> in (19) because I tends toward a constant (∼ 45◦ ) as R<sup>f</sup> increases (**Figure 7A**). The decrease in F/gD, in turn, is caused by the local intensity becoming more and more contaminated by non-dipole field components with arbitrary sign that occasionally balance out as heat flux increases, so that on average F < gD. We have found these ratios are roughly the same for τ = 50 , 100 , and 500 kyr, indicating that this bias does not depend sensitively on the length of time averaging.

On the other hand, the axial dipole moment ratio bADM increases with reversal frequency from bADM ∼ 1 at R<sup>f</sup> = 0 up to <sup>b</sup>ADM <sup>∼</sup> 1.4 at <sup>R</sup><sup>f</sup> <sup>=</sup> 8 Myr−<sup>1</sup> (**Figure 11**). In contrast to bDM, the increasing trend of bADM with R<sup>f</sup> is caused by gAD decreasing faster than F in (20) as a function of R<sup>f</sup> . This bias can be attributed to the GAD assumption that the field is purely axial when it is not and non-axial dipolar field components contributing to increase local intensity, so that on average F > gAD.

Combining these two constraints, pV DM < pDM from **Figure 10** and pV ADM > pA DM from **Figure 11**, gives

$$\sqrt{\mathbf{g}\_{1,0}^2(1+3\cos^2(\theta\_k))} < \overline{F} < \sqrt{(\mathbf{g}\_{1,0}^2 + \mathbf{g}\_{1,1}^2 + h\_{1,1}^2)(1+3\cos^2(\overline{\theta}))}\tag{21}$$

where we have converted the observed I to θ using the GAD assumption in (17). At mid-latitudes where the observed paleoco-latitude is similar to the true co-latitude (**Figure 7A**), (21)

θk (see legend), solid lines are the GAD inclination, computed from (17) with θk . (B) Residual between observed and GAD inclination. (C) Observed running average magnetic intensity F (circles) from (7), colors correspond to site co-latitude θk (see legend), and solid lines are the GAD intensities. (D) Residual between observed and GAD intensity.

reduces to

$$\sqrt{g\_{1,0}^2} < \sqrt{\frac{2}{5}} \overline{F} < \sqrt{g\_{1,0}^2 + g\_{1,1}^2 + h\_{1,1}^2} \quad \text{for } \overline{\theta} \approx \theta\_k \approx 45^\circ. \tag{22}$$

The inequality in (22) implies that local fields of arbitrary sign combine to decrease F slightly less than the full dipole field (l = 1, m = 0, ±1) but to increase it slightly more than the axial dipole (l = 1, m = 0) alone. In this sense, the time averaged intensity at mid-latitudes is bracketed by the GAD and full dipole field intensity.

Lastly, in an attempt to extract a scaling law that describes these effects we assume a bi-linear form for the dipole moment bias b as a function of observed (or assumed) reversal frequency Rf , and observed (or inferred) paleo-co-latitude θ,

$$b\_{\rm DM} = b\_{\rm DM0} + b\_{\rm DM1}R\_f + b\_{\rm DM2}\theta + b\_{\rm DM3}R\_f\theta \tag{23}$$

$$b\_{ADM} = b\_{ADM0} + b\_{ADM1}R\_f + b\_{ADM2}
\theta + b\_{ADM3}R\_f \theta \quad \text{(24)}$$

where the fit coefficients extracted from the synthetic observations in **Figures 10**, **11** are shown in **Table 2**. This scaling law could potentially be used to adjust estimates of the dipole moment from local measurements of intensity and inclination. Next we will apply this bias scaling law to published paleointensity data.

#### 4. PINT PALEOINTENSITY DATABASE

The PaleoINTensity (PINT) database of Biggin et al. (2009) is a compilation of absolute paleointensity measurements using Thellier-type experiments with each site mean produced from at least 3 individual measurements and a standard deviation that is not more than 25% of the mean. The database (downloaded from http://earth.liv.ac.uk/pint/) was last updated in 2015 to include quality QPI factors (Biggin and Paterson, 2014), contains a total of 362 dated paleointensity measurements, 147 of which are virtual dipole moments (VDM) and 215 virtual axial dipole moments (VADM) (**Table 3**). The data span 522–3458 Ma. The difference between the PINT VDM and VADM data are the time over which the paleomagnetic inclination is averaged: VDM's use site mean inclinations while VADM's use study mean inclinations

over much a longer time range. Since both quantities involve a paleomagnetic inclination we will adjust both as if they were

VDM's, as in (15). **Figure 12** shows VDM and VADM from the PINTQPI paleointensity database for all QPI values. We apply an inverse distance squared smoothing (Algeo, 1996; Driscoll and Evans, 2016) to the dataset to produce a running mean µ(t) V(A)DM profile (black line), standard deviation σ(t) (dark gray), and standard error e(t) (light gray) with smoothing time scale of 30 Myr. We see no obvious trend in **Figure 12** and the paleointensity appears more or less flat over this time span with significant variability.

Some of these data, however, may be misleading. In an attempt to avoid moments derived from non-ideal analyses we consider a subset of the data with QPI ≥ 3 in **Figure 13** (similar to Biggin et al., 2015), and with QPI ≥ 4 in **Figure 14**. Also in **Figures 13**, **14** we have applied the dipole moment adjustments found in **Figure 10** by dividing all VDM's and VADM's by 0.8, which are rough estimates of the moment bias we found from synthetic observations of numerical dynamos. The adjusted values are marked by filled circles with a vertical line connecting them to their original (pre-adjusted) values. The same inverse distance squared smoothing is then applied to the adjusted values.

For QPI ≥ 3 a jump in V(A)DM occurs around 1.3 Ga (**Figure 13**), which was interpreted as a signature of inner core nucleation (ICN) by Biggin et al. (2015). The smoothed average then drops back down to normal, pre-jump intensities, which is not predicted from dynamo models of ICN (Aubert et al., 2009; Driscoll, 2016; Landeau et al., 2017). Using an even more stringent criteria of QPI ≥ 4, which includes even less data, does not show a peak around 1.3 Ga (**Figure 14**). By increasing the QPI cutoff from 3 to 4 a pivotal set of intensities by Thomas (1993) are removed at 1.3 Ga, which results in a relatively flat or possibly slowly decreasing moment over Precambrian time. The time averaged Precambrian dipole moment (combined VDM and VADM) is in the range 58 − 62 ZAm<sup>2</sup> in **Table 3**, slightly higher than the estimate by Sprain et al. (2018). The other striking observation from **Figures 12**–**14** is the number and length of time gaps in the data. In particular there are 4 paleointensities in the Neoproterozoic with QPI ≥ 3 and none with QPI ≥ 4.

Excluding so much data may seem overly harsh but there is reason to be careful. It is well known that multidomain components can produce concave-up demagnetization curves (Arai plots) (e.g., Shcherbakova et al., 2000; Chauvin et al., 2005; Biggin et al., 2007; Tauxe and Yamazaki, 2015; Smirnov et al., 2016; Sprain et al., 2018), which could result in an over estimate of the primary intensity if not heated sufficiently. As discussed by Smirnov et al. (2016) and Sprain et al. (2018), examples of this effect may include the Keweenawan rocks dated at 1.1 Ga, which exhibit anomalously high paleointensities and puzzling asymmetries between the normal and reversed polarity sections (Pesonen and Halls, 1983). These asymmetries could be created by rapid changes in paleolatitude or paleo-secular variation (e.g., Swanson-Hysell et al., 2009; Sprain et al., 2018). The Keweenawan intensities of Pesonen and Halls (1983) are also 2–3 times higher than the Tudor Gabbros (Yu and Dunlop, 2001), Abitibi dykes (Macouin et al., 2003), and central Arizona diabase sheets (Donadini et al., 2011), all of which were emplaced within about 50 Myr of the Keweenawan. The period of high paleointensity found by Pesonen and Halls (1983) may

TABLE 2 | Table of coefficients in (23) and (24), where n refers to the coefficient number.


be a real transient, which would be consistent with modern levels of intrinsic geodynamo variability, or in fact may be an artifact (Yu and Dunlop, 2001; Valet, 2003). Similarly, a series of high paleointensity measurements from the Gardar Basalts in southern Greenland (Thomas, 1993; Thomas and Piper, 1995) may have misinterpreted several low temperature multidomain magnetizations as a primary single domain paleointensity. In light of these possible issues we find QPI ≥ 4 is sufficient to filter out the contentious data. We note that the paleointensity data are heterogeneous in the various experimental methods and types of rocks used, which complicates their interpretation, and that these issues with paleointensity measurements are an additional complicating factor to consider along with the observed bias found in the synthetic observations above.

Several limitations of our attempt to adjust or unbias the PINT data by the results of our synthetic analysis above should be mentioned. In adjusting the PINT V(A)DM values we divided the VDM and VADM values by a constant 0.8, whereas a more precise adjustment would use the fit in (19,20) with an estimate of the reversal frequency R<sup>f</sup> and site paleolatitude θ. Although this approach will be difficult for Precambrian data where the paleolatitude is often not independently known and the reversal frequency record is discontinuous (e.g., Pavlov and Gallet, 2010;


Columns are number of data N and dipole moment p in [ZAm<sup>2</sup> ] for all data, VDM only, and VADM only. Dipole moments divided by 0.8 labeled as "Adjusted."

Biggin et al., 2011; Gallet et al., 2012). Correlating polarity ratio, which is easier to obtain, with polarity reversal frequency may be a way to extend the record (e.g., Driscoll and Evans, 2016). More readily our predicted biases could be applied to the better constrained paleointensity record over the last 200 Myr, where the reversal frequency and paleolatitude are better known, in order to reexamine the prediction of an inverse relationship between reversal frequency and paleointensity (e.g., Driscoll and Olson, 2011; Sprain et al., 2016).

#### 5. CONCLUSIONS

We have generated synthetic magnetic observations from numerical dynamos that span a range of dynamical regimes and reversal frequencies to investigate time averaged magnetic field orientations and intensities. The range of dynamo regimes found, from stable non-reversing dipolar regimes to reversing non-dipolar regimes, are driven by models that span a factor of 10 increase in bottom boundary heat flux. We find that running averages over 20 − 40 kyr are needed to obtain stable

paleomagnetic poles with <sup>α</sup><sup>95</sup> <sup>&</sup>lt; <sup>10</sup>◦ , and over 40 − 150 kyr for α<sup>95</sup> < 5 ◦ . To obtain stable paleomagnetic poles we find that models with higher heat flux and more frequent polarity reversals require longer time averages, similar to previous studies (McMillan et al., 2001; Davies and Constable, 2014). However, we also find that obtaining stable intensities requires longer time averaging than directions.

Surprisingly we find that 500 kyr running averages of local field intensity and inclination produce underestimates of VDM by as much as 50% and overestimates of VADM by as much as 150% with increasing basal heat flux and reversal frequency. These biases are caused by the running averaged local field intensity F having an intermediate intensity between the RMS axial dipole and full dipole intensities. Similar to Lhuillier and Gilder (2013) we find significant dependence on site latitude with high (low) latitudes producing high (low) moment estimates compared to mid-latitudes. These biases remain even for averages of running time averages over 500 kyr. We derive a scaling law for paleointensity bias as a function of reversal frequency and site latitude that could be applied in the future to records where both are known.

Our method for computing time-averaged synthetic observations is not unique and the details of this time averaging likely influence the conclusions. When the underlying field itself is changing on 10<sup>3</sup> − 10<sup>6</sup> yr time scales the length of time over which the field is averaged will change the estimates. Also, for non-linear quantities like inclination the order of operations might be important. We chose to compute a running average of the inclination over time, but computing inclination from a running average of the underlying field vectors may produce different results. The length of the running time average might also be important. Our choice here to compute running averages over 500 kyr is well suited to addressing the GAD assumption that is used to infer slow (Myr) tectonic motions, but shorter time scales may be more applicable for investigating secular variation, reversal rates, and variability of the intensity. Ultimately the averaging scheme should be specific to the question being asked. The same non-uniqueness of the synthetic observations also applies to the site latitudes chosen. We sampled the dynamos at 6 regularly spaced latitudes and a single longitude. Do observations spread out randomly over Earth's surface do a better job of sampling the underlying stationary field? This should be explored in future studies. In particular, the influence of site latitude and reversal frequency (or dipolarity) on the inferred paleo-latitude is important during so-called "snowball Earth" events where glacial deposits are interpreted to have occurred at equatorial latitudes (e.g., Evans and Raub, 2011). In the future, synthetic observations from numerical dynamos could explore in greater detail how the method of time averaging and site location affects their ability to recover the underlying global magnetic field.

Finally, applying the paleointensity adjustments found by the dynamo models to the PINT paleointensity record produces little change to the overall trend of a relatively flat intensity over the last 3.4 Ga. Filtering data with quality factors ≥ 3 produces a tentative jump in paleo-dipole moment around 1.3 Ga, but disappears for quality factors ≥ 4, where a possible secular decrease in intensity in seen from 2.7 to 1 Ga. A more careful analysis could be applied during periods when the reversal frequency, or some proxy for reversal frequency (such as secular variation or polarity ratio), and site latitude is known. Correcting for intensity and inclination bias may be important for identifying trends and events (like inner core nucleation) in the paleointensity record. Future analysis could be extended to synthetic observations at more locations on the surface, quantifying secular variation, and identifying higher order magnetic components. Models at lower Ekman number and with a range of boundary conditions should also be analyzed to investigate if these results are sensitive to this region of parameter space.

#### AUTHOR CONTRIBUTIONS

PD conceived the project, produced simulations and data analysis, and wrote the paper. CW assisted with software setup, data analysis, and writing.

#### REFERENCES


#### FUNDING

The authors thank The Carnegie Institution for Science for funding.

#### ACKNOWLEDGMENTS

We thank Andrew Biggin (University of Liverpool) for maintaining the PINT database used here and for a thorough review that significantly improved this paper. This article benefited from discussions with Chris Davies, Richard Bono, and Florian Lhuillier. Rayleigh (https://github.com/geodynamics/ Rayleigh), the dynamo code used here, has been developed by Nicholas Featherstone (University of Colorado) with support by the National Science Foundation through the Computational Infrastructure for Geodynamics (CIG), and was supported by NSF grants NSF-0949446 and NSF-1550901. Carnegie's Memex high performance computational cluster was used to generate all numerical dynamos. The data can be made available upon request.


Dipole (GAD) hypothesis in the Precambrian. Precambrian Res. 244, 33–41. doi: 10.1016/j.precamres.2013.10.009


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Driscoll and Wilson. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Recent Advances in Chinese Archeomagnetism

Shuhui Cai 1, 2, 3 \*, Lisa Tauxe<sup>1</sup> , Greig A. Paterson3, 4, Chenglong Deng2, 3, 5, Yongxin Pan3, 4, 5 , Huafeng Qin2, 3 and Rixiang Zhu2, 3, 5

*<sup>1</sup> Scripps Institution of Oceanography, University of California, San Diego, San Diego, CA, United States, <sup>2</sup> State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, China, 3 Institutions of Earth Science, Chinese Academy of Sciences, Beijing, China, <sup>4</sup> Key Laboratory of the Earth and Planetary Physics, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, China, <sup>5</sup> University of Chinese Academy of Sciences, Beijing, China*

The geomagnetic field is one of Earth's fundamental properties with a history of ∼3.5 Gyr. The field, generated in Earth's core is a window to the deep interior of Earth and may have played a key role in evolution of life on our planet. Materials on Earth's surface that contain magnetic minerals can record information about the geomagnetic field in which they formed. Fired archeological materials (e.g., pottery, brick, and burnt clay) are favorable recorders of the field, and have been widely employed to recover geomagnetic variations over periods of hundreds to thousands of years. The longevity of Chinese civilization and the abundant nature of archeological artifacts make Chinese archeomagnetism a promising source of data. The main work of Chinese archeomagnetism was carried out in the 1980s and 90s, followed by a break of more than a decade; in the 2010s activity resumed. In this paper, we review the development of Chinese archeomagnetism, including a summary of previous work, recent progress, remaining issues and future studies with the aim of promoting an understanding of archeomagnetic work in China and to guide the way for future studies. Here, we compile published data, including some data discovered in old publications that have not yet been included in paleomagnetic databases. We also establish the first, albeit preliminary, archeomagnetic reference curves (with 42 declination / inclination pairs and 76 / 192 archeointensities) for the geomagnetic field in China (ArchInt\_China.1a / ArchInt\_China.1b, ArchDec\_China.1, ArchInc\_China.1), which can be used for global comparison of the field and regional archeomagnetic dating.

#### Edited by:

*Joshua M. Feinberg, University of Minnesota, United States*

#### Reviewed by:

*Chuang Xuan, University of Southampton, United Kingdom Nobutatsu Mochizuki, Kumamoto University, Japan*

> \*Correspondence: *Shuhui Cai shc252@ucsd.edu*

#### Specialty section:

*This article was submitted to Geomagnetism and Paleomagnetism, a section of the journal Frontiers in Earth Science*

> Received: *01 September 2017* Accepted: *30 October 2017* Published: *14 November 2017*

#### Citation:

*Cai S, Tauxe L, Paterson GA, Deng C, Pan Y, Qin H and Zhu R (2017) Recent Advances in Chinese Archeomagnetism. Front. Earth Sci. 5:92. doi: 10.3389/feart.2017.00092* Keywords: archeomagnetism, China, geomagnetic reference curves, the Holocene, archeomagnetic dating

#### INTRODUCTION

Fired archeological artifacts, such as pottery, brick, burnt clay, furnace fragments and metallurgical slags, are favorable materials for recording the evolution of the geomagnetic field, owing to their suitable magnetic characteristics as well as their abundance and relative temporal continuity. Archeomagnetism contributes greatly to recovering the secular variation of the geomagnetic field during the Holocene, which has applications for exploring the geodynamo in Earth's interior (Tarduno et al., 2015; Terra-Nova et al., 2016; Davies and Constable, 2017) and establishing various global models [e.g., the CALS series (Korte and Constable, 2011; Korte et al., 2011; Constable et al., 2016), ARCH3k.1 (Korte et al., 2009), ARCH10k.1 (Constable et al., 2016), pfm9k (Nilsson et al., 2014), and SHA.DIF.14k (Pavón-Carrasco et al., 2014)]. Archeomagnetic studies can also be used to solve archeological issues, such as dating an artifact by comparing its recorded geomagnetic intensity and/or direction to a local geomagnetic reference curve (Aitken, 1990; Pavón-Carrasco et al., 2011; Carrancho et al., 2017; Peters et al., 2017), or testing the synchronicity of archeological units by comparing the geomagnetic information extracted from them (Carrancho et al., 2016). Archeomagnetic studies even have potential applications for exploring the relationship between positions of virtual geomagnetic poles and historical records of aurorae (Liritzis, 1988).

Archeomagnetic studies originated in France (Folgerhaiter, 1899; Chevallier, 1925; Thellier, 1938; Thellier and Thellier, 1959), followed by studies from other countries, including Japan (Hirooka, 1971; Sakai and Hirooka, 1986), the U.K. (Aitken et al., 1981) and Bulgaria (Kovacheva, 1980). Recent archeomagnetic studies are mostly concentrated in Europe (Gallet et al., 2002; Gómez-Paccard et al., 2012; Tema et al., 2012; Genevey et al., 2013; Hervé et al., 2013a,b; Kovacheva et al., 2014) and the Middle East (Ertepinar et al., 2012, 2016; Gallet et al., 2015; Shaar et al., 2016; Ben-Yosef et al., 2017), with a few publications from other areas such as Mexico (Guerrero et al., 2016), Africa (Mitra et al., 2013; Tarduno et al., 2015; Kapper et al., 2017) and Asia (Yu et al., 2010; Hong et al., 2013; Venkatachalapathy et al., 2013). China constitutes a huge part of Eastern Asia and has a civilization that spans thousands of years leaving abundant archeological artifacts. Archeomagnetic studies in this region are essential and feasible. In this paper, to promote an understanding of current archeomagnetic studies in China and to help guide future work, we summarize the development of Chinese archeomagnetism and establish the first archeomagnetic reference curves of geomagnetic field variations in China.

#### ADVANCES OF CHINESE ARCHEOMAGNETISM

Archeomagnetic studies in China were first carried out in the 1960s by Deng and Li (1965), who retrieved a few paleointensity and inclination data points from the Beijing area. This was followed by a paper, with the by-line of Paleomagnetism Laboratory of the Institute of Geology, Chinese Academy of Sciences (named the author "IGCAS" hereafter), which reported a number of paleointensity results from the Jokhang Temple in Lhasa, Tibet (IGCAS, 1977). A great quantity of Chinese archeomagnetic studies were carried out in the 1980s by Wei et al. (1980, 1981, 1982, 1983, 1984, 1986, 1987) from the Institute of Geophysics, Chinese Academy of Sciences, which constitute the majority of archeomagnetic data in China. A few results were then published in the 1990s (Tang et al., 1991; Yang et al., 1993a,b; Shaw et al., 1995, 1999; Batt et al., 1998; Huang et al., 1998). After that, archeomagnetic studies in China ceased for ∼15 years before revival in the current decade (Cai et al., 2014, 2015, 2016, 2017).

The quality of the twentieth century Chinese archeomagnetic data, especially paleointensity data, is uneven and sometimes hard to assess, which is due to lack of modern experimental techniques and/or relaxed selection criteria. The descriptions of experimental techniques for the paleointensity publications are summarized in Table S1, including a range of Thelliertype thermal methods (Thellier and Thellier, 1959; Coe, 1967; Yu et al., 2004), the Shaw method (Shaw, 1974) as well as microwave based methods (Walton et al., 1993). We also summarize if the studies reported basic statistical descriptions (σB, the paleointensity standard deviation as a measure of the consistency of sister specimens); if authors conducted partial thermal remanent magnetization (pTRM) checks to monitor alteration during experimental heating (Coe et al., 1978); if they considered the effect of TRM anisotropy (Aitken et al., 1981, 1988), by either applying an anisotropy correction (Veitch et al., 1984) or aligning the laboratory field to the direction of natural remanent magnetization (NRM); if they considered the possible bias caused by cooling rate effects (Dodson and McClelland-Brown, 1980; Halgedahl et al., 1980) and carried out a cooling rate correction (Genevey and Gallet, 2002). The data published in the twentieth century either have no or only partial quality controls, making unambiguous estimation of their reliability difficult. In contrast, the data published in the current decade are obtained from more rigorous modern experimental techniques with stringent data selection and openly available measurement data; these should therefore be more robust.

Most of the published data from China are archived in the GEOMAGIA50.v3 database (http://geomagia.gfz-potsdam. de) (Brown et al., 2015) with the exception of data published by Deng and Li (1965); IGCAS (1977); Wei et al. (1980), and Yang et al. (1993a,b) and the recent work of Cai et al. (2014, 2015, 2016, 2017). The former three are absent from the database probably because they were originally published in Chinese. Yang et al. (1993a) did not include their data list in the paper. However, results from Shaw et al. (1999) were obtained from the same batch of samples as in Yang et al. (1993a), but with different experimental techniques. Data from Yang et al. (1993b) were left out for some unknown reason. The last four papers were published in or after 2014 and data therein have not been included in the GEOMAGIA50.v3 database yet, but are already in the MagIC database (https://www2.earthref.org/ MagIC). In this paper, we collated the data from the publications not included in any database and listed them in Table S2 for the convenience of future work; these have also been uploaded into the MagIC database (https://earthref.org/MagIC/16283). The geographic locations of all published archeomagnetic data (both direction and intensity) are plotted in **Figure 1**.

#### COMPILATION OF PUBLISHED ARCHEOMAGNETIC DATA FROM CHINA

The archeointensity data published in recent years by Cai et al. (2014, 2015, 2017) were selected with different criteria, making the data quality inconsistent. In this paper, we combined all the measurement data from each paper and reanalyzed the data with the "Thellier Auto Interpreter" function incorporated in the Thellier GUI software (Shaar and Tauxe, 2013) included in the PmagPy software package (Tauxe et al., 2016). The strict selection criteria suggested by Cromwell et al. (2015) and named "CCRIT"

by Tauxe et al. (2016) were adopted. These criteria are established based on paleointensity study on modern lava flows in Hawaii where the historical field can be reproduced with reasonable accuracy (Cromwell et al., 2015). The parameters of CCRIT are listed in Table S3. The definition of each parameter is described in Cromwell et al. (2015) while the detailed explanation can be found in Shaar and Tauxe (2013) and Paterson et al. (2014) and references therein. We only provide a brief reminder of each statistics here: β is the normalized standard deviation of the slope of selected data points; DANG is the angle of the best-fit line deviated from the origin; MADfree is the unanchored maximum angular deviation of selected NRM data points; FRAC is the fraction of remanence used for calculating the paleointensity; SCAT is a Boolean that defines the allowed degree of scatter of the selected data points (including pTRM checks); | −→k ′ | is the absolute value of curvature of the data points used for determining the best-fit line (Paterson, 2011; Paterson et al., 2015); Gap Max is the maximum length of the normalized vector differences between consecutive NRM steps along the chosen segment; Nmin is the minimum number of accepted specimens; σ is the one-sigma standard deviation of site-mean intensity. The new results reanalyzed with CCRIT are listed in Table S4 (sample means) and Table S5 (specimen results with statistic values). The reanalyzed paleointensity data (Int\_new pub) as well as the data published in the twentieth century are plotted in **Figure 2A**. We calculated two Chinese archeointensity reference curves of the geomagnetic field: one (ArchInt\_China.1a) with only the reanalyzed data published in the 2010s (76 in total), assuring equal quality of data used for calculating the reference curve; the other (ArchInt\_China.1b) with both the reanalyzed data and selected old data published in the twentieth century (192 in total). Since limited statistic parameters were reported

reanalyzed archeointensity data published by Cai et al. (2014, 2015, 2017). "Dec\_new pub" and "Inc\_new pub" represent directional data published by Cai et al. (2016). The shaded area of each curve represents the one standard deviation coverage interval. All the directional data were relocated to the center of China (35◦N, <sup>105</sup>◦ lE).

in the old publications (Table S1), we can only use the most general selection criteria when selecting the old data. Only those data with age sigma less than 500 yr and standard deviation of the intensity less than 4 µT or 10% (the same requirement for sample mean as in CCRIT) were included. The reference curve was calculated with a parametric bootstrap and running average technique, following Cai et al. (2017) and Gallet et al. (2015). The procedure is: resample 1,000 times at each data point considering uncertainties of both age and virtual axial dipole moment (VADM) and then calculate the running average of the new dataset with a time window of 200 y shifted by 10 y (only time intervals including more than three data points were calculated). The one-sigma standard deviation (orange / light blue shadow in **Figure 2A**) is calculated as well. ArchInt\_China.1b is generally similar to ArchInt\_China.1a except for two obvious differences: 1) ArchInt\_China.1b smooths out the field low at ∼2200 BCE in ArchInt\_China.1a; and 2) ArchInt\_China.1b is higher than ArchInt\_China.1a between ∼100–600 CE because of high intensity data published by Wei et al. (1982, 1986). Special attention should be paid to these two time periods in the future work to resolve the source of these differences.

The published directional data in China are plotted in **Figures 2B,C**. Only data with both declination and inclination are included. In order to reduce the difference caused by locations, we relocated all the directional data to the center of China (35◦N, 105◦E) using the VGP relocation method (Noel and Batt, 1990): first, calculate the VGP with declination, inclination, site latitude and longitude; and then calculate the declination and inclination at the relocated location from the VGP assuming a dipolar field. All the new and old published data and the relocated data are presented in Cai et al. (2016). The directional data are scarce because in-situ materials are less likely to be preserved. The experimental techniques and influence factors of determining geomagnetic directions are less complicated than those of paleointensities, meaning that directional data are less prone to suffer from large biases away from the correct value. Actually, all the 95% confidence intervals (α95s) of the old published data are less than 10◦ , except one is 12.6◦ . Therefore, we calculated the preliminary Chinese archeodirection reference curves of the geomagnetic field (ArchDec\_China.1 and ArchInc\_China.1) with all the published data (42 declination / inclination pairs in total). The reference curves were calculated using the same technique as used for calculating the archeointensity reference curve. The running average data for all four reference curves are attached in Table S6.

We also calculated reference curves for adjacent areas of Japan and Korea using the same technique above (Figure S1). The Japanese reference curves were calculated with data from GEOMAGIA50.v3. The intensity data were selected with the same criteria as for the old published intensity data in China while for directional data only those with <sup>α</sup><sup>95</sup> less than 10◦ were included. Since there are only two directional data points between 6000 BC and 400 CE after selection, the Japanese directional reference curves start from 400 CE. The Korean intensity / directional curves were calculated with data from Hong et al. (2013) and Yu et al. (2010). The Korean intensity curve and predictions at the center of China (35◦N, 105◦E) of global models

[pfm9k.1a (Nilsson et al., 2014) and ARCH10k.1 (Constable et al., 2016)] are generally consistent with our new archeointensity curves except certain time periods, for example, the extremely low intensity at ∼2,200 BCE (Figure S1). The Japanese intensity curve agrees well with our new curve before ∼2,800 BCE and after ∼100 CE but deviates from our curve between them. The directional data are too sparse to discuss the consistency among different models before ∼1,100 BCE. After then the global models are consistent well with our new reference curves because they all mainly rely on the present published data (Figure S1A). Both the Japanese declination and inclination curves agree well with our new curves and the global models. However, the Korean declination curve fits our curve between ∼1,100 BCE and ∼500 CE but departs from our curve after then while the inclination curve deviates far away from the Chinese curve most of the time except after ∼1,200 CE (Figures S1B,C). The difference between 500 and 1,000 CE is probably caused by the absence of Korean data during that period. Discrepancies during other time periods can either due to geomagnetic local field anomalies or lack of robust data from both areas.

Our new archeomagnetic reference curves are preliminary due to the lack of a large number of robust data, for example, the paucity of paleointensity data at certain times (i.e., before ∼1,300 BCE, ∼1,000 BCE-500 BCE, 1 CE-500 CE) and the overall scarcity of directional data. More high-quality data are required in the future to enhance the precision and resolution of the reference curves and to promote global comparison.

#### CURRENT ISSUES AND THE FUTURE OF CHINESE ARCHEOMAGNETISM

The Chinese archeointensity data published during the twentieth century are scattered because of a lack of modern experimental techniques and/or relaxed selection criteria—their reliability is difficult to assess. The data published in recent years have a higher precision, but are insufficient in both space and time to define a robust reference curve. Most of the data are from Eastern China and data from other regions are lacking (**Figure 1A**). Similarly, data from certain time periods (e.g., before ∼1,300 BCE, ∼1,000 BCE-500 BCE, 1 CE-500 CE) are sparse (**Figure 2A**). The directional data from China, especially those with full directions (both declination and inclination) are quite rare (**Figures 2B,C**) because of the difficulty of preserving in-situ materials since they were fired. Data with full vector information (both direction and paleointensity) of the field are certainly more deficient. Furthermore, dating of the published data mainly relies on archeological context and the age resolution of some data is a limiting factor (e.g., data around 1,300 ± 300 BCE in **Figure 2A**).

Filling these spatial and temporal data gaps is a frontier for Chinese archeomagnetism. Future workers should carefully focus on obtaining more high-quality data. Considering the sparseness of directional data, special attention should be paid to collecting more oriented in-situ materials. Fortunately, relics as kilns or hearths are not rare and a number of them have not been moved since they were fired last. All we need to do is cooperating with the archeologists and collecting oriented samples from these sites before they were destroyed. The enhancement of data quality must include two aspects: (1) increasing the data precision by adopting modern experimental techniques and stringent selection criteria and (2) reinforcing age constraints by combining multiple dating techniques (e.g., archeological context, radiocarbon dating and stratigraphic information)—for example, dating techniques adopted in Cai et al. (2016). We also strongly encourage the deposition of original measurement data (both inside and outside China) in an accessible database such as MagIC. This can allow older datasets to be carefully reassessed and avoids the blanket rejection of older data.

With the accumulation of new reliable data, the Chinese archeomagnetic reference curves can be updated and become a more precise tool for archeological and Earth scientists. To achieve these ambitious goals, more attention should be paid to enhancing the cooperation between paleomagnetists and archeologists. For example, archeologists can assist paleomagnetists in collecting abundant excellent samples while paleomagnetists can lend their expertise in archeological applications (e.g., estimating firing temperature of artifacts, testing the synchronicity of various archeological units, determining if artifacts from an archeological unit are in-situ, among others).

#### SUMMARY

Archeomagnetic research is one of the most efficient methods to explore the detailed secular variation of the geomagnetic field during the Holocene and has applications to geodynamics, global modeling, establishing regional reference curves, archeomagnetic dating and other archeological and Earth science issues. In this paper, we have outlined the current state of archeomagnetic studies in China and the future challenges that this discipline faces—namely a scarcity of high quality data measured with modern standard of reliability. Nevertheless, by reanalyzing the most recently acquired datasets, we provide a number of consistent high-precision archeointensity data

#### REFERENCES


from China. We have established the first, albeit preliminary, Chinese archeomagnetic reference curves (ArchInt\_China.1a / ArchInt\_China.1b, ArchDec\_China.1, ArchInc\_China.1). These can be used for regional geomagnetic comparison and archeomagnetic dating and form a basis on which future studies can expand by targeting time periods and regions that are currently underrepresented. Future work that focuses on filling these spatial and temporal gaps should make great efforts to obtain more high-quality (both high data precision and well age constraints) data, which can be achieved by enhancing the cooperation between paleomagnetists and archeologists to the benefit of both communities.

#### AUTHOR CONTRIBUTIONS

SC, LT, and GP wrote the paper; RZ, CD, YP, and HQ helped designing research and sampling as well as revising the paper.

#### FUNDING

This work was supported by the NSFC grants 41504052, 41621004, 41374072 and L1524016, the CAS priority project XDB18010203. LT acknowledges support from NSF Grants No. EAR1520674 and EAR1547263.

#### ACKNOWLEDGMENTS

We thank all the archeologists we have worked with (Guiyun Jin, Jianming Zheng, Xuexiang Chen, Fei Xie, Wei Chen) for providing the samples.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2017.00092/full#supplementary-material

to the archeological and volcanic database. Earth Planets Space 67, 83. doi: 10.1186/s40623-015-0232-0


limekiln at pinilla del valle site (Madrid, Spain). Archaeometry 59, 373–394. doi: 10.1111/arcm.12245


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Cai, Tauxe, Paterson, Deng, Pan, Qin and Zhu. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## The First Catalog of Archaeomagnetic Directions From Israel With 4,000 Years of Geomagnetic Secular Variations

Ron Shaar <sup>1</sup> \*, Erez Hassul 1,2, Kate Raphael <sup>1</sup> , Yael Ebert <sup>1</sup> , Yael Segal 3,4, Ittai Eden<sup>1</sup> , Yoav Vaknin<sup>5</sup> , Shmuel Marco<sup>3</sup> , Norbert R. Nowaczyk <sup>6</sup> , Annick Chauvin<sup>7</sup> and Amotz Agnon1,2

<sup>1</sup> The Institute of Earth Sciences, The Hebrew University of Jerusalem, Jerusalem, Israel, <sup>2</sup> Neev Center for Geoinfomatics, The Hebrew University of Jerusalem, Jerusalem, Israel, <sup>3</sup> Department of Geophysics, Tel Aviv University, Tel Aviv, Israel, <sup>4</sup> Department of Marine Chemistry, Israel Oceanographic and Limnological Research, National Institute of Oceanography, Haifa, Israel, <sup>5</sup> Department of Archaeology and Ancient Near Eastern Cultures, Tel Aviv University, Tel Aviv, Israel, <sup>6</sup> Helmholtz Centre Potsdam, GFZ German Research Centre for Geosciences, Potsdam, Germany, <sup>7</sup> Univ Rennes, CNRS, Géosciences Rennes, UMR 6118, Rennes, France

#### Edited by:

Yohan Guyodo, UMR7590 Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie (IMPMC), France

#### Reviewed by:

Anita Di Chiara, Lancaster University, United Kingdom Chuang Xuan, University of Southampton, United Kingdom

> \*Correspondence: Ron Shaar ron.shaar@mail.huji.ac.il

#### Specialty section:

This article was submitted to Geomagnetism and Paleomagnetism, a section of the journal Frontiers in Earth Science

> Received: 24 May 2018 Accepted: 28 September 2018 Published: 25 October 2018

#### Citation:

Shaar R, Hassul E, Raphael K, Ebert Y, Segal Y, Eden I, Vaknin Y, Marco S, Nowaczyk NR, Chauvin A and Agnon A (2018) The First Catalog of Archaeomagnetic Directions From Israel With 4,000 Years of Geomagnetic Secular Variations. Front. Earth Sci. 6:164. doi: 10.3389/feart.2018.00164 The large and well-studied archaeological record of Israel offers a unique opportunity for collecting high resolution archaeomagnetic data from the past several millennia. Here, we initiate the first catalog of archaeomagnetic directions from Israel, with data covering the past four millennia. The catalog consists of 76 directions, of which 47 fulfill quality selection criteria with Fisher precision parameter (k) ≥ 60, 95% cone of confidence (α95) < 6 ◦ and number of specimens per site (n) ≥ 8. The new catalog complements our published paleointensity data from the Levant and enables testing the hypothesis of a regional geomagnetic anomaly in the Levant during the Iron Age proposed by Shaar et al. (2016, 2017). Most of the archaeomagnetic directions show <15◦ angular deviations from an axial dipole field. However, we observe in the tenth and ninth century BCE short intervals with field directions that are 19◦ -22◦ different from an axial dipole field and inclinations that are 20◦ -22◦ steeper than an axial dipole field. The beginning of the first millennium BCE is also characterized with fast secular variation rates. The new catalog provides additional support to the Levantine Iron Age Anomaly hypothesis.

Keywords: paleomagnetism, archaeomagnetism, Israel, levantine Iron-Age Anomaly, geomagnetic field, geomagnetic secular variations

#### INTRODUCTION

Despite decades of intense paleomagnetic research, many details of geomagnetic secular variations have still remained elusive. It is well accepted that secular variations average out globally to an axial dipole field over long geological timescales. Yet, many aspects concerning the spatial and temporal characteristics of secular variations remain unclear, especially when dealing with periods preceding direct human observational data. For example, while regional deviations of field direction from an axial dipole field are widely recognized, neither the degree limits nor the lifetime of these deviations are fully known. Today, the largest deviation from an axial dipole field, between 20◦ and 30◦ , occupies a confined area in the southern Atlantic associated with a low field intensity anomaly termed "South Atlantic Anomaly" (SAA) (Thebault et al., 2015). The question whether the SAA is a typical secular variations characteristic or, instead, a unique geomagnetic phenomenon is yet to be tested. Equivalently, it is not fully understood if rates and amplitudes of secular variations measured during the past few centuries (Jackson et al., 2000) also represent the characteristic behavior of the geomagnetic field in earlier periods. To fill these gaps in knowledge there is a growing need for reliable and precise paleomagnetic datasets in sub-millennial temporal resolution from periods preceding direct measurements of the geomagnetic field.

Archaeomagnetic data from in-situ archaeological structures, such as ovens, furnaces, kilns, and burnt buildings provide an excellent opportunity to capture the direction of the ancient field. When these structures cooled from high temperatures they acquired thermoremanent magnetization (TRM) parallel to the ambient field, thus preserving an instantaneous recording of the ancient field. In many cases, the age of the TRM can be precisely dated using radiocarbon, historical constraints, archaeological correlation, indicative pottery, coins, or a combination of these methods. In this perspective, the long, continuous, well-studied archaeological record of Israel offers a unique opportunity for archaeomagnetic research.

The global role of archaeomagnetic data from Israel is illustrated in **Figure 1A**, which shows a map of the published archaeomagnetic directional data available in the GEOMAGIA50 database (Korhonen et al., 2008; Brown et al., 2015) from the past four millennia. Our study area is located in an important geographic area that extends the densely scattered data from Europe to the southeast. To date, only several archaeomagnetic directions from Israel were published in journal articles (Aitken and Hawley, 1967; Segal et al., 2003; Shaar et al., 2016; Shahack-Gross et al., 2018). However, there are considerable unpublished data that are available only as unpublished theses (Segal, 2003; Hassul, 2015). The purpose of this work is to gather and compile all the available data from these sources and provide the first catalog of archaeomagnetic directions from Israel. To this end, we have collected the archaeological and the chronological information from the above sources, added new data from 15 additional sites, gathered all the raw paleomagnetic measurement data (if it exists) and translated them to a community standard MagIC format (Tauxe et al., 2016). The combined data were then re-analyzed using identical standards and selection criteria. The resulting catalog includes new secular variations data spanning the past four millennia.

#### METHODS

#### Sites and Locations

The Israeli archaeomagnetic catalog is assembled from a collection of several sources: Two unpublished Masters theses: (A) Segal (2003) that also includes two sites published in Segal et al. (2003) and (B) Hassul (2015); Two published articles: (C) Shaar et al. (2016), and (D) Shahack-Gross et al. (2018); and (E) New data from 15 structures labeled hereafter "this study". Here, we also revise and augment the directions previously reported in Shaar et al. (2016) with new measurements. Therefore, the paleomagnetic interpretations reported here are slightly different from those previously published and replace the previous interpretations. In the catalog we follow the standard paleomagnetic hierarchy nomenclature and define "location" as a collection of structures from the same area (i.e., an archaeological site), "site" as a single structure (i.e., cooling unit), "sample" as an oriented piece from a given site, and "specimen" as the part from the sample that was measured.

Sites can be classified under one of the three categories shown in **Figure 2**. Cooking ovens (tabuns, **Figure 2A**) are rounded structures, typically about 1 m in diameter, frequently found in domestic settings. Although the ovens used fire as a heating source, burnt remains that could be used for radiocarbon dating usually did not survive. Hence, the age of the ovens are typically dated by the age of living stratum in which it was found. Burnt walls (**Figure 2B**) are in-situ remains of large mudbrick structures, which were burnt as a whole during historical destruction events (e.g., Shahack-Gross et al., 2018). In the case of a large conflagration there may be a large amount of burnt organic material that can be directly dated and crosscorrelated with known historical military campaigns, leading to high precision dating, sometimes with an uncertainty of several years (e.g., Tel Megiddo; Finkelstein and Piasetzky, 2009), Tel Hazor (Sandhaus, 2013; Zuckerman, 2013), Tel 'Eton (Faust, 2008), Bethsaida (Arav, 2014), and Lachish (Ussishkin, 1990)]. Furnaces and kilns (**Figure 2C**) are large industrial structures that were used to manufacture ceramics. The kilns can be dated from the type of the ceramics and the typology of other finds, such as coins.

The **Supplementary Material** provides a short summary of the archaeological context and the location of each of the structures. We note that some sites were collected from old and presently inactive excavations, where a revised inspection of the exact archaeological contexts and the corresponding ages is not always possible. Also, in some previously studied sites there might be an ambiguity regarding their precise age because the dates of the corresponding archaeological strata have been refined throughout the years as more archaeo-chronological data have been accumulated. Therefore, we clearly stress the need for a detailed review of the archaeological contexts and the ages in this catalog. This effort requires a thorough and detailed archaeological study that is beyond the scope of this research, but we make such a future investigation possible with the information given in the **Supplementary Material**.

#### Sampling and Lab Procedures

Oriented samples were obtained either by drilling standard paleomagnetic cores, 1" in diameter, with a portable electrical drill or by hand samples. In the case of hand samples, orientations of flat surfaces were measured and marked in the field before detaching the oriented samples from the structure (site). In some cases, the material was hardened in the field with epoxy before being measured and removed. Specimens were prepared from the hand samples by sawing small cubes from the sample and gluing them inside non-magnetic paleomagnetic plastic sampling boxes. Samples were oriented in-situ with a Brunton compass prior to detachment and collection for both the cores (using Pomeroy or ASC orientation device) and the hand samples. A declination correction was added to the azimuth measurements

for all sites, except those of Segal (2003) from which we do not have the original measurement data. The latter results in a possible declination offset of the Segal (2003) dataset by 2–3 degrees.

The sampling, the formation, and the condition of the structures (sites) have a considerable effect on the uncertainty of the paleomagnetic directions. Man-made archaeological structures can collapse, break apart, incline or tilt, especially when the archaeological layer had been buried under a heavy overburden before being excavated. Therefore, extra care should be taken during paleomagnetic sampling of archaeological objects. To enable comparison between sites in the catalog, we assign to each site a "site formation quality index" (Qi) with fourscale grading indices (**Table 1**). The highest score (Qi = 1) is granted when the entire periphery of an oven or kiln was sampled or when several bricks from at least two walls in a burnt structure were sampled. A medium score (Qi = 2) is given when only a segment of a wall, oven, or a furnace was sampled. A low score (Qi = 3) is given to samples with high orientation uncertainty, and Qi=4 indicates poorly oriented samples.

#### Paleomagnetic Measurements and Data Analysis

The sample set of Segal (2003) was measured in the paleomagnetic lab in the Geophysical Institute of Israel using a 2G cryogenic magnetometer, and in the paleomagnetic laboratory at the University of Rennes, France, using a spinner magnetometer and a Leti cryogenic magnetometer. Part of the sample set of Hassul (2015) was measured in the paleomagnetic laboratory at the Helmholtz Centre Potsdam, GFZ Germany using a 2G cryogenic magnetometer. The majority of the data were measured in the paleomagnetic laboratory at the Institute of Earth Sciences, the Hebrew University of Jerusalem using a 2G cryogenic magnetometer. Specimens from all sites underwent progressive demagnetizations with Alternating Field (AF).

All the raw measurement data, except the dataset of Segal (2003) were translated into the community standard MagIC format (Tauxe et al., 2016) and merged into a single measurement file. All the data with the exception of Segal (2003) were re-analyzed, including previously published data, using the Demag GUI program, which is part of the PmagPy software package (Tauxe et al., 2016). Paleomagnetic directions of specimens and site means were calculated using the principal component analysis technique (Kirschvink, 1980) and Fisher statistics (Fisher, 1953). The interpretations follow a fairly strict set of selection criteria listed in **Table 2**, accepting only specimens with MAD (Kirschvink, 1980) ≤ 5, DANG (Tauxe and Staudigel, 2004) ≤ 5, and sites with n (number of specimen per site) ≥ 8, k (Fisher, 1953) ≥60, and α<sup>95</sup> ≤ 6. The measurement data and the interpretations (except Segal's dataset) are available in the MagIC database (https://www2.earthref.org/ MagIC).

symbols show data from Segal (2003). Locations of sedimentary cores data shown in Figure 4 are shown with asterisks with the same color code as in Figure 4.

FIGURE 2 | Types of archaeomagnegic structures used in this study: (A) Cooking ovens (tabuns) from Herodium (left) and Abel Beth Maacah (right). (B) Burnt structures from Tel 'Eton (left) and Tel Hazor (right). (C) Ceramic kilns from Yavneh. For a figures and archaeological information of all sites see Supplementary Material.

TABLE 1 | Site formation Quality Index.


RESULTS AND DISCUSSION

The catalog consists of 76 paleomagnetic sites collected from 33 different locations (archaeological sites) in Israel. Forty seven sites passing the acceptance criteria in **Table 2** are listed in **Table 3**, and their locations are shown in **Figure 1B**. The declinations and the inclinations of the accepted sites are shown in **Figures 3A,B**, where the inclination error bar is the α<sup>95</sup> value and the declination error bar is calculated using the equation: △ D = sin−<sup>1</sup> sin <sup>α</sup><sup>95</sup> cosI , where D is the declination error and I is the inclination. Sites that do not pass the criteria in **Table 2** are listed TABLE 2 | Acceptance criteria.


<sup>a</sup>Number in brackets refer to the following references: [1]: Kirschvink (1980); [2] Tauxe and Staudigel (2004); [3] Fisher (1953).

in **Table 4**. In the catalog we distinguish between sites that have all their measurement data available in the MagIC database and can be downloaded and re-interpreted by others using different criteria (labeled "this study" in **Tables 3**, **4** and marked with filled symbols in **Figures 1**, **3**) and sites of which we have only the site's mean parameters (labeled "Segal, 2003" in **Tables 3**, **4** and marked with open symbols in **Figures 1**, **3**).

From **Figures 3A,B** it can be seen that some periods have several coeval sites that show overlapping directions and demonstrate internal consistency and cross correlation between archaeological locations. These include: the fourteenth century BCE with two sites (Tel Megiddo and Tel Rehov); the thirteenth century BCE with two sites (Tel es-Safi and Tel Hazor); the beginning of the eighth century BCE with two sites (Tel Azekah, Tel Hazor,) the end of the eighth century with two sites (Bethsaida, Tel 'Eton), the third century BCE with two sites (Bnei Brak), and the second century BCE with two sites (Tel Shimron). An exceptional period with inconsistent directions is the tenth to ninth century BCE, which includes 11 sites from Tel Megiddo and Tel es-Safi, with non-overlapping directions. We interpret the latter as a time interval with fast changes in the geomagnetic field, and discuss this result in section Non-axial Dipole Field During the Iron Age Anomaly below.

#### Comparison With Global Models

The continuous curves in **Figure 3** show the predictions of three global spherical harmonic models of the geomagnetic field that use archaeomagnetic data as a data source: ARCH10K.1 (Constable et al., 2016), pfm9k.1b (Nilsson et al., 2014), and SHA.DIF.14k (Pavon-Carrasco et al., 2014). The pfm9k.1b model that is largely based on sedimentary data is smoother than the other two models. To first order, the archaeomagenticbased models nicely predict the trends in the direction of the Levant's geomagnetic field, but with a lower amplitude. During the past 2.5 millennia only two time intervals do not fit the models: the first century BCE and the seventh century CE. In earlier periods several misfits are observed in pfm9k.1b and SHA.DIF.14k between the eighteenth century BCE and seventh century BCE, especially in periods with high positive declination and high (steep) inclination. Model ARCH10K.1 show the best fit to our data. It is likely that with the new data in the current catalog, misfits will be minimized in future geomagnetic models that will be refined by the new data.

#### Comparison With Sedimentary Data

There are several advantages of archaeomagnetism over sedimentary magnetic data. The archaeomagnetic TRM does not suffer from inclination shallowing, lock-in depth, and post-depositional effects associated with depositional remanent magnetization (DRM). Hence, paleomagnetic directions from archaeomagnetism can sometimes provide better precision and age control than sedimentary data. However, sedimentary magnetism provides continuous datasets spanning much larger time intervals than archaeomagnetism. Both types of records are available in Israel and we compare them in **Figure 4** that shows data from four Holocene cores available in the GEOMAGIA50 database (Brown et al., 2015). **Figures 4A–C** show data from three piston cores raised from the Holocene Dead Sea (Frank et al., 2007a,b). These cores were obtained without azimuthal orientation and Frank et al. (2007a,b) corrected their declination profile by setting the mean declination value of the core top to zero. Their age model is based on a large number of radiocarbon measurements (Migowski et al., 2004). Given the uncertainty in the age models of the Dead Sea cores, there is fairly good agreement between the archaeomagnetic and the sedimentary data, especially for the past two millennia. The high inclination values observed in the ninth century BCE are not observed in the sedimentary data. This could be a result of inclination shallowing, post depositional magnetization, and smoothing of the sedimentary data. As the Holocene Dead Sea sediments are dominated by authigenic greigite (Ron et al., 2006; Frank et al., 2007b; Thomas et al., 2016; Ebert et al., 2018) a complicated magnetic acquisition mechanism is expected. Thus, differences between the archaeomagnetic and the sedimentary data are likely. **Figure 4D** shows data from a core taken in the Birket Ram maar lake in northern Israel, which was dated using only two radiocarbon ages (Frank et al., 2002). Here, the trends in the inclination and the declination profiles agree with the archaeomagnetic data, but the temporal resolution in Birket-Ram is much lower than at the Dead Sea. In summary, we observe a reasonable correlation between the sedimentary and the archaeomagnetic data highlighting the potential of combining these two types of records into a single joint master secular variation curve for the Levant. Yet, owing to the large uncertainties in both the sedimentary magnetic acquisition mechanism and the sedimentary age models this challenge requires a more detailed investigation.

#### Non-axial Dipole Field During the Iron Age Anomaly

Shaar et al. (2016, 2017) hypothesized a positive local geomagnetic anomaly in the Levant between 1050 BCE to 750 BCE which they termed "The Levantine Iron Age Anomaly" (LIAA). The paleointensity data from the Levant supporting the LIAA hypothesis (**Figure 5A**) show high field values between the mid-Eleventh century BCE and the eighth century BCE and two geomagnetic spikes (virtual axial dipole moment, VADM > 160 ZAm<sup>2</sup> ) (Ben-Yosef et al., 2009, 2017; Shaar et al., 2011, 2016). Fourteen sites in **Figure 3**, from Tel Megiddo, Tel es-Safi, Tel Azekah, Tel Hazor, and Bethsaida cover the interval between 900 and 750 BCE. Of these sites, 12 show inclinations above 60◦ , while two sites from Tel Megiddo show exceptionally high inclinations of 73◦ (mgq05t2, Shaar et al., 2016) and 71◦ (QTMB; Hassul, 2015; this study), considerably higher than the expected geocentric axial dipole (GAD) inclination in Jerusalem (51◦ ). We note that one site in Segal (2003) dataset (Ceramic kiln from Kfar Menachem, see **Supplementary Material**) showed an even higher inclination of 81◦ around this time. Yet, this site did not pass the selection criterion for the age as the age of the kiln was not supported by any direct dating method, only by a correlation with nearby archaeological sites. The declinations during the Levantine Iron Age Anomaly interval show large scatter in the ninth century with declinations ranging from −3 ◦ to 23◦ , and closely grouped values around 5◦ in the eighth century. The angle between the archaeomagnetic directions and the axial dipole field (**Figure 3C**) show exceptionally high values between 19◦ and 22◦ around the ninth century in Tel Megiddo. These are from two sites from Shaar et al. (2016) (mgq05t1, mgq05t2) and from two sites from this study (QTMA, QTMB). Today, angular deviations from the axial dipole exceeding 20◦ occur only in the southern


TABLE

3


Paleomagnetic

 means of sites passing acceptance

 criteria.

(Continued)



hemisphere in the area affected by the South Atlantic Anomaly (Thebault et al., 2015; Figure 7 in Shaar et al., 2016).

To inspect how irregular is the angular deviation from axial dipole field during the Iron Age Anomaly, we gathered all the global archaeoamgnetic and volcanic data from the past four millennia from the GEOMAGIA50 database (shown in **Figure 1A**). We applied the same acceptance criteria to the global data, and plot on **Figure 6** the angular deviation from axial dipole. Only 1% show angular deviations larger than 19◦ , demonstrating that the high inclination high declination episode in the tenth–ninth century BCE is a unique phenomenon. Altogether the new archaeomagnetic data covering the LIAA show large azimuthal deviations from an axial dipole (**Figures 3C**, **6**) and steep inclinations (**Figures 3B**, **5B**), corroborating the LIAA hypothesis of Shaar et al. (2016, 2017).

#### Archaeo-Chronological Applications

The forty seven archaeomagnetic data points reported here comprise a significant first step toward a master archaeomagnetic curve for the Levant, from which archaeomagnetic dating can be developed, using an approach similar to that described in Lanos (2004) and Pavon-Carrasco et al. (2014). Still, even in the absence

of a continuous secular variation curve, archaeo-chronological insights can be obtained by merely comparing archaeomagnetic data from structure whose age is not tightly constrained, with the available data. The burnt structure MGDF from Shahack-Gross et al. (2018) is an example of this approach. For dating the destruction that caused the fire at this archaeological level, Shahack-Gross et al. (2018) compared the paleomagnetic direction of MGDF with the available archaeomagnetic directions of periods with known destructions in Megiddo (Late Bronze III, Late Iron I, Iron IIA, and Iron IIB) and found that only one period fit to the archaeomagnetic data. Another promising potential of archaeomagnetism is the opportunity to use archaeomagnetic time-series of fast variation in the geomagnetic field, such as the Iron Age, to complement radiocarbon data, in particular for periods when the radiocarbon calibration curve is flat. For a more detailed review on addressing the Iron Age chronology problems with archaeomagnetism see Stillinger et al. (2016, 2018) and Herve and Lanos (2017).

#### SUMMARY AND CONCLUSIONS

We initiate here the first catalog of archaeomagnetic directions from Israel, with 76 sites that have been collected and analyzed since 2003. From this catalog, 47 sites pass a set of strict selection criteria and serve as a basis for future development of a master archaeomagnetic secular variation curve for the Levant. The measurement data of this catalog can be approached from MagIC database (https://www2.earthref. org/MagIC) and the archaeological information is provided in the **Supplementary Material** of this article. Some misfits with the global models prior to the second half of the first millennium BCE imply that more archaeomagnetic data are needed to refine our knowledge of fast secular variations during early periods in the Holocene. The most prominent feature in the new catalog is the high inclinations and high declination interval in the ninth century BCE. During this interval two sites, one published here, and another published in Shaar et al.

(2016), show steep inclinations exceeding 70◦ , while the axial dipole inclination in Jerusalem area is 51◦ . Four other sites, two published here, and two published in Shaar et al. (2016), show angular deviations from an axial dipole between 19◦ and 22◦ . Only <1% of the published archaeomagnetic data from the past four millennia show angular deviations exceeding 19◦ suggesting an unusual field behavior. The large angular deviation occurred during a period with extremely high field intensity in the Levant, providing additional support to the Levantine Iron Age Anomaly hypothesis of Shaar et al. (2016, 2017) of a regional high field anomaly at the beginning of the first millennium BCE.

#### AUTHOR CONTRIBUTIONS

EH, RS, KR, YE, YS, YV, SM, AA, and IE collected the samples. EH, RS, YE, YS, YV, SM, NN, IE, and AC made or conducted the measurements. All authors contributed to data analyses and manuscript preparation. RS compiled the data, made the final interpretations, figures and tables.

#### FUNDING

This study was supported by Israel Science Foundation grant 1364/15 to RS and 1181/12 to AA, German-Israeli Foundation (GIF) Young Scientists program grant I-2398-301.8/2015 to RS and Ministry of National Infrastructure, Energy, and Water Resources grant 215-17-027 to RS. Work by IE and YE was supported by Alpha–Research Program in the Sciences of the Future Scientists, Center for the Advancement of the Gifted and Talented in Israel (https://www.madaney.net/en/homepage).

#### ACKNOWLEDGMENTS

We kindly thank all the archaeologists who helped in the sampling process and provided their most appreciated assistance: Israel Finkelstein (Tel Megiddo), Ruth Shahack Gross (Tel Megiddo), Amihai Mazar (Tel Rehov), Uri Davidovich (Tel Rehov), Aren Maeier (Tel es-Safi), Oded Lipschitss (Tel Azekah), Sharon Zuckerman (Tel Hazor), Shlomit Bechar (Tel Hazor), Amnon Ben-Tor (Tel Hazor), Rami Arav (Bethsaida), Eli Yannai (Yavne1), Liat Ziv (Yavne2), Ron Beeri (Bnei Brak), Avraham Faust (Tel Eton), Arthur Segal and Michael Eisenberg (Susita-Hippos), Nava Panitz-Cohen (Abel Beth Maacah), Naama Yahalom-Mack (Abel Beth Maacah), Oren Tal (Apollonia), Gideon Hadas (Ein Gedi), Doron Ben-Ami (Givati), Boaz Zissu (Jamjum), Alla Nagorski (Qatra), Itzick Shai (Tel Burna), Mario Martin (Tel Shimron). We wish to thank Yotam Asscher from the Israeli Antiquity Authority (IAA) for his help. We acknowledge the long standing support of

#### REFERENCES


the Israeli Antiquity Authority in the Israeli archaeomagnetic project. We thank Prof. John Hall for his review of the final manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2018.00164/full#supplementary-material


Levantine Iron Age chronology. Near East Archaeol. 81, 141–144. doi: 10.5615/neareastarch.81.2.0141


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Shaar, Hassul, Raphael, Ebert, Segal, Eden, Vaknin, Marco, Nowaczyk, Chauvin and Agnon. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Archeomagnetic Intensity Spikes: Global or Regional Geomagnetic Field Features?

Monika Korte<sup>1</sup> \* and Catherine G. Constable<sup>2</sup>

*<sup>1</sup> GFZ German Research Center for Geosciences, Potsdam, Germany , <sup>2</sup> Scripps Institution of Oceanography, University of California, San Diego, San Diego, CA, United States*

Variations of the geomagnetic field prior to direct observations are inferred from archeo- and paleomagnetic experiments. Seemingly unusual variations not seen in the present-day and historical field are of particular interest to constrain the full range of core dynamics. Recently, archeomagnetic intensity spikes, characterized by very high field values that appear to be associated with rapid secular variation rates, have been reported from several parts of the world. They were first noted in data from the Levant at around 900 BCE. A recent re-assessment of previous and new Levantine data, involving a rigorous quality assessment, interprets the observations as an extreme local geomagnetic high with at least two intensity spikes between the 11th and 8th centuries BCE. Subsequent reports of similar features from Asia, the Canary Islands and Texas raise the question of whether such features might be common occurrences, or whether they might even be part of a global magnetic field feature. Here we use spherical harmonic modeling to test two hypotheses: firstly, whether the Levantine and other potential spikes might be associated with higher dipole field intensity than shown by existing global field models around 1,000 BCE, and secondly, whether the observations from different parts of the world are compatible with a westward drifting intense flux patch. Our results suggest that the spikes originate from intense flux patches growing and decaying mostly in situ, combined with stronger and more variable dipole moment than shown by previous global field models. Axial dipole variations no more than 60% higher than observed in the present field, probably within the range of normal geodynamo behavior, seem sufficient to explain the observations.

#### Edited by:

*Greig A. Paterson, Institute of Geology and Geophysics (CAS), China*

#### Reviewed by:

*Chuang Xuan, University of Southampton, United Kingdom Philip Livermore, University of Leeds, United Kingdom*

> \*Correspondence: *Monika Korte monika@gfz-potsdam.de*

#### Specialty section:

*This article was submitted to Geomagnetism and Paleomagnetism, a section of the journal Frontiers in Earth Science*

> Received: *07 November 2017* Accepted: *16 February 2018* Published: *07 March 2018*

#### Citation:

*Korte M and Constable CG (2018) Archeomagnetic Intensity Spikes: Global or Regional Geomagnetic Field Features? Front. Earth Sci. 6:17. doi: 10.3389/feart.2018.00017* Keywords: archeomagnetism, paleomagnetism, intensity spike, secular variation, global magnetic field models

#### 1. INTRODUCTION

Variations of the geomagnetic field prior to direct observations can be inferred from archeo- and paleomagnetic experiments. Investigations of past secular variation are relevant to understanding the geodynamo process in Earth's core. Seemingly unusual variations not seen in the present-day and historical field are of particular interest to constrain the full range of core dynamics. These include "geomagnetic intensity spikes," characterized by very high field intensity values that are associated with rapid secular variation rates. They were first noted in the Levantine area by Ben-Yosef et al. (2009) and Shaar et al. (2011) at around 900 BCE. Subsequent studies from Turkey (Ertepinar et al., 2012) and Georgia (Shaar et al., 2013) confirm high regional intensities. Based on a careful quality assessment of new and previously published Levantine data Shaar et al. (2016) interpret the observations as an extreme regional geomagnetic high with at least two intensity spikes between the 11th and 8th centuries BCE.

Reported high intensity values from other parts of the world at similar times, however, raise the question whether similar features occurred more commonly or if they are even part of a global magnetic field feature. No clear definition of intensity spikes exists. High intensity values potentially related to spikes have recently been reported from the Canary Islands (de Groot et al., 2015; Kissel et al., 2015), Korea (Hong et al., 2013), eastern China (Cai et al., 2016), and Texas, United States of America (Bourne et al., 2016).

Livermore et al. (2014) investigated how core flow models might explain a localized spike and conclude that the reported rates of secular variation from the Levant are not compatible with commonly accepted core-flow dynamics under the assumption of no radial magnetic diffusion. However, the careful reinterpretation of the Levantine results by Shaar et al. (2016) revised the local maximum secular variation rates to slightly lower values as it eliminated some of the highest intensity values reported previously. Davies and Constable (2017) recently tested the hypothesis that the Levantine spike could be caused by a narrow localized intense flux patch at the core-mantle boundary (CMB). They conclude that the spike signature at the surface must span at least 60◦ in longitude and suggest that the source is a flux patch first growing in place and then migrating northand westward.

Numerous global spherical harmonic based reconstructions of the geomagnetic field evolution through the past millennia have been developed over the last decade (see Constable and Korte, 2015; Korte et al., 2017, for reviews). None of them shows sharp intensity spikes. This is not surprising as both spatial and temporal resolution of such models are clearly limited compared to what is known from the present-day field, due to heterogeneous, sparse global data distribution and uncertainties in data and dating. Moreover, many of the recently reported high intensity values have not yet been included in global models.

Here we use spherical harmonic modeling to test the hypotheses that the Levantine and other potential spikes might be related to dipole field intensity around 1,000 BCE that is higher and more variable than previously thought, and whether the observations from different parts of the world are compatible with a westward drifting intense flux patch.

In section 2 we first provide an overview of recently published high intensity archeomagnetic data discussed in the context of geomagnetic spikes and other published data used in this study. We compare secular variation and potential westward drift rates around the reported spikes to maximum rates observed in high resolution field models based on recent magnetic observatory and satellite observations in section 3. We then (section 4) derive a suite of new global models to test if what we consider a physically reasonable model can explain the spike observations by invoking a higher dipole moment than seen in previous models and/or westward drift of an intense flux patch around 1,000 BCE, without degrading the fit to the other available archeomagnetic observations. We discuss our findings in section 5 before we briefly summarize our conclusions.

#### 2. DATA AND COMPARISON OF REPORTED GEOMAGNETIC SPIKE OBSERVATIONS

Three compilations of data are considered in this study:


The term "geomagnetic intensity spike" was coined by Ben-Yosef et al. (2009), who reported among other data two very high archeointensity results that yield values higher than 200 ZAm<sup>2</sup> in terms of virtual axial dipole moment (VADM), and indicating very rapid field changes. Recently, Shaar et al. (2016) re-interpreted all data from the studies by Ben-Yosef et al. (2008a,b, 2009) and Shaar et al. (2011) using an automatic, and hence more objective, interpretation of the experimental results from the PmagPy Thellier GUI software by Shaar and Tauxe (2013) and applying very strict quality control in the acceptance criteria. They also find several high intensity values, but none exceeding 200 ZAm<sup>2</sup> . We consider only the revised Levantine data set by Shaar et al. (2016) and deleted the values from the Ben-Yosef et al. (2008a,b, 2009) and Shaar et al. (2011) studies from the global data set taken from GEOMAGIA50.v3. More recent work by Ben-Yosef et al. (2017), which had not yet been published when we carried out our study, further corroborates the Levantine data set without exceeding intensity amplitudes found by Shaar et al. (2016).

All recent studies forming our spike data set are listed in **Table 1**. Except for the work from Ertepinar et al. (2012) these studies had not been included in GEOMAGIA50.v3 in October 2016 and the data have been taken directly from the publications. Most of the data derive from archaeological material, with the exceptions of volcanic rocks in case of the Canary Islands and relative intensity results from cave sediments from Texas.

A limited number of sediment records have been selected from the compilation by Korte et al. (2011) to improve global data coverage. Selection criteria, in particular in areas with multiple available records, were availability of the full vector information and temporal resolution or sedimentation rate. Relative paleointensity (RPI) records do not always reliably reflect the geomagnetic signal (see e.g., Tauxe and Yamazaki, 2015), so as an additional simple criterion to ensure RPI reliability we only chose records, reflecting the general dipole moment trend with higher values between 1,000 BCE and 1,000 CE than at earlier or later times, which is clearly reflected in intensity variations on multi-millennial periods (see e.g., Figure 20 of Constable and Korte, 2015). These records, listed in **Table 2**, have been used with the uncertainty estimates provided by Panovska et al. (2015) following the method given by Panovska et al. (2012). A recently published new record from Lake Mavora, New Zealand (Turner et al., 2015) has been added. Uncertainty estimates for both this and the Texan cave sediment record included in our spike data set have also been determined following Panovska et al. (2012). **Figure 1** shows the locations of the data from all three sets from the time interval 2,000–0 BCE. All spike data come from a latitudinal band between 28 and 42◦N.

To obtain an overview of potential spikes observed in the different regions we compare the recently published spike data among the different regions and to the other archeomagnetic data from nearby (**Figure 2**). The other data come from approximately 1,500 km radii around the spike data. The best documented spikes at present clearly are the two Levantine spikes at about 980 and 740 BCE. The Chinese spike appears

TABLE 1 | Highest intensity values potentially related to geomagnetic spikes, converted to VADM, and their ages from East to West.


\**Scaled relative paleointensity value from cave sediments, see text for details.*

TABLE 2 | Sediment records used in this study with their average sedimentation rate (SR).


considerably older, but its age uncertainties are significant and partly overlap with the 980 BCE Levantine spike data. Similarly, the Canary Island spike appears clearly younger, but age uncertainties extend to the times of the 800 BCE Levantine spike.

The Chinese and Canary Island data sets most clearly show how age uncertainties hamper the detailed investigation of high intensity spikes: Several data are only known to lie within the same age range and are assigned the same central age. Small data uncertainties may suggest that the data come from different ages within the age interval, but any time series or spatial interpolation will tend to consider them at the same time and average the values, not even trying to fit the maximum values closely.

The few new Korean data around 1,200 BCE are high, but do not clearly characterize a spike. The Texan cave sediments, on the other hand, give an exceptionally high value at about 890 BCE, which exceeds 400 ZAm<sup>2</sup> when the whole RPI time series is scaled by a global model. The feature shows more than double the intensity seen elsewhere, yet is supported by only three data points so far. It should be regarded with caution until confirmed by future data. Previous absolute intensity data from the region are sparse.

The relatively abundant other data around the Levant, Turkey and Georgia indicate another factor complicating the derivation of a clear spatio-temporal picture of the spike: the scatter in regional data is rather large (the occasionally also large age uncertainties are not even shown in **Figure 2**) with uncertainties not always overlapping, at least at one sigma level. Older studies did not consider tests and quality criteria developed only recently, yet rejecting all data without this information depletes the available information unreasonably. Most previous data are below the spike maxima, but a few individual very high intensity data are noted at different times, in the examples shown here in particular in previously published data from China. As shown in the work by Shaar et al. (2016) such values might not survive critical re-evaluation, but it also clearly seems unwarranted to discard them without examination. Unfortunately the raw data required for a re-assessment as performed by Shaar et al. (2016) is hardly ever available for older data. For this study we take all published data at face value and do not attempt to judge their quality in any way.

Concerning regional spike observations we also note that high intensity values reaching up to 200 ZAm<sup>2</sup> and thus reminiscent of the presently discussed spikes have been found at 890 BC in Hawaiian lavas by Pressling et al. (2006, 2007).

Calibrated relative paleointensities from the sediment records considered in this study are included in **Figure 2**. Interestingly, both South American records (esc, tre) show a signature reminiscent of reported intensity spikes around 1,050 BCE, i.e., between the Asian and Levantine highs, and the Italian Lago di Mezzano (mez) record seems to reflect the two Levantine spikes.

The spatial distribution of all new spike reports permit different interpretations. Firstly, considering the age uncertainties there might have been only two spikes of large spatial extent between about 1,000 and 750 BCE. Secondly, there might have been at least four regionally confined spikes between about 1,400 and 500 BCE. Thirdly, considering that

high intensities are observed earlier in Korea and China than in the Levant and later on the Canary Islands, all at similar mid-latitudes, westward drift of magnetic flux structure might be involved. **Table 3** shows simple estimates of the order of core flow velocities required under the frozen flux assumption for this interpretation. The estimates were obtained by considering large-circle distances between the approximate centers of the regions, time differences between the highest intensity values as given in **Table 1** and considering that with westward drift the spike should have occurred earlier in Korea than China (well possible within the Chinese age uncertainties), and 3481 km as radius of the outer core to obtain velocities at the CMB. The reported Canary island spike is younger than the Texas one, so they are incompatible with a westward drift explanation of a single flux patch causing one of the Levantine spikes. With 0.35 and 0.11◦ /year or 21 and 7 km/year at the CMB, for a feature drifting from Asia to the Levant and from there to the Canary Islands, respectively, these rates are well within core flow velocities deduced from the present-day field (e.g., Hulot et al., 2002; Livermore et al., 2014) which can reach values in the order of 0.9◦ /year or 55 km/yr. Typical present-day values for the large-scale flow are lower, in the order of 15 to 20 km/year (e.g., Holme, 2015), but the high-latitude circumpolar jet deduced by Livermore et al. (2017) to explain the accelerating westward drift of intense flux lobes over Canada and Siberia observed from satellite data reaches maximum speeds of 40 km/year. Based on simple considerations of paired records only the Levantine and Texan spike observations seem incompatible with a westward drifting flux patch, and even that would only require core flow velocities 33% higher than inferred from high resolution present-day field models.

#### 3. INSIGHTS FROM GLOBAL GEOMAGNETIC FIELD MODELS

In order to see how unusual the spikes appear in the context of known secular variation we perform some comparisons with global archeomagnetic and high resolution present-day field models. Shaar et al. (2016) pointed out that present-day regional deviations from the axial dipole moment range from 50 to 150%, while the Levantine spike exceeded the axial dipole moment at that time by a factor of two. This is reflected by global minimum and maximum VADM values tabulated for the International Geomagnetic Reference Field (IGRF Thébault et al., 2015) for 2015 in **Table 4**. The standard deviation obtained from a global grid of VADM values is 15 ZAm<sup>2</sup> , but note that they are not normally distributed. Despite the fact that the lowest values are <50% and the highest ones do not reach 150% the globally averaged VADM is slightly higher than the (tilted) dipole moment (DM). We do not know the actual DM strength at the times of the archeomagnetic spikes. Despite displaying the same general trends global archeo- and paleomagnetic models show differences in more detailed structures, and recovering the absolute DM strength is surprisingly challenging as it critically depends on the, globally seen, rather sparse absolute intensity information (Panovska et al., 2015). At times of their highest DM, within two centuries from the reported spikes, the three Holocene models ARCH10k.1, CALS10k.2, and HFM.OL1.A1 (Constable et al., 2016) reach maximum local VADM values of up to 143 ZAm<sup>2</sup> . These models include some of Levantine spike data, but all of them clearly underestimate the associated high intensity values. The lower VADM standard deviation and closer similarity of globally averaged VADM and DM compared to IGRF 2015 are indicative of the lower spatial resolution of these models. Considering all available data VADMs for a 300 year interval around the spikes yields maximum values higher than the revised Levantine results (**Table 4**) (which come from Hawaii in 890 BCE; Pressling et al., 2006, 2007) and an average value clearly higher than the maximum model global averages. The averaged data VADM is closer to the archeomagnetic model average if the recently published spike data are not included. All this might suggest that the actual DM was higher during the occurrence time(s) of the spikes, which might make the magnitude, although not necessarily the rates

bars indicate one standard error in magnetic data and age uncertainties as provided in the studies. Age errors of sediments and other archeomagnetic/volcanic data have been omitted for clarity of figures. All sediment data were calibrated by the ARCH10k.1 global model (see text). The gray intervals in (A) mark the approximate times of the discussed spikes by highest intensity values and their age uncertainties. Note the different ordinate scale in the bottom left panel.

of change, reconcilable with present-day VADM deviation from (axial) DM.

Not considering the Texas cave sediment data, the fastest rates of intensity changes related to archeomagnetic spikes seem to come from the Levant and China. In the former region, the fastest change seems related to the spike around 740 BCE, with age uncertainty interval 800–732 BCE. Assuming that the data all are correct within their data and age uncertainties, these values suggest a minimum rate of intensity or VADM change of 31 µT/century or 61 ZAm<sup>2</sup> /century, respectively. This is clearly lower than the rates estimated by Livermore et al. (2014) based on the originally published Levantine data. The Chinese values around 1,300 BCE have a rather large age uncertainty of ±300 years, giving an interval of 1,600 to 1,000 BCE. However, assuming that a low value dated at 1,328 ± 28 BC is correct and the highest values occurred at 1,000 BC the minimum required rate of change, again also considering data uncertainties, amounts to 15 µT/century in intensity or 30 ZAm<sup>2</sup> /century. These rates are tabulated in **Table 5** together with several estimates based on models of the present and historical field, with gufm spanning the interval 1,590–1,990 based on historical and observatory geomagnetic data (Jackson et al., 2000). All values are given as change per century, but note that they have been calculated as first differences over the given time intervals, i.e., representing average changes over 5–400 years intervals. DM and axial dipole coefficient changes are not given for 400 years, because the axial dipole coefficient in gufm has been extrapolated back between 1,840 and 1,590 due to the lack of historical absolute intensity information. Nearly instantaneous and average changes are remarkably similar over the past ca. 200 years for field intensity and DM including similar rates of axial dipole change. Larger differences in VADM than F result from the latitude dependence in the conversion. The lower values of dF/dt over 400 years might indicate that steady monotonic rates of change did not persist over this time scale at any location on Earth. The Chinese field change is only slightly higher than and thus appears compatible with recent variation rates, but the Levantine field change indeed appears more than twice as fast.

#### 4. CAN GLOBAL FIELD MODELS REPRODUCE RAPID SPIKES?

Presently available global geomagnetic field models do not reproduce rapid intensity variations, at least partly because they are strongly smoothed in space and time as a consequence of limited data coverage and, compared to present day observations, large data and dating uncertainties. In order to test how closely the spike data can be explained by physically reasonable models we derived a large number of test models. The modifications among these models mainly include, individually or in combination, variations of the data basis and weighting of sub-sets of data to enforce a closer fit to the spikes, and background models favoring westward drift over growth and decay of flux concentrations. Note that all these models clearly are not considered to be the best representations of the global field, but are designed to test two hypotheses: firstly, might the spikes be related to higher and more variable axial dipole contributions

TABLE 3 | Required core flow velocities if spike observations were related to westward drifting features.


TABLE 4 | Global VADM statistics from models and data compilations.

Model or data set IGRF 2015 ARCH10k.1 860 BCE CALS10k.2 930 BCE HFM.OL1.A1 930 BCE All data 1,100–800 BCE No spike data 1,100–800 BCE DM (ZAm<sup>2</sup> ) 77.2 107.4 96.7 99.1 – – Avg. VADM (ZAm<sup>2</sup> ) 80.4 107.3 95.4 97.8 115.9 107.1 σ VADM (ZAm<sup>2</sup> ) 14.9 11.4 11.7 11.0 27.9 27.5 Min. VADM (ZAm<sup>2</sup> ) 39.7 86.0 73.4 78.8 36.3 36.3 Max. VADM (ZAm<sup>2</sup> ) 116.3 142.6 123.5 125.3 200.7 200.7

All our spherical harmonic models with time-dependent coefficients based on cubic B-splines follow the modeling strategy of the CALSxk and ARCHxk models (for details see, e.g., Korte et al., 2009; Panovska et al., 2015). They span the interval 4,000 BCE to 1,990 CE, but here we limit the analysis to 2,000– 0 BCE. Linearisation was done by 15 iterations without outlier rejection. Space and time regularizations minimizing Ohmic dissipation (Gubbins, 1975) and the second derivative of the radial field were applied (see e.g., Jackson et al., 2000; Korte et al., 2009). The strength of regularizations was chosen by comparison to present-day geomagnetic power spectra of main field and secular variation. The influence of this regularization is discussed in more detail in section 5. We present our findings based on seven representative cases out of a larger range of test models studied. Their different characteristics are summarized in **Table 6**.

TABLE 5 | Maximum rates of geomagnetic field intensity change in models compared to minimum estimates required to explain the observed spikes in China and the Levant.


*Change of field intensity expressed as F and VADM, global dipole moment (DM) and axial dipole coefficient g<sup>0</sup> 1 .*

#### TABLE 6 | Characteristics of test models.


*Data types a and v stand for archeomagnetic and volcanic, respectively. SH, Southern hemisphere; sed, sediments; For details on notes see text.*

Model A is a rather standard model based on archeomagnetic and volcanic data only, similar to e.g., SHA.DIF.14k (Pavón-Carrasco et al., 2014) or A\_FM (Licht et al., 2013) and in particular differing from ARCH10k.1 (Constable et al., 2016) only in the addition of the spike data, weighting of data (see below) and resolution.

Model C additionally includes the Hall's cave (Texas) record and, in order to improve the global data coverage, the southern hemisphere sediment records.

Model AS includes all data, but all sediment and other archeomagnetic and volcanic data are downweighted compared to the spike data, in this case to 0.67 for directions and only 0.27 for all intensities. The weighting scheme used here is one example from a number of different schemes that we tried of giving more weight to the spike data in order to push the model to primarily fit these data.

For model M all data used were weighted similarly to model AS, but the main spikes' data (Levant, China, Canary Islands) were modified within their age uncertainties to avoid the contemporaneity of high and low values and represent a possible temporal field evolution within the age error interval and in relation to temporally adjacent values. In detail, these changes, depicted in **Figure 3**, aimed for least rapid variations and comprise the following: (i) The 980 BCE Levantine spike was broadened as far as possible by shifting the 977 and 974 BCE values to their youngest and the lower values at 985 BCE to their oldest possible age within the given age uncertainties. All values from the second spike given as 740 BCE were ordered by intensity and distributed to form a uniform slope toward the lower boundary of the highest values within their age uncertainty interval, with lowest values at 800 and highest at 732 BCE. There is a data gap between 732 and 621 BCE, allowing for a similar or less steep slope of this spike on its younger side. (ii) For China, motivated by two low VADM values with small uncertainties at 1328 and 1,256 BCE and models like AS that suggest an intensity maximum closer to 1,000 than 1,300 BCE we shifted all the 1,300 BC values to form a uniform slope toward the highest

value at 1,000 BCE, the youngest end of the given age uncertainty interval. Again, a data gap between that time and 697 BCE allows for a similarly gentle slope of the spike toward the younger side. (iii) Five values were shifted within their age uncertainties in the Canary Island data set, namely from 667 to 867 BCE, from 575 to 450 BCE, from 579 to 440 BCE, from 615 to 420 BCE and from 590 BCE to 730 (see **Figure 3**).

data are not shown.

Model WW is used as a background model to test if the data are compatible with westward drifting flux features. It was designed using SH degrees and orders 1–4 from ARCH10k.1 at 980 BCE together with degrees and orders 5 to 10 from gufm at an arbitrary time (here 1,810 CE) for higher spatial resolution. This time-slice, showing a strong normal flux patch under the Levant area, was then rotated to give a westward drift of 0.32◦ /year over the whole model time, resulting in a full rotation of the flux pattern around the Earth in 1.125 kyr.

Models Cww and Mww use the same data and weighting as C and M, respectively, but the regularization is applied around background model WW. The regularization in these cases trades off between fit to the data and the westward drifting features shown by the WW model, i.e., if the westward drift is not preserved in the resulting model it is not compatible with (parts of) the data.

Models Mh and Mhww are higher resolution versions of models M and Mww. The only differences are relaxed regularization constraints, i.e., weaker regularization factors somewhat subjectively chosen to better fit the data, in particular the spike signatures, at the cost of containing more (secular variation) spectral power compared to the present-day field than the other models. In the case of Mhww this means that the imposed drift is stronger than in Mww.

Model Dip is a pure, constant axial dipole included here for comparison of deviation of the data from any of the models. Its strength of −42 µT was chosen empirically as giving the smallest root mean square (RMS) misfit for the archeomagnetic and volcanic data in the time interval 2,000–0 BCE.

**Table 7** lists the root mean square misfit of the models to the original data set, weighted by the data uncertainties. Note that here we take a different approach for weighting from previous models which were aimed at best reflecting the geomagnetic field evolution all over the globe. Instead of setting relatively large but in our opinion more realistic minimum data uncertainties, we here only consider uncertainties smaller than <sup>α</sup><sup>95</sup> <sup>=</sup> 1.7◦ (equal to 1 ◦ in inclination) and 1 µT as unrealistically small and use these values as minimum uncertainties. Missing uncertainty estimates were set to <sup>α</sup><sup>95</sup> <sup>=</sup> 3.4◦ and 5µT, values that commonly were used as minimum uncertainty estimates in previous models. For sediments we kept previous minimum and missing uncertainty assignments of α<sup>95</sup> = 6 ◦ and 5µT.

The choice of sediment records and in particular the size of their estimated uncertainties leads to this data type being fit nearly within their uncertainties in all components, while archeomagnetic and volcanic data show much larger normalized misfits. This does not mean that individual features of sediment records are fit particularly well, as can be seen in the two southern hemisphere records included in **Figure 4**. In particular intensity maxima around 1,200 BCE in South America (record esc) and slighly earlier in New Zealand (record mav), which might have supported a global source of the Asian spike, are not represented by the models. However, the sediment uncertainties in fact are such that these data also are nearly fit to a weighted rms of 1 by a constant axial dipole (**Table 7**). For archeomagnetic and volcanic data all models, except the ones with strong enforced westward drift (WW, Mhww), show smaller misfits than a constant dipole e.g., ranging from 3.60 to 4.74 compared to 6.14 for all components, or 3.74 to 4.37 compared to 7.10 for intensity. The westward drifting background model WW is less compatible with all data than a constant axial dipole. If we allow more spatial complexity and in particular temporal variation we can find models (Mh, Mhww) that fit the Levantine intensity data more closely (**Table 8**), in particular visually appearing to describe the perceived spikes better in intensity (**Figure 4**). However, in the case of enforced westward drift (model Mhww compared to Mww) this comes at the expense of higher rms misfit to data from, on average, all other regions, in particular including all the rest of the world (**Tables 7**, **8**). In the unconstrained case (model Mh compared to M) the average misfit to the archeomagnetic and volcanic data is smaller, but the opposite is true for the sediments. These findings do not change notably if we consider both data and age uncertainties when evaluating the misfit.

Additionally, three ensembles of 9,000 bootstrap models were created based on the data underlying model C by varying the archeomagnetic and volcanic data by data and age uncertainties. Two of the ensembles use a westward drifting background model, with weaker and stronger regularization, respectively. The three ensembles thus have been labeled "no ww," "ww," and "ww2" for no, weakly and strongly enforced westward drift, respectively. More details about these models are given in the Supplemental Material.

Within the bootstrap ensemble models we find realizations of the data that can be fit closer by models with similar spatial and temporal complexity to the original data set (**Table 8** and Table S1, model C compared to no ww models and model Cww compared to ww models). Only in the case of more strongly enforced westward drift does model Mhww show a closer fit to the original Levantine data than found in the ensemble for the same region (model Mhww, **Table 8** compared to models ww2 in Table S1). Although our ensembles might be too small to cover the absolute best fit cases, this rather seems to be a consequence of lower temporal resolution in the ensemble models than model Mhww (see secular variation spectra in **Figure 6** compared to Figure S3). Our ensembles do not include any models that give close to minimum misfit for all three spike regions in the same realization (Figure S1). Again this does not necessarily reflect too small ensemble sizes, but might also be an indication that incompatible data exist in the data set even when data and age uncertainties are considered.


*RMS is the weighted root mean square misfit to all normally weighted input data, a,v for archeomagnetic and volcanic data, sed for sediment data.*

FIGURE 4 | Comparison of VADM from several models to regional data. The solid lines show predictions from models (A) C, green; A, blue; AS, brown; M, red; and (B) Cww, green; Mh, blue, Mww, brown; Mhww, red; see Table 6. Black dots with data and age error bars are recently published spike data sets, gray dots with data error bars only are from six regions with ∼1,500 km radius and lakes Escondido (Argentina, esc) and Mavora (New Zealand, mav). Colored dots mark model predictions for the actual data locations, while the solid line is the continuous prediction for the center of a region.

Korte and Constable Archeomagnetic Intensity Spikes

TABLE 8 | Weighted root mean square (rms) misfit of the global models to the three modified regional archeomagnetic and volcanic data sets and to the rest of these data types.


Models from the ensemble representing the best compromise fit to the spike data from the Levant, China and the Canary Islands were chosen as realizations with smallest misfit to both the Chinese and Canary Island data from the 50 models with minimum misfit to the Levantine data (Figure S1). These provide better fits to the Levantine and Chinese data than the comparable original models, except for the ww2 version (**Table 8** and model versions 4 in Table S2).

Considering the spectral power distribution the standard models show a lack of main field resolution starting from SH degree 6 (**Figures 6A,C**). Models built on westward drifting background models show similar spatial complexity to presentday models. All models indicate that the dipole moment was clearly higher than over the past 115 years, and probably played a role for the observed features by varying more strongly than at present (**Figures 6B,D**). Interestingly, this applies more to the younger of the two Levantine spikes in several of the models (**Figure 5A**). The two higher resolution models display a much stronger secular variation in large scales than all other models (**Figure 6D**). Models with double the rate of axial dipole change observed since 1,840 show indications of the Levantine spikes (**Figure 5A**) without fitting the maximum values. This can be achieved with comparable axial dipole acceleration as seen in the historical field (see gufm variability in the **Figure 5**), although large values would have to persist on somewhat longer timescales. The minimum misfit bootstrap ensemble models give similar results (Figure S3). A very close (visual) fit of the Levantine spikes is only obtained with higher temporal resolution models, which display more than five times larger rates of axial dipole change and require nearly three times higher acceleration of the axial dipole than the maximum found in the historical field.

Time-longitude plots of the radial field component at the CMB for 32◦N are shown in **Figure 7**. None of the models suggests that westward drift of a strong flux patch connects the spike observations from Asia and the Levant. Only when fit to the global data is compromised for a closer agreement with the westward drifting background model (model Mhww) is this movement largely preserved. All other models suggest that a strong flux patch related to the Levantine spikes grows and weakens mostly in place, perhaps moving slightly west during a weakening phase around 600 BCE. Similarly, timelatitude plots for the Levant area (35◦E) do not indicate any clear latitudinal drifts of flux patches (**Figure 8**). Some of the models (C,A,AS,Cww) suggest that some slight northward drift and extension of the patch might have occurred. Again the best fit bootstrap ensemble models show similar results (Figures S4, S5).

#### 5. DISCUSSION

Our modeling results indicate that the dipole moment, particularly the axial dipole strength is more closely correlated with spike features than previously thought, by being both higher and more variable than predicted by previous models. Geomagnetic field models gufm (Jackson et al., 2000) and the International Geomagnetic Reference Field IGRF (Thébault et al., 2015) indicate that the power in dipole variation has changed notably over the past 115 years and was lower by about 75% in 1900 than 2015 (**Figure 6**). Having no good reason to believe that maximum possible field variation occurs today it seems quite conceivable that the dipole variability was higher by 53% at times in the past, which is enough to find global field models that account for more than 80% of the observed overall spike intensity (model M in **Figures 4**, **5**).

Note that the models and in particular their power spectra depend strongly on the applied regularizations in space and time, which were chosen such that main field power of degrees two and higher appears close to a white spectrum dropping off for higher degrees at the CMB and that low secular variation degrees (high degrees for the two higher resolution models) are of similar magnitude to the present-day field. We considered these assumptions reasonable to obtain physically realistic fields. Our models do not include the option that the geodynamo process was drastically different from what it is today. The shape of the secular variation spectrum is strongly determined by the regularization norm, which essentially damps short-term variations more strongly than longer term variations by some function of spherical harmonic degree. We used the Ohmic heating norm in space and second derivative of the radial field in time, and did not find any norm that provided a closer resemblance of the spectral shape to the present-day field among

a range from simple power laws to other field quantities like horizontal gradients. Using global basis functions and a global regularization it is not possible to obtain arbitrarily detailed structure in certain regions that are much better constrained by dense data with low uncertainties than others. The main field and secular variation spectra (**Figure 6**) show that all models that are not artificially constrained by westward drift lack small scale (∼SH degree 5 and higher) spatial and temporal resolution. Considering that in the present-day field SH degrees of 5 and higher contribute to the field intensity at the Earth's surface by up to 3.7 µT regionally this can explain part of the missing amplitude fit to the spikes without having to invoke drastic differences in geodynamo processes.

This conclusion is supported by the fact that models with strongly westward drifting patches do nearly fit the amplitudes and rates of change of the spikes. The spatially and temporally localized small-scale structure of these models might represent the field morphology producing the spikes. The fact that these models show much higher secular variation power combined with worse fit to other globally distributed data might again be a consequence of the global nature of basis functions and regularization. Nevertheless, the fact that a full westward drift of a strong flux patch from China all the way to the Levant is not preserved in a model with comparatively weak regularization toward a westward drifting background suggests that this field behavior is not compatible with some other data.

Results from an ensemble of 9,000 models with underlying data varied within their magnetic and age uncertainties mainly confirm these findings. The highest intensity values and perceived rapid changes are fit less well by the minimum misfit models out of the ensemble than by some individual realizations (Figure S6 and **Figure 4**). Interestingly, the minimum misfit ensemble models suggest a slightly different interpretation of the

Chinese high intensity data around 1,200 BCE. In contrast to the individual models that in general show one maximum around that time most of the ensemble realizations displayed in Figure S6 show two maxima between 1,600 and 1,000 BCE. Nine thousand realizations do not sample every possible combination of data points and the coeval data will tend to be distributed to both sides within the large age uncertainties (Figure S7), so this result should be considered with caution. For other regions the individual and ensemble models shown give broadly similar predictions (Figure S5 and **Figure 4**). The maximum in the south American data of the same age as the Chinese high intensity values is not reflected by any of the models. An inspection of how the data underlying the minimum misfit ensemble models were shifted (Figure S7) does not reveal any clear patterns, except that the best fit to the Chinese data is obtained when the highest intensity occurred at the somewhat younger age of about 1,100 BCE instead of 1,300 BCE.

#### 6. CONCLUSIONS

We have compared observations of so-called archeomagnetic intensity spikes reported from several regions around the world and investigated potential causes. Age differences between spike observations in Asia, the Canary Islands and the Levant are compatible with a westward drifting intense flux patch as a common cause, but global field models constrained specifically

to test this hypothesis indicate that this explanation is not compatible with other field observations. On the contrary, these and additional test models suggest that the observed spikes are linked to individual flux patches, more or less growing and decaying in situ with only slight westward drift.

The spikes are clearly not global phenomena, but higher dipole moment than shown by previous models in particular around the times of the two Levantine spikes (980 and 740 BCE) facilitates explaining the spikes as normal non-dipole, regional intensity variation on top of a higher background field. Our test models suggest that this can be reconciled with other available data. However, this would require rather high dipole variability, with rates of change greater by more than 50% than directly observed in the last 180 years. This might still be within the normal range of geodynamo processes.

The recently published archeomagnetic and volcanic data reporting on spikes have been obtained following extremely strict quality criteria and, in particular in the case of the Levantine data could be confined to very low uncertainties in both data and age. The full available global data compilation contains other high intensity values at different ages and locations that might reflect similar geomagnetic field features, which are not as well constrained. Future re-assessments of previously published data in combination with new high quality data will indicate how rare or common archeointensity spike occurrences were in the past. It will then be interesting to test whether high dipole strength and variability are as relevant in other cases as they seem to be

#### REFERENCES


for the Levantine spikes. However, as long as the worldwide data distribution remains strongly inhomogeneous a clear description of such features in geomagnetic field models can only be obtained by models that do not rely on global basis functions in the way available spherical harmonic models do.

#### AUTHOR CONTRIBUTIONS

MK intially designed the study and carried out the analyses, adapting it iteratively following discussions of obtained preliminary results with CC. MK drafted the first version of the manuscript which was subsequently improved and agreed upon by contributions from both authors.

#### FUNDING

The collaboration was greatly facilitated by an extended research visit of MK to UCSD funded by the Alexander von Humboldt-Foundation. CC acknowledges NSF funding under grant number EAR 1623786.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2018.00017/full#supplementary-material

extreme behaviors of the geomagnetic field. Proc. Nat. Acad. Sci. U.S.A. 114, 39–44. doi: 10.1073/pnas.1616976114


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Korte and Constable. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Paleointensity Study on the Holocene Surface Lavas on the Island of Hawaii Using the Tsunakawa–Shaw Method

Yuhji Yamamoto\* and Ryo Yamaoka

*Center for Advanced Marine Core Research, Kochi University, Kochi, Japan*

Investigating volcanic paleointensity during the Holocene is important for linking archeointensity and sedimentary paleointensity. Across the globe, the island of Hawaii is one of the most studied subaerial locations. Many published data from Hawaii are accessible in the paleointensity databases, but it is necessary to reassess these data because they were determined with experimental protocols not incorporating a modern test for multi-domain particles and with relatively loose selection criteria. To obtain a new paleointensity dataset based on coercivity spectra rather than blocking temperature spectra, we applied the Tsunakawa–Shaw method to Holocene surface lavas collected from 34 sites on the island of Hawaii. In total, 135 successful results were obtained after applying the specimen-level selection criteria, yielding 22 site-mean Tsunakawa–Shaw paleointensities (TS dataset) that fulfilled the site-level selection criteria. We compared the TS dataset with the IZZI dataset, which includes recently reported blocking-temperature-based paleointensities determined by the IZZI Thellier method with the stringent criteria CCRIT. The cumulative distribution function (CDF) curve of the TS dataset, except for three sites, almost overlaps that of the IZZI dataset, and the medians coincide with each other, 42.9 µT (*N* = 19) for the TS dataset and 43.5 µT (*N* = 28) for the IZZI dataset. The coincidence of the CDF curves suggests equivalent reliability of the Tsunakawa–Shaw method and the IZZI Thellier method. The Holocene paleointensity variation in Hawaii is thought to be reliably characterized by both the TS dataset and the IZZI dataset: overall, the paleointensity throughout the Holocene is suggested to be higher than the present-day field. It is also inferred that there are possible decadal- and centennial-timescale large-intensity variations between ∼1,800 and ∼2,000 cal BP, and between ∼3,000 and ∼3,500 cal BP. The latter variation might be inferred as a global-scale dipolar phenomenon, as consistent paleointensity results are reported from archaeological sources in the Levant by the IZZI Thellier method.

Keywords: paleointensity, Tsunakawa-Shaw method, coercivity, Hawaii, holocene

#### INTRODUCTION

Investigating volcanic paleointensity during the Holocene, namely the last 10 kyr, is important for linking archeointensity (paleointensity from archeological materials) and sedimentary paleointensity because we can utilize multiple materials to establish a reliable master curve of the paleointensity variation for that period. As pointed out by Cromwell et al. (2017), the island

Edited by:

*Ron Shaar, Hebrew University of Jerusalem, Israel*

#### Reviewed by:

*Anita Di Chiara, Lancaster University, United Kingdom Nadia Solovieva, University College London, United Kingdom*

> \*Correspondence: *Yuhji Yamamoto y.yamamoto@kochi-u.ac.jp*

#### Specialty section:

*This article was submitted to Geomagnetism and Paleomagnetism, a section of the journal Frontiers in Earth Science*

> Received: *31 January 2018* Accepted: *16 April 2018* Published: *01 May 2018*

#### Citation:

*Yamamoto Y and Yamaoka R (2018) Paleointensity Study on the Holocene Surface Lavas on the Island of Hawaii Using the Tsunakawa–Shaw Method. Front. Earth Sci. 6:48. doi: 10.3389/feart.2018.00048* of Hawaii is one of the most studied subaerial locations in the world for Holocene volcanic paleointensity. Many published data from the island are accessible in paleointensity databases, such as GEOMAGIA (Korhonen et al., 2008; Brown et al., 2015) and MagIC (Tauxe et al., 2016). For example, if we apply a sitelevel selection criteria with (1) a minimum of three successful results for a site (N ≥ 3) and (2) successful results providing a site mean with a standard deviation < 15% of the mean (σ ≤ 15%), 55 site-mean volcanic paleointensities obtained by the Thellier-type method with pTRM checks can be selected from the GEOMAGIA50.v3 database (Brown et al., 2015): 48 site-means from surface lavas (GEOM-SL dataset) and another 7 site-means from Hawaiian Scientific Drilling Project (HSDP) cores (GEOM-HSDP dataset). The GEOM-SL dataset consists of 3 data (Coe et al., 1978), 5 data (Tanaka and Kono, 1991), 11 data (Mankinen and Champion, 1993), 1 data (Cottrell and Tarduno, 1999), 3 data (Chauvin et al., 2005), 18 data (Pressling et al., 2006), and 7 data (Pressling et al., 2007). The GEOM- HSDP dataset is comprised of 3 data (Teanby et al., 2002) and 4 data (Laj et al., 2002).

It is appropriate to reassess these data because they were determined with experimental protocols not incorporating a modern test for multi-domain (MD) particles and with relatively loose selection criteria. Since the study by Levi (1977), MD particles have been known to yield inaccurate paleointensity estimates, and their relative content in a specimen can be characterized by bulk hysteresis data. Paterson et al. (2017) found an evident stability trend in hysteresis data for both geological and archeological materials, named bulk domain stability (BDS). They concluded that specimens having lower BDS values, namely higher contents of MD particles, result in more curved Arai diagrams yielding more inaccurate results. Regarding selection criteria, one of a set of the latest stringent criteria is that adopted by Cromwell et al. (2015), later named "CCRIT" (Tauxe et al., 2016). Cromwell et al. (2015) demonstrated its robustness by reanalyzing the original measurement-level data obtained from the Kilauea 1960 lava by Yamamoto et al. (2003) with the CCRIT and observed an improvement of the site mean from 49.0 ± 9.6 µT (N = 17) to 39.1 ± 5.0 µT (N = 8), which is reasonably consistent with the expected field intensity of 36.2 µT.

Cromwell et al. (2017) recently obtained 22 new site-mean paleointensities by the IZZI Thellier method with the CCRIT from newly collected glassy volcanic materials from Holocene surface lavas on the island of Hawaii. They compared these data to the previously published data based on cumulative distribution function (CDF) diagrams together with six sitemean paleointensities of the same quality that were obtained from historical surface lavas on the island of Hawaii by Cromwell et al. (2015). They found that the CDF curve of their dataset (IZZI dataset; median of 43.5 µT, N = 28) was shifted to lower paleointensity values than those of the published Thellier-type data (median of 54.5 µT, N = 74). Unfortunately, they could not conclude a concrete cause for the shift because of inaccessibility of the original measurementlevel data, but Holocene paleointensity in Hawaii is thought to be more reliably characterized by the IZZI dataset because it was based on the modern technique and the stringent criteria. In order to assess reliably of the IZZI dataset, it is necessary to obtain an independent paleointensity dataset with modern techniques and/or selection criteria. In this study, we applied the Tsunakawa–Shaw method (Tsunakawa and Shaw, 1994; Yamamoto et al., 2003) to surface lavas collected from the island of Hawaii, most of which are Holocene ones, to obtain a new paleointensity dataset based on coercivity spectra rather than on blocking temperature spectra.

#### MATERIALS AND METHODS

#### Samples

The surface lavas used in this study were collected as standard one-inch paleomagnetic cores from 34 sites (HA31, 33–63, 65, and 66) on the island of Hawaii (**Figure 1** and **Table 1**). Six of these sites are historical lava flows (HA31, 1840 CE; HA58, 1881 CE; HA59, 1855 CE; HA60, 1935 CE; HA61, 1960 CE; HA62, 1955 CE), whereas the other 28 sites are older lava flows with <sup>14</sup>C ages reported in Rubin et al. (1987) (the radiocarbon laboratory ID number by the United States Geological Survey (USGS) is indicated for each site in **Table 1**). Conventional radiocarbon ages of these sites are between 230 and 9,540 year BP, except for two sites (HA55, 14,015 year BP; HA63, 23,840 year BP). Tanaka and Kono (1991) and Tanaka et al. (1995) previously reported Coe–Thellier paleointensities from the seven sites of HA31, 33, 35, 36, 48, 56, and 60, the site means of which range between 35.0 and 61.8 µT.

#### Methods

#### Rock Magnetic Experiments

A chip sample of ∼5 mg was used for acquisition of a thermomagnetic curve for each site. Using a magnetic balance (NMB-89, Natsuhara Giken, Osaka, Japan), the sample was heated from room temperature to 700◦C and then cooled to 50◦C, with constant application of a field of 500 mT in vacuum (1–10 Pa). Similar chip samples were subjected to hysteresis measurements for each site to determine the hysteresis parameters of saturation magnetization (Ms), saturation remanent magnetization (Mrs), coercive force (Bc), and coercivity of remanence (Brc). The measurements were conducted by vibrating sample magnetometer (MicroMag 3,900 VSM, Lake Shore Cryotronics Inc., USA).

#### Tsunakawa–Shaw Paleointensity Experiment

We applied the Tsunakawa–Shaw method (Tsunakawa and Shaw, 1994; Yamamoto et al., 2003) to 169 specimens cut from one-inch cores of all sites to determine absolute paleointensities based on coercivity spectra (this method mainly utilizes alternating field (AF) demagnetizations). For most of the specimens, we used an automated spinner magnetometer with an AF demagnetizer (DSPIN, Natsuhara Giken) to measure remanent magnetizations and to conduct stepwise AF demagnetizations and impart ARMs with the peak AF of 180 mT. A spinner magnetometer (ASPIN-A, Natsuhara Giken) and an AF demagnetizer (DEM-8601C, Natsuhara Giken) with the peak AF of 120 mT were used for some specimens. ARMs were imparted by a DC bias field of 50 µT with the peak AFs, setting the bias field direction approximately parallel to the NRM and laboratory TRM directions. To impart

laboratory TRMs, the specimens were heated to 610◦C in a vacuum of 1–10 Pa, maintained at that temperature for 15 (TRM1) and 30 (TRM2) minutes, and then cooled to room temperature for ∼3 h using a thermal demagnetizer with a built-in DC field coil (TDS-1, Natsuhara Giken). Throughout this process, a DC field of 20–60 µT was applied. Prior to stepwise AF demagnetization of each remanence, lowtemperature demagnetization (LTD; Ozima et al., 1964) was conducted on a specimen by soaking it in liquid nitrogen for 10 min and then leaving it at room temperature for 30 min in a zero field. The detailed procedures are described in Yamamoto and Tsunakawa (2005).

Based on the results obtained, NRM–TRM1<sup>∗</sup> and TRM1– TRM2<sup>∗</sup> diagrams were constructed for each specimen, where TRM1<sup>∗</sup> and TRM2<sup>∗</sup> denote the TRM imparted in the first (TRM1) and second (TRM2) heating, as corrected using the technique of Rolph and Shaw (1985). Paleointensity was estimated from the linear segment of the NRM–TRM1<sup>∗</sup> diagram when the ARM correction was judged to be valid based on the linear segment of the TRM1–TRM2<sup>∗</sup> diagram. We adopted the following selection criteria, which are similar to those used with the Tsunakawa–Shaw method in recent paleointensity studies (e.g., Yamamoto et al., 2010, 2015; Yamazaki and Yamamoto, 2014):



*Site, site number; Lat., Long., latitude and longitude of each sampled site;*

*<sup>14</sup>C ID, USGS Radiocarbon Laboratory ID number referred from Rubin et al. (1987); Age, conventional radiocarbon age with its standard deviation by Rubin et al. (1987);*

*Intensity, mean paleointensity with its standard deviation;*

*VADM, virtual axial dipole moment if the site-mean fulfills the site-level selection criteria; Nint, number of specimens used to calculate the site-mean paleointensity;*

*slopeA*1*, mean slopes in the ARM0–ARM1 diagrams for each site*

*Js-T, classification of the thermomagnetic curve.*


segment is unity within experimental errors (1.05 ≥ slope<sup>T</sup> ≥ 0.95) as proof of the validity of the ARM correction.

#### RESULTS

#### Rock Magnetic Experiment

The thermomagnetic curves were of the three types of A, B, and C. Type A was recognized for the curves from three sites (HA34, 44, and 56), which are characterized by a single Curie temperature (Tc) phase at 490–560◦C (**Figure 2A**). The type B curves show two T<sup>c</sup> phases at 100–250◦C and 500–550◦<sup>C</sup> (**Figure 2B**), and these were seen in the curves from 21 sites (HA31, 35, 37, 38–42, 45, 46, 48, 51–55, 57–59, 63, and 66). Type C is categorized as an intermediate type between type A and type B (**Figure 2C**), and it was obtained from 10 sites (HA33, 36, 43, 47, 49, 50, 60–62, and 65). The main magnetic carriers for the three types are considered to be titanomagnetites of different Ti contents reflecting different degrees of deuteric oxidation.

Ratios of the hysteresis parameters, Brc/B<sup>c</sup> and Mrs/M<sup>s</sup> , (Table S1) generally lie on the bulk domain stability (BDS) trend (Paterson et al., 2017), as shown by **Figure 3**. Paterson et al. (2017) introduced the BDS value, which can be calculated easily from the ratios, and showed based on their compilations of historical and laboratory-controlled paleointensity and hysteresis parameter dataset that specimens with BDS values lower than 0.10 are less likely to yield a meaningful paleointensity. BDS values of the present specimens range between 0.121 and 0.830 (average = 0.578, standard deviation = 0.146; Table S1), so they fulfill the minimum reliability criterion.

#### Tsunakawa–Shaw Paleointensity

After applying the selection criteria, 135 successful results were obtained from a total of 169 specimens of the 34 sites (success rate of 80%). The resultant paleointensities range between 13.6 and 81.6 µT (Table S1), and 91% of these were yielded from more than 50% of the total extrapolated NRM fractions (f<sup>N</sup> ≥ 0.50). The other 27 results were rejected mainly due to non-unity slopes of the linear segments in the TRM1–TRM2<sup>∗</sup> diagrams. Representative successful and rejected results are illustrated in **Figure 4**. Among the results, those from 22 sites fulfilled the sitelevel selection criteria (site-level success rate of 65%): (1) the successful results were yielded from more than three individual specimens for each site (N ≥ 3) and (2) they resulted in a site-mean paleointensity with a standard deviation < 15% of the site mean (σ ≤ 15%). These 22 site-mean Tsunakawa–Shaw paleointensities (TS dataset) were derived with a relatively small amount of ARM corrections (≤25% suggested from slopeA1 values in **Table 1**). They are not significantly smaller than the present-day field around Hawaii of 34.6 µT (IGRF2015; Thébault et al., 2015), except for the three sites of HA35 (1,400 year BP), HA55 (14,015 year BP), and HA66 (9,540 year BP) (**Figure 5** and **Table 1**).

#### DISCUSSION

#### Direct Comparison of the Site-Mean Paleointensities Between the TS Dataset and the Previously Published Coe–Thellier Data

Some site means of the TS dataset can be compared directly with the site-mean Coe–Thellier paleointensities of the same quality (N ≥ 3 and σ ≤ 15%) published in previous studies, namely the data from the GEOM-SL dataset (section Introduction). Tanaka and Kono (1991) reported Coe–Thellier paleointensities from sites HA31, 33, 48, and 56 using sister specimens of the present study. Mankinen and Champion (1993) published one site mean from the lava with the USGS radiocarbon laboratory ID number of W3881, which is considered to be the same lava as that at site HA56. Pressling et al. (2006) also reported two site means from the lavas with the USGS radiocarbon laboratory ID numbers of W4047 and W4006, which are thought to be the same lavas as at sites HA39 and 52. **Figure 6** illustrates the comparison between the Tsunakawa–Shaw and the Coe–Thellier paleointensities for these site means. At the one standard deviation level, the differences between the two site-mean paleointensities are within ± 10% for sites HA31, 33, 39, and 48. On the other hand, the site-mean Coe–Thellier paleointensities are 32–48% higher than the site-mean Tsunakawa–Shaw paleointensities for sites HA52 and 56.

Yamamoto et al. (2003) observed a similarly higher value from the Kilauea 1960 lava (expected intensity of 36.2 µT): their Coe– Thellier experiment resulted in a site-mean paleointensity of 49.0 ± 9.6 µT (N = 17), which is 24% higher than that of 39.4 ± 7.9 µT (N = 8) obtained by their Tsunakawa–Shaw experiment. They suggested that a possible cause of the higher result in the Coe–Thellier experiment is contamination of thermochemical remanent magnetization (TCRM) into NRM, and Yamamoto (2006) later investigated the possibility and reported supporting results. Another possible cause is MD-like remanence because the Coe–Thellier site-mean paleointensity was improved to 39.1 ± 5.0 µT (N = 8), which almost coincides with the Tsunakawa– Shaw site-mean paleointensity, by reanalysis applying the CCRIT (Cromwell et al., 2015). These results indicate that the Coe– Thellier method with ordinary selection criteria can sometimes result in high paleointensities compared with those of the Tsunakawa–Shaw method and the Coe–Thellier method with CCRIT. Interestingly, it is reported in Cromwell et al. (2017) that the published Thellier-type data resulted in a higher median of 54.5 µT (N = 74) than the 43.5 µT (N = 28) median value from the IZZI Thellier data with CCRIT for Holocene Hawaiian lavas.

#### Holocene Paleointensity Variation in Hawaii

As reviewed in section Introduction, Cromwell et al. (2017) found that the CDF curve of the IZZI dataset (median of 43.5 µT, N = 28) was shifted to lower paleointensity values than those of the published Thellier-type data (median of 54.5 µT, N = 74) for the Holocene. Because the IZZI dataset was obtained by the modern technique and the stringent criteria, paleointensity in Hawaii during Holocene is thought to be more reliably characterized by the IZZI dataset. It is suggested that the paleointensity in Hawaii during Holocene is higher than the present-day field around Hawaii of 34.6 µT (IGRF2015; Thébault et al., 2015) but is not high enough to result in the median of 54.5 µT. This thought is based on the comparison made on the blocking-temperature-based paleointensity data in Cromwell et al. (2017). We have obtained a new independent coercivity-based paleointensity dataset, the TS dataset, for comparison to assess further the thought.

The IZZI dataset covers the time period between present and 6,500 cal BP. After excluding the three site means from HA40 (8,550 year BP), HA55 (14,015 year BP), and HA66 (9,540 year BP) to have the same time coverage,

we compared the TS dataset (present to 4,340 year BP) with the IZZI dataset on a CDF diagram (**Figure 7**). The CDF curve of the TS dataset almost overlaps the curve of the IZZI dataset and the medians coincide: 42.9 µT (N = 19) for the TS dataset and 43.5 µT (N = 28) for the IZZI dataset. The Tsunakawa–Shaw method estimates paleointensity based on coercivity spectra, whereas the IZZI Thellier method uses blocking-temperature spectra to yield paleointensity. The two datasets are not necessarily from the same lava flows (**Figure 1**). However, considering these independences, the coincidence of the CDF curves suggests equivalent reliability of the Tsunakawa–Shaw method and the IZZI Thellier method, which was already demonstrated for smaller datasets by Yamamoto et al. (2010) and Ahn et al. (2016). Cromwell et al. (2017) argued that total TRM-based paleointensity methods, including the Shaw method, are likely to underestimate the expected geomagnetic field strength in general because they do not test for SD-like behavior. The Tsunakawa– Shaw method, on the other hand, places much emphasis on the selective removal of MD-like components by lowtemperature and AF demagnetizations and also on correction to compensate possible magnetostatic interactions among the magnetic assemblages (e.g., Yamamoto and Tsunakawa, 2005; Yamamoto and Hoshi, 2008) and anisotropy of remanences (e.g., Yamamoto et al., 2015) using ARM. These improvements have not been incorporated into most of the Thellier-type experiments.

The Holocene paleointensity variation in Hawaii is thought to be reliably characterized by both the TS dataset and the IZZI dataset: they show generally lower paleointensities than the GEOM-SL dataset and do not show a good match with the global model of CALS10k.2 (Constable et al., 2016; **Figure 8**). The lower paleointensities are obvious if each data is binned with a 1,000 year interval (**Figure 9**). In the figure, the ages are indicated as calendar ages. For the TS dataset, the ages were recalculated using the CALIB 7.1 program (Stuiver et al., 2017) with the IntCal13 radiocarbon calibration curve (Reimer et al., 2013); for the IZZI dataset, the ages, directly referred from Cromwell et al. (2015, 2017), are based on the IntCal04 radiocarbon calibration curve (Reimer et al., 2004); for the GEOM-SL dataset, the ages were directly referred from the GEOMAGIA50.v3 database (Brown et al., 2015).

Although both the TS dataset and the IZZI dataset resulted in lower paleointensities than the GEOM-SL dataset, it is suggested that the paleointensity throughout the Holocene is overall higher than the present-day field (**Figure 8**). One of the prominent features suggested by Cromwell et al. (2017) is the possible decadal- and centennial-timescale large-intensity variation between ∼1,800 cal BP (∼ 67 µT) and ∼2,000 cal BP (∼35 µT). We recognize the paleointensities of ∼68 µT (HA33) and ∼50 µT (HA38) at around ∼1,900 cal BP in the TS dataset, which are supportive of this feature. Similar possible variation might also be inferred between ∼3,000 and ∼3,500 cal BP by the IZZI dataset and the newly obtained Tsunakawa-Shaw

**75**

FIGURE 5 | Resultant Tsunakawa–Shaw paleointensities according to site. Individual specimen-level paleointensities are indicated by red circles for all sites. Site-mean paleointensities are also shown by black squares with their one standard deviation (1σ) if they were yielded from more than three individual specimens (*N* ≥ 3) and the associated standard deviation was < 15% of the site mean (σ ≤ 15%).

10% high (red), equivalent (black), and 10% low (blue).

paleointensity data of ∼58 µT (HA50) at around ∼3,100 cal BP. The corresponding virtual axial dipole moment (VADM) is calculated to be 130 ZAm<sup>2</sup> , which is consistent with the VADMs between ∼2,900 and ∼3,000 cal BP reported from archaeological sources in the Levant by the IZZI Thellier method (Figure 5 in Shaar et al., 2016). Levant is almost opposite to Hawaii on the

globe and this variation might be inferred as a global-scale dipolar phenomenon. Though the two variations do not appear to match with the global model of CALS10k.2 (Constable et al., 2016; **Figure 8**), the model is compiled with a dominance of sediment data (Constable et al., 2016) and is not thought to reflect possible rapid variations in paleointensity.

#### CONCLUSIONS

2016) is also indicated for reference.

We applied the Tsunakawa–Shaw method to surface lavas collected from 34 sites, most of which are Holocene ones, on the island of Hawaii. The results of the rock magnetic experiments indicated that the main magnetic carriers are titanomagnetites of various Ti contents due to varying degrees of deuteric oxidation and that the specimens fulfilled the minimum reliability criterion for the BDS value. After applying the specimen-level selection criteria, 135 successful results were obtained from a total of 169 specimens, yielding 22 site-mean Tsunakawa–Shaw paleointensities (TS dataset) that fulfilled the site-level selection criteria (N ≥ 3 and σ ≤ 15%). They are not significantly smaller than the present-day field around Hawaii of 34.6 µT, except for three sites. Six site means of the TS dataset can be compared directly with site-mean Coe–Thellier paleointensities of the same quality published in previous studies. At the one standard deviation level, the differences between the two sitemean paleointensities are within ± 10% for four sites. For the other two sites, the site-mean Coe–Thellier paleointensities are 32–48% higher than those of the TS dataset. Previous studies indicate that the Coe–Thellier method with ordinary selection criteria can sometimes result in higher paleointensities than the Tsunakawa–Shaw method and the Coe–Thellier method with CCRIT.

For the same time period, the CDF curve of the TS dataset almost overlaps that of the IZZI dataset and the medians coincide [42.9 µT (N = 19) for the TS dataset; 43.5 µT (N = 28) for

#### REFERENCES

Ahn, H. S., Kidane, T., Yamamoto, Y., and Otofuji, Y. (2016). Low geomagnetic field intensity in the Matuyama Chron: palaeomagnetic study of a lava the IZZI dataset]. The coincidence of the CDF curves suggests equivalent reliability of the Tsunakawa–Shaw method and the IZZI Thellier method. The Holocene paleointensity variation in Hawaii is thought to be reliably characterized by both the TS dataset and the IZZI dataset: they show generally lower paleointensities than the GEOM-SL dataset and do not show a good match with the global model of CALS10k.2 (Constable et al., 2016). Overall, the paleointensity throughout the Holocene is suggested to be higher than the present-day field. We recognize the paleointensities of ∼68 and ∼50 µT at ∼1,900 cal BP in the TS dataset, which are supportive of the possible decadal- and centennial-timescale large-intensity variation between ∼1,800 cal BP (∼67 µT) and ∼2,000 cal BP (∼35 µT) observed by Cromwell et al. (2017). Similar possible variation might also be inferred between ∼3,000 and ∼3,500 cal BP by the IZZI dataset and the newly obtained Tsunakawa-Shaw paleointensity data of ∼58 µT (HA50) at around ∼3,100 cal BP. The corresponding VADM of 130 ZAm<sup>2</sup> is consistent with the VADMs between ∼2,900 and ∼3,000 cal BP reported from archaeological sources in the Levant by the IZZI Thellier method (Shaar et al., 2016), and this variation might be inferred as a global-scale dipolar phenomenon. Though the two variations do not appear to match with the global model of CALS10k.2 (Constable et al., 2016; **Figure 8**), it is not thought to reflect possible rapid variations in paleointensity.

#### AUTHOR CONTRIBUTIONS

YY designed research; YY and RY conducted the measurements and data analyses; YY wrote the paper.

#### FUNDING

This study was partly supported by the Kochi University Research Project Research Center for Global Environmental Change by Earth Drilling Sciences and JSPS KAKENHI Grant Number 15H05832.

#### ACKNOWLEDGMENTS

We thank Masaru Kono and Hidefumi Tanaka for providing the samples, and Yukako Nabeshima, Shinsuke Yagyu, and Nana Yoshikane for assistance with the measurements. Constructive comments by the associate editor Ron Shaar and the two reviewers improved the manuscript. We would like to thank Editage (www.editage.jp) for English language editing.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2018.00048/full#supplementary-material

sequence from Afar depression, East Africa. Geophys. J. Int. 204, 127–146. doi: 10.1093/gji/ggv303

Brown, M. C., Donadini, F., Korte, M., Nilsson, A., Korhonen, K., Lodge, A., et al. Constable (2015). GEOMAGIA50.v3: 1. General structure and modifications to the archeological and volcanic database. Earth Planets Space 67:83. doi: 10.1186/s40623-015-0232-0


archeomagnetic data from Tel Megiddo and Tel Hazor, Israel. Earth Planet. Sci. Lett. 442, 173–185. doi: 10.1016/j.epsl.2016.02.038


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Yamamoto and Yamaoka. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## A New View of Long-Term Geomagnetic Field Secular Variation

Steve P. Lund\*

*Department of Earth Sciences, University of Southern California, Los Angeles, CA, United States*

This study carries out a statistical analysis of high-resolution PSV records for the last ∼70 ka from three different regions of the Earth. We consider directional and intensity variability in each region on time scales of 10<sup>3</sup> -10<sup>5</sup> years in order to evaluate long-term PSV. We then compare those results with more traditional long-term PSV statistical studies averaged over ∼10<sup>6</sup> years. Three replicate PSV records from one region (subtropical North Atlantic Ocean) were averaged at overlapping 3 and 9 ka intervals. Variability in both scalar inclination and declination variability and vector angular dispersion are significant and coherent among the three records. The vector dispersion is relatively low for most of the time but contains two relatively narrow intervals (∼30–42 and 60–65 ka) of high dispersion. (Vector dispersion in all records was calculated after removing directions with true excursional VGPs, VGPs < 45◦ N). We have carried out a comparable statistical analysis on two other PSV records from other parts of the Earth (Chile margin; Philippines/Indonesia). The results for these three regions are comparable in their overall style of variability. The scalar directional variability from the Philippines/Indonesia is quite different in detail from the other two regions, as might be expected, but the scalar directional variability between the Western Hemisphere regions is remarkably consistent considering their distance from one another. This may be associated with them being on the same longitude swath and having some coherent dynamo activity occurring along that path. Three magnetic field excursions occur in the study interval. All three excursions are associated with the two highest vector dispersion intervals. Paleointensity records from the three regions were subjected to the same statistical analysis as the directions. These records are all coherent in their pattern of variability. The similarity in paleointensity variability on a global scale is expected even though the detailed scalar directional variability is not coherent on a global scale. The pattern of intensity variability is strongly correlated with the pattern of vector dispersion and excursions on a global scale—high (low) intensity is associated with low (high plus excursions) vector dispersion. The fact that regional directional variability is always larger than "normal" during low intensity/excursional intervals, even though the effect of true excursional directions was removed, suggests that we need to reevaluate what field variability was like during low intensity/excursional intervals on a global scale and how/why it was different from today's field (last 10<sup>4</sup> years).

Keywords: secular variation, excursions, paleointensity, dynamo, paleomagnetism

#### Edited by:

*Christopher Davies, University of Leeds, United Kingdom*

#### Reviewed by:

*Geoffrey Cromwell, United States Geological Survey, United States Peter Aaron Selkin, University of Washington Tacoma, United States*

> \*Correspondence: *Steve P. Lund slund@usc.edu*

#### Specialty section:

*This article was submitted to Geomagnetism and Paleomagnetism, a section of the journal Frontiers in Earth Science*

> Received: *14 December 2017* Accepted: *04 April 2018* Published: *04 May 2018*

#### Citation:

*Lund SP (2018) A New View of Long-Term Geomagnetic Field Secular Variation. Front. Earth Sci. 6:40. doi: 10.3389/feart.2018.00040*

#### INTRODUCTION

Magnetic compass readings have been associated with the Earth as their source since perhaps the thirteenth Century AD; experiments by, first, Petrus Peregrinus (1269) and then Alexander Needham (1600) put that theory on firm ground (see historical summary in Merrill et al., 1998). The fact that geomagnetic field measurements at the Earth's surface have changed over time (secular variation) was first noted by the Chinese by the thirteenth Century AD and later by Henry Gellibrand in Europe in 1634 (Merrill et al., 1998). We now have a good sense of historical secular variation (HSV) of the geomagnetic field for the last 400 years or so (e.g., Bullard et al., 1950; Yukutake, 1967; Yukutake and Tachinaka, 1968; Bloxham and Gubbins, 1985; Jackson et al., 2000). We also understand, to some degree, the source of that magnetic field through dynamo activity in the Earth's outer core (see Merrill et al., 1998 for overview).

Evidence for prehistoric secular variation, based on paleomagnetic measurements, was first noted by Chevaliier (1925) using lava flows from Mount Etna. Thellier and Thellier (1951, 1952) measured directional paleomagnetic secular variation (PSV) in archeological materials (baked hearths and pottery) from Europe. The earliest directional PSV records from sediments were derived from deep-sea sediment cores and exposed varved lake sediments (McNish and Johnson, 1938; Johnson et al., 1948). Thellier (1937a,b) was the first to determine paleointensity variations in archeological materials.

The magnetic field at any one location has a mean direction after averaging secular variation for some period of time—several thousand years or longer. The actual amount of time necessary to average the full range of local secular variation is still not well established (e.g., Merrill et al., 1998). However, paleomagnetic studies suggest that the average field direction at a locality during stable polarity (the period between global polarity reversals) is generally that of a dipole located at the Earth's center and aligned with the rotation axis; this hypothesis is termed the geocentric axial dipole (GAD) hypothesis (see Merrill et al., 1998 for overview).

Merrill et al. (1998) noted that "the interpretation of paleomagnetic results has always depended upon the fundamental hypothesis that the time-averaged geomagnetic field is that of a geocentric axial dipole." But, a large number of statistical PSV studies have made it clear that there are systematic (but usually subtle) differences in local field directions from the GAD when averaged over millions of years (e.g., Johnson and Constable, 1996; McElhinny et al., 1996; Merrill et al., 1998; Johnson et al., 2008) or shorter intervals (e.g., Lund et al., 1988, 2016a). The unresolved problem is what causes these non-GAD field directions and how do they vary over time-scales less than a million years.

This paper focuses on a new set of replicate, high-resolution PSV records recovered from lake and deep-sea sediment cores for intervals of 50–100,000 years that demonstrate notable long-term variations in the statistical properties of local field. These new statistical results are developed and then compared with traditional results to suggest that a more fundamental understanding of geomagnetic field variability and the underlying dynamo source activity can be discerned from the high-resolution PSV data.

#### TRADITIONAL ANALYSIS OF LONG-TERM PSV

The normal method of long-term PSV analysis is to statistically average a spatially regional group of lava-flow data, which are uncorrelated in time and cover a long time interval (>10<sup>5</sup> - 10<sup>6</sup> years). The expectation is that individual data points are randomly placed in their age distribution and no two data points are strongly serially correlated in both age and direction. These spot averages are then compared on a global scale using spherical harmonic analysis or latitude/longitude transects (e.g., Johnson and Constable, 1996; McElhinny et al., 1996; Merrill et al., 1998; Harrison, 2007; Johnson et al., 2008). In this manner 1I (the difference between local average inclination and that expected for the GAD), and angular dispersion of field vectors or VGPs (**Figures 1A,B**, solid dots) can be determined on a global scale (see summary in Merrill et al., 1998).

Another method of analysis is to consider long sequences of successive PSV data points from deep-sea cores (e.g., Opdyke, 1972; **Figure 2**) or dated lava flows (e.g., McDougall et al., 1977) and look for long-term variations in statistical patterns at one site. One shortcoming is a lack of detailed age control and, thus, a need to average over long time intervals (>10<sup>5</sup> years) to minimize age uncertainty. Another shortcoming is the lack of certainty about the PSV variability itself. In sediments, that is due to low sedimentation rates (such as the data in **Figure 2**) and sediment averaging that preclude detailed PSV. In the case of lava flow sequences, unknown data intervals between successive flows preclude detailed PSV. In both cases, statistical studies are carried out over long intervals (>10<sup>5</sup> years) to be sure that enough recorded PSV is averaged.

The key element in all of these studies is that there is not strong evidence for detailed directional field variability as is common in studies of historical secular variation or Holocene PSV (e.g., Creer et al., 1983; Lund, 1996). The resolution of such studies is typically a million years or longer. Also, there is only limited ability to assess paleointensity variability associated with the vector variability.

#### NEW LONG PSV TIME SERIES

There are now sufficient, published full-vector PSV time-series (inclination, declination, paleointensity) for the last glacial cycle (0–71 ka) from deep-sea and lake sediment sequences to begin to assess long-term PSV on a global scale using high-resolution, serially-correlated data. Appendix 1 summarizes the criteria we use to distinguish what constitutes high-resolution PSV records. **Figure 3** shows the locations of long PSV records from 6 different regions (A–F). All of these PSV records come from paleomagnetic studies of deep-sea or lake sediment sequences. There are replicate PSV records within each of these regions with

(A) and VGP (B) angular dispersion and selected more recent (∼10<sup>3</sup> years) results. The data come from a mixture of lava-flow and sediment studies. Long-term averages are indicated by solid dots. ∼3 and 10–20 ka averages (all latest Quaternary in age) are indicated by open circles and squares respectively. Smooth curves estimate the general pattern of variability for each average vs. latitude. VGP data from all studies flip Southern Hemisphere data to the Northern Hemisphere for averaging. Gray diamonds indicate statistical averages for our studied PSV data sets (Table 1) (modified from Merrill et al., 1998; Lund et al., 2016a).

good-quality age control, to permit correlation and analysis of directional-waveform and relative-paleointensity variability.

Region A (Bering Sea/Gulf of Alaska) contains multiple long PSV records from IODP Sites U1343-1345 from the Bering Sea (Takahashi et al., 2011; Lund et al., 2016b) and U1418-1419 from the Gulf of Alaska (Jaeger et al., 2014). Region B (SW USA) contains records from Mono Lake, California (Lund et al., 1988), Pyramid Lake, Nevada (Lund et al., 2017a), and Summer Lake, Oregon (Negrini et al., 2000; Zic et al., 2002). Region C (subtropical North Atlantic Ocean) contains records from the Bermuda Rise (89-9P, ODP Site 1062) (Keigwin et al., 1998; Schwartz et al., 1998; Lund S. et al., 2001; Lund S.P. et al., 2001) and the Blake-Bahama Outer Ridge (88-10P, JPC-14, ODP Sites 1060, 1061) (Keigwin et al., 1998; Schwartz et al., 1998; Lund S. et al., 2001; Lund S.P. et al., 2001). Region D (equatorial west Atlantic Ocean) contains records from the Demerara Rise (9GGC, 25GGC, CDH11, CDH14) (Huang et al., 2014; Lund et al., 2017b) and the Amazon Fan (CDH83, CDH86) (Nace et al., 2014). Region E (SE Pacific Ocean/Chile Margin) contains records from ODP Sites 1233 and 1234 (Mix et al., 2003; Lund et al., 2007). Region F (Philippines/Indonesia) contains records from MD97-2134 (Blanchet et al., 2006; Lund et al., 2018) and MD98-2181 (Stott et al., 2002; Lund et al., 2018).

This paper will consider, in detail, the statistical pattern of directional field variability from three Region C PSV records and compare it to individual PSV records from Regions E and F. Paleointensity variability will also be compared with all the records from these three regions. **Figures 4**, **5** show the inclination and declination variability for three PSV records from Region C (subtropical North Atlantic Ocean) over the last 70 ka. These three records span a distance of ∼1,500 km and all contain individual inclination and declination records that display good serial correlation. The individual inclination and declination features can be correlated among these three records over the entire time interval.

#### STATISTICAL ANALYSIS OF LONG-TERM DIRECTIONAL PSV TIME SERIES

Holocene PSV records from around the world suggest that PSV is largely a centennial-millennial scale process with longerterm trends. The traditional long-term PSV studies averaged uncorrelated directional data over ∼10<sup>6</sup> years time scales or longer. Well-dated high-resolution PSV time series permit statistical analysis over significantly shorter time intervals. The argument is that the dynamo process in the Earth's outer core that generates the geomagnetic field is operating on much shorter (∼10<sup>2</sup> -10<sup>4</sup> years) time scales and, if we want to properly understand the space/time pattern of dynamo activity there is merit in considering statistical averages at those shorter time scales.

There is no magical time scale over which to average the field to study dynamo activity, but previous PSV studies suggest that simple PSV directional features or waveforms (e.g., Lund, 2007) that define the pattern of PSV variability between successive extremes in inclination or declination (**Figures 4**, **5**) operate on time scales of ∼0.7–1.5 ka (e.g., summaries of regional PSV studies in Creer et al., 1983) and repeating patterns of more complex directional waveforms have been noted by Lund et al. (1988) and Lund (1996) with repeating time intervals of 2.5–3.5 ka. On that basis, we have started by averaging

780 ka) and an average inclination close to axial dipole expectation. This record is a low-resolution (∼2 cm/kyr sediment accumulation rate) that does not resolve high-resolution PSV features. But, it does describe the overall near GAD pattern of field during the current Brunhes normal chron and is representative of the available long-term sediment PSV records before the advent of more high-resolution sediment PSV studies such as we present here. Modified from Opdyke (1972).

the three Region C PSV records at 3 ka intervals with 1.5 ka overlap between successive averages. We have averaged scalar inclination, declination, vector direction, and angular dispersion of directions. The vector dispersions were calculated after removing directions that have excursional VGPs (VGP latitudes less than 45◦ N; see Merrill et al., 1998; Johnson et al., 2008 for more detailed discussion). The results of this analysis are plotted in **Figure 6** (open circles). The vector inclination variability is not significantly different from the scalar inclination variability plotted in **Figure 6**. We then carried out a second round of statistical averaging over a longer time average −9 ka with a 4.5 ka overlap between successive averages. Those results are also shown in **Figure 6** (closed circles).

Both the 3 and 9 ka averages show a strong serial correlation between successive values. This suggests that a shorter averaging interval is not needed to begin to see the temporal range of PSV

variability. The 9 ka averages always have maximum or minimum values synchronous with 3 ka averages. This suggests that there is significant dynamical variability in the range of ∼3–9 ka. Even so, there is also clear evidence for a long, >60 ka, trend in inclination that is indicated by a dashed line in **Figure 6**. Older PSV records from these sites (Keigwin et al., 1998) show that the average inclinations of the three records (∼45◦ ) are not significantly different from the long-term inclination averages at these sites for the entire Brunhes Epoch, so this trend does disappear in longer time intervals.

closed bars. Times of known global magnetic field excursions are indicated by vertical dashed lines.

The range of directional variability in 3 ka averages is dramatic. Individual 3 ka intervals may have inclination and declination averages as much as 30◦ different from other intervals. The 9 ka variability still sees individual intervals with inclination and declination averages as much as 20◦ different from other intervals. These values are consistent with the differences noted at mid Northerly latitudes around the world in historical maps of spatial inclination and declination variability (e.g., Merrill et al., 1998). Thus, the 3 and 9 ka statistical variability may simply reflect the overall spatial/temporal dynamics of dynamo activity. The detailed pattern of scalar inclination and declination variability is also consistent among the three records as indicated by inclination features A–I and declination features A–J noted in **Figure 6**. As noted in other studies (Lund, 1996, 2007), it makes sense that the dynamical variability seen in **Figure 6** should be consistent among these three PSV records that are located less than a few thousand km apart in Region C.

The vector statistical variability that results from this study is similar in pattern to the scalar variability and more directly comparable to long-term PSV studies noted above. The pattern of scalar inclination variability is not significantly different from the pattern of vector inclination variability. Thus, one could think of the varying inclination averages as varying 1I anomalies (site inclination minus site axial dipole expectation, **Figure 1**). It is clear that the 1I anomaly can vary by as much as 20◦ over ∼9 ka intervals. The directional angular dispersion can also vary by almost a factor of two, with most time spent in relatively low dispersion, ∼10◦–12◦ .

These values can be compared with the overall (∼10<sup>5</sup> years) statistical averages of the entire three time series. The results are summarized in **Table 1** and plotted in **Figure 1**. The ∼10<sup>5</sup> years length of these PSV records produce long-term statistical results that are generally comparable to the traditional long-term (∼10<sup>6</sup> ) PSV results. Our results suggest, however, that these long-term averages smear out the variability at 10<sup>3</sup> -10<sup>4</sup> years that probably reflects dynamo activity.

We can also compare the statistical results from Region C with similar results for Regions E and F (**Figure 3**). We have carried out a comparable statistical analysis on MD98-2181 from Region F and ODP Site 1233 from Region E. The results are shown in **Figure 7**. The PSV data for these records have been previously published and compared with other records from the same region (Region E: Mix et al., 2003; Lund et al., 2007) (Region F: Stott et al., 2002; Blanchet et al., 2006; Lund et al., 2017a) to verify the PSV pattern and chronology.

All three regions display an overall sense of directional variability comparable to that of Region C. It is interesting to note, however, that the detailed pattern of scalar inclination and declination variability is very similar between Regions C and


TABLE 1 | Long-term statistics of selected PSV time series.

Region E, as noted by inclination features A–I and declination features A–J (**Figure 7**). Both regions also note the same longterm trend in inclination (dashed trend lines in **Figure 7**). However, the variability in Region F is quite different in its detailed pattern from the other two regions. One possible explanation for this is that Regions C and E are along the same longitude swath (**Figure 3**) that aligns with the historical location of the North American high-altitude flux lobe (Bloxham and Gubbins, 1985). If that flux lobe is evidence of an outer-core wide pattern of dynamo activity (e.g., convection roll; Gubbins and Bloxham, 1987) along this longitude swath, then it may be that our scalar directional variability in Regions C and E are part of a coherent regional dynamo pattern. The variability in Region F, located along a very different longitude swath, could then be associated with a different longitudinally-controlled dynamo pattern.

We can also evaluate the vector dispersions (or, in principle, the VGP dispersions) associated with each PSV record. **Figures 6**, **7** show that there are two intervals of enhanced vector dispersion in the last 70 ka—∼30–42 and ∼60–65 ka (black horizontal bars in **Figures 6**, **7**). This pattern is the same in all three studied regions. The pattern does not go hand-in-hand with the scalar patterns of inclination and declination variability. Also, it appears that the general "level" of vector dispersion, estimated by horizontal lines in **Figures 6**, **7**, is quite a bit lower in each record than the two, relatively narrow, intervals of enhanced vector dispersion. Also, the general "level" of vector dispersion is significantly lower than the long-term average vector dispersion for each record noted in **Table 1**. It seems that this lower level of vector dispersion may be the more proper parameter for considering the "level" of dynamo activity in a region over time. The higher long-term averaged dispersion (**Table 1**) is then associated with the average of two different "states" of the local dynamo activity—one that is lower in dispersion and operates for most of the time and a second "state" that is higher in dispersion, but operates only for a more limited portion of the time.

#### RELATIONSHIP OF LONG-TERM DIRECTIONAL PSV TO EXCURSIONS

There are three excursions reported repeatedly during the last glacial cycle: the Mono Lake Excursion (34 ka; Liddicoat and Coe, 1979), the Laschamp Excursion (41 ka; Bonhommet and Babkine, 1967; Bonhommet and Zahringer, 1969) and an unnamed excursion in Stage 4 (61 ka; Nowaczyk and Baumann, 1992; Nowaczyk et al., 1994). The excursion times are shown as vertical dotted lines in **Figures 6**, **7**, **9**–**11**. These excursions are all associated with intervals of notably higher field vector dispersion (**Figures 6**, **7**) and low paleointensity (**Figures 9**–**11**).

The higher vector dispersion, in each case, occurs even though the actual excursional field directions in all-time series were deleted from the statistical analyses using a 45◦ latitude cutoff for VGPs (e.g., Johnson et al., 2008). This suggests that field variability surrounding excursions is routinely more highamplitude even when individual directions are not excursional. The evidence is not certain, however, that higher amplitude dispersion really precedes excursion intervals, an observation previously noted for the Mono Lake Excursion (Liddicoat and Coe, 1979; Lund et al., 1988). Regions C, E, and F share the relationship of higher amplitude dispersion associated with excursions even though the overall pattern of directional variability in Region F is significantly different from the other two regions. One possibility is that the low paleointensity associated with all the excursional intervals is the driving force in causing higher amplitude dispersions and as well as localized excursions in selected regions.

There is no indication that directional statistical parameters in Regions C, E, or F behave in an anomalous manner before excursional intervals. This might indicate that there are no distinctive changes in the local dynamo convection process before excursions. Rather, it may be that the common factor is the diminished field intensity.

#### RELATIONSHIP OF LONG-TERM DIRECTIONAL PSV TO PALEOINTENSITY

**Figure 8** shows the relative paleointensity variability of the three Region C PSV records. There is a strong sense of centennialto millennial-scale variability in these records that is coherent among them. They have 47 features in common between 15 and 60 ka. The directional and paleointensity spectra of these three PSV records are summarized in Lund (2007). Furthermore, **Figure 9** shows the relative paleointensity variability between Regions C, E, and F. There appears to be a similar coherent pattern of variability among all three regions on at least a millennial scale. This is consistent with other estimates that paleointensity is a coherent process on a global scale (e.g., Channell et al., 2009).

The Regions C, E, and F paleointensity records in **Figure 9** have been averaged in 3- and 9-ka intervals just like the directional records in **Figures 6**, **7**. The results are shown in **Figure 10**. As expected, all three regions share a common pattern of paleointensity variability during the last glacial cycle. These three high-resolution paleointensity records can be compared

with the PISO-1500 compilation of Channell et al. (2009) in **Figure 11**, based on 13 stacked paleointensity records from around the World (average sediment accumulation rate of 11 cm/kyr). The 3–9 ka averaged results from **Figure 10** are consistent with the stacked PISO-1500 records (labeled features 2–9 in **Figures 10**, **11**).

There is a general sense in the paleomagnetic community (e.g., Merrill et al., 1998) that local directional field variability is inversely proportional to field intensity with low intensity associated with larger field variability and excursions. The three regions in **Figure 10** do show low paleointensity during the three excursion intervals (dashed lines) with anomalously high directional dispersion (**Figures 6**, **7**). Similarly, low directional dispersions in all three regions are clearly associated with high intensity intervals around 18 and 48 ka. This does suggest, that regional vector and paleointensity variability are

inversely correlated globally even though regional inclination and declination variability (noted above) are not. Correlation of inclination and paleointensity variability at three IODP sites In the Bering Sea (Lund et al., 2017b) for the entire Brunhes Epoch (last 780 ka) corroborate this pattern.

#### DISCUSSION

Traditional long-term PSV studies provide a sense of averaged field characteristics on a global scale (SHA or latitude/longitude transects) with averaging over ∼10<sup>6</sup> years. This study estimates averaged field characteristics within limited regions on timescales of 10<sup>3</sup> -10<sup>4</sup> years. The question is how do these estimates compare and what might we conclude about dynamo activity from their differences?

The PSV data from Region C (subtropical Atlantic Ocean) show strong directional variability on time scales of 10<sup>3</sup> -10<sup>4</sup> years. The scalar inclinations and declinations display ∼30◦

variations over 3-ka averaging and ∼15◦ variations over 9 ka averaging. The source of that variability must be dynamo activity, with most coming from dynamo activity directly below Region C in the outer core. The fact that all three PSV records in Region C show the same statistical pattern of long-term averaged field suggests that a similar degree of spatial coherence in dynamo activity operates in the outer core. The 9-ka averages are comparable to averaging, at any locality, the entire Holocene variability noted in many recent SHA studies (e.g., Constable et al., 2016). This indicates that the Holocene is not sufficient to characterize the entire space/time pattern of dynamo activity. It also suggests that there is significant regional-scale dynamo activity operating on a 10<sup>4</sup> years scale (as suggested by Lund et al., 1988). It may be that such PSV analyses may provide a new perspective on the actual space/time pattern of fluid flow and magnetic-flux regeneration that occurs regionally in the outer core.

The PSV data from Regions E and F corroborate the overall style of long-term directional field variability noted in Region C. It is likely that these three regions, together, reflect the overall long-term style of directional variability expected in any region on a global scale over 10<sup>3</sup> -10<sup>4</sup> years (e.g., Creer et al., 1983; Merrill et al., 1998; Lund, 2007). The detailed pattern of longterm directional variability for Regions C and E is consistent between them, but neither is consistent with Region F. Previous studies have suggested that directional variability in different regions should be different (e.g., Creer et al., 1983; Lund, 2007). So, the fact the Region F is different from Region C or Region E should not be a surprise. The surprise is that Regions C and E appear to be similar. One explanation for that is that the regional-scale dynamo activity is correlated between those two regions in the outer core. Some previous conceptual models of dynamo activity have focused on "convection rolls" as one source for dynamo activity; these convection rolls might operate in longitudinally-constrained swaths of the outer core between northern hemisphere and southern hemisphere high-latitude flux lobes (Bloxham and Gubbins, 1985). Regions C, D (work in progress), and E might all be associated with such a zone of coherent dynamo activity.

The long-term variability in field-vector dispersion is comparable in all three studied regions, even though the details of directional PSV are different in each region (or at least in Region F vs. Regions C and E). The long-term, 10<sup>5</sup> years average, vector dispersion for each region (**Table 1**) is comparable to 10<sup>6</sup> years averages (**Figure 2**), as well. But, the overall pattern of vector dispersion in each region is characterized by significantly lower vector dispersions for ∼75% of the time and significantly high vector dispersion for ∼25% of the time. This suggests that the vector dispersion on 10<sup>3</sup> -10<sup>4</sup> years time scales is much <10<sup>6</sup> years averages predict. This has been noted before (e.g., Lund, 1985; Lund et al., 2017b). It suggests that dynamo activity operates with a lower range of variability for most time, but is punctuated by shorter times of higher variability.

Long-term intensity variability on 10<sup>3</sup> -10<sup>4</sup> years timescales is largely synchronous in all studied regions. This has been noted previously in global studies of paleointensity (e.g., Channell et al., 2009). The intervals of lowest intensity (features 6 and 7 ′ /8 in **Figure 10**) occur at times of largest vector dispersion noted above. This relationship of higher (lower) directional field variability associated with low (high) paleointensity has been noted in other paleomagnetic studies (e.g., Merrill et al., 1998; Lund et al., 2016b). This indicates that vector and intensity variability are global in their pattern even though the details in style of vector variability differ from region to region. This suggests some kind of "connectedness" in the energy content of dynamo activity on a global scale, even though the detailed pattern of dynamo activity may be different from region to region. One speculation is that fluid-flow and flux-regeneration in each region generates toroidal field that connects regions and gives a coherence to the overall energy state of the core dynamo, but region-specific aspects of the fluid-flow and flux regeneration creates poloidal field that is largely region-specific and uncorrelated on a global scale.

The lowest intensity intervals with largest vector dispersion are also times of magnetic field excursions. True excursional directions do not always occur in all of our studied regions. But, enhanced vector dispersion does occur in all excursion intervals. (Remember, the vector dispersions were calculated after all true excursional directions with VGP latitudes <45◦ N were removed.) This suggests that directional excursions are a characteristic of low intensity intervals. It appears that high vector dispersion is also a characteristic of low intensity intervals. As such, it seems that low intensity is the fundamental feature and driving force (?) for excursions and high vector dispersions. In as much as low vector dispersion dominated for ∼75% of our study interval (last ∼70 ka), it may be that one can argue the field operates in a normal mode for most time and a separate, distinctive low-intensity (excursional) mode for a smaller portion of time. The low-intensity state may intrinsically yield a more complicated fluid-flow and flux regeneration pattern than highintensity intervals. This may produce even smaller regional-scale variability in both directions and intensity.

#### CONCLUSIONS

This study carried out a statistical analysis of high-resolution PSV records for the last ∼70 ka from three different regions of the Earth. We consider directional and intensity variability in each region on time scales of 10<sup>3</sup> -10<sup>5</sup> years in order to evaluate longterm PSV. We then compare those results with more traditional long-term PSV statistical studies averaged over ∼10<sup>6</sup> years.

Three replicate PSV records from one region (C, subtropical North Atlantic Ocean) were averaged at overlapping 3 and 9 ka intervals. Variability in both scalar inclination and declination variability and vector angular dispersion are significant and coherent among the three records. The vector dispersion is relatively low for most of the time but contains two relatively narrow intervals (∼25–45 and 60–65 ka) of high dispersion. (Vector dispersion in all records was calculated after removing directions with true excursional VGPs, VGPs < 45◦ N).

We have carried out a comparable statistical analysis on two other PSV records from other parts of the Earth (Region E, Chile margin; Region F, Philippines/Indonesia). The results for these two records are comparable to Region C in their overall style of variability. The scalar directional variability from Region F is quite different in detail from those of Regions C and E, as might be expected, but the scalar directional variability between is Regions C and E is remarkably consistent considering their distance from one another. This may be associated with them being on the same longitude swath and having some coherent dynamo activity occurring along that path.

#### REFERENCES


There is evidence for three magnetic field excursions in the studied interval—Mono Lake Excursion (∼34 ka), Laschamp Excursion (∼42 ka), and an unnamed excursion in OIS Stage 4 (∼61 ka). These three excursions all occur in the lowest intensity intervals and highest vector dispersion intervals in all three regions. It seems likely that low intensity intervals in dynamo activity are distinctly correlated with enhanced vector dispersion and occurrence of excursions around the World.

Three paleointensity records from Region C were subjected to the same statistical analysis as the directions. These records are also coherent in their pattern of variability. We also carried out a statistical analysis on the paleointensity records from Regions E and F. They also displayed the same pattern and coherence in variability that is noted in Region C. That similarity in paleointensity variability on a global scale is expected even though the detailed scalar directional variability is not coherent on a global scale. The pattern of intensity variability is strongly correlated with the pattern of vector dispersion on a global scale—high (low) intensity is associated with low (high) vector dispersion. The excursions all occur in intervals of the lowest paleointensity in all regions.

#### AUTHOR CONTRIBUTIONS

The author confirms being the sole contributor of this work and approved it for publication.

#### FUNDING

This work has been supported by several grants from the USA National Science Foundation (EAR1547605, NSF-OCE0962385, OCE-0082698, EAR-9980410, EAR-95 26940).


Brazilian North margin and paleoceanography of the western tropical Atlantic during the Late Quaternary. Palaeog. Palaeoclim. Palaeoecel. 415, 3–13. doi: 10.1016/j.palaeo.2014.05.030


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Lund. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

#### APPENDIX 1

Data quality required for high-resolution, long-term PSV studies of sediment sequences.

Minimum requirements:


Preferred additional requirements:


## Short-Term Magnetic Field Variations From the Post-depositional Remanence of Lake Sediments

Andreas Nilsson<sup>1</sup> \*, Neil Suttie<sup>1</sup> and Mimi J. Hill <sup>2</sup>

*<sup>1</sup> Department of Geology, Lund University, Lund, Sweden, <sup>2</sup> Department of Earth, Ocean and Ecological Sciences, University of Liverpool, Liverpool, United Kingdom*

Paleomagnetic records obtained from lake sediments provide important constraints on geomagnetic field behavior. Secular variation recorded in sediments is used in global geomagnetic field models, particularly over longer timescales when archeomagnetic data are sparse. In addition, by matching distinctive secular variation features, lake sediment paleomagnetic records have proven useful for dating sediments on various time scales. If there is a delay between deposition of the sediment and acquisition of magnetic remanence (usually described as a post-depositional remanent magnetization, pDRM) the magnetic signal is smoothed and offset in time. This so-called lock-in masks short-term field variations that are of key importance both for geomagnetic field reconstructions and in dating applications. Understanding the nature of lock-in is crucial if such models are to describe correctly the evolution of the field and for making meaningful correlations among records. An accurate age-depth model, accounting for changes in sedimentation rate, is a further prerequisite if high fidelity paleomagnetic records are to be recovered. Here we present a new method, which takes advantage of the stratigraphic information of sedimentary data and existing geomagnetic field models, to account for both of these unknowns. We apply the new method to two sedimentary records from lakes Kälksjön and Gyltigesjön where <sup>14</sup>C wiggle-match dating floating varve chronologies provide an independent test of the method. By using a reference magnetic field model built from thermoremanent magnetization data, we are able to demonstrate clearly the effect of post-depositional lock-in and obtain an age-depth model consistent with other dating methods. The method has the potential to improve the resolution of sedimentary records of environmental proxies and to increase the fidelity of geomagnetic field models. Furthermore, it is an important step toward fully explaining the acquisition of post-depositional remanence, which is presently poorly understood.

Keywords: paleomagnetism, geomagnetic field, post-depositonal remanent magnetization, lock-in, age-depth models

#### INTRODUCTION

Sediments record geomagnetic field variations through acquisition of a detrital remanent magnetization (DRM). Measurements of the acquired magnetization, i.e., paleomagnetic data, contain information about past geomagnetic field behavior on a range of time scales. Geomagnetic field models constructed from sedimentary paleomagnetic data and remanent magnetizations

#### Edited by:

*Christopher Davies, University of Leeds, United Kingdom*

#### Reviewed by:

*Andrew Philip Roberts, Australian National University, Australia Shuhui Cai, Scripps Institution of Oceanography, University of California, San Diego, United States*

> \*Correspondence: *Andreas Nilsson andreas.nilsson@geol.lu.se*

#### Specialty section:

*This article was submitted to Geomagnetism and Paleomagnetism, a section of the journal Frontiers in Earth Science*

> Received: *21 December 2017* Accepted: *03 April 2018* Published: *19 April 2018*

#### Citation:

*Nilsson A, Suttie N and Hill MJ (2018) Short-Term Magnetic Field Variations From the Post-depositional Remanence of Lake Sediments. Front. Earth Sci. 6:39. doi: 10.3389/feart.2018.00039* obtained from archeological artifacts and lava flows (hereafter referred to as archeomagnetic data) provide a global picture of the geomagnetic field and its evolution both at Earth's surface and at the core-mantle boundary (Korte et al., 2009, 2011; Licht et al., 2013; Nilsson et al., 2014). The superior geographical distribution offered by sedimentary paleomagnetic data compared to archeomagnetic data make sedimentary records essential for such global field reconstructions. Another advantage of sedimentary data is that they normally yield continuous records, often encompassing most of the Holocene, as opposed to archeomagnetic data that provide spot readings concentrated within the last 2000 years. Stratigraphic control in a sediment core also enables construction of precise age-depth models (e.g., Bronk Ramsey, 2009; Blaauw and Christen, 2011) using a variety of dating techniques. For example, varve counting can potentially provide annual resolution (e.g., Stanton et al., 2010; Striberger et al., 2011; Mellström et al., 2013).

Whilst sedimentary paleomagnetic data are useful, potential effects of post-depositional remanent magnetizations (pDRM) that can lead to both delay and smoothing of the recorded geomagnetic signal (e.g., Roberts and Winklhofer, 2004), present huge challenges for their application in geomagnetic field modeling and dating. The classical pDRM concept describes a process where the magnetization is acquired gradually as sediments compact until magnetic particles cannot further align with the geomagnetic field and become locked into place by the sediment matrix (Irving and Major, 1964; Kent, 1973; Hamano, 1980; Otofuji and Sasajima, 1981). This process has been described by so-called "lock-in" functions, which model the fraction of pDRM acquired as a function of depth below the surface or the mixed layer (e.g., Meynadier and Valet, 1996; Channell and Guyodo, 2004; Roberts and Winklhofer, 2004; Suganuma et al., 2011; Roberts et al., 2013a). Other efforts have questioned this classical pDRM acquisition concept arguing that flocculation of sediments prevents substantial postdepositional movement of grains within pore spaces (Katari et al., 2000). Alternative sediment mixing models have instead demonstrated the potential role of bioturbation in pDRM acquisition (Mao et al., 2014; Egli and Zhao, 2015). In addition, the presence of magnetotactic bacteria either living at the sediment/water interface or in the sediments, giving rise to either bio-depositional or bio-geochemical magnetizations, is becoming more recognized (e.g., Tarduno et al., 1998; Heslop et al., 2013; Roberts et al., 2013b; Larrasoaña et al., 2014; Mao et al., 2014). Despite uncertainties regarding the processes involved, there is mounting empirical evidence for natural remanent magnetization (NRM) acquisition at depths >10 cm below the sediment/water interface (Lund and Keigwin, 1994; Sagnotti et al., 2005; Suganuma et al., 2011; Mellström et al., 2015; Simon et al., 2018). These processes, hereby collectively referred to as lock-in, all act as a natural smoothing filter that produces a time lag between the age of the magnetization and the time of deposition of the sediments carrying the magnetization. This attenuation of high frequency geomagnetic field variations limits the possible resolution of paleomagnetic data, which depending on accumulation rates could mean effective resolutions on the order of 100–1000 years, which reduces the information content of the record (Roberts and Winklhofer, 2004).

Understanding the nature of lock-in is critical if geomagnetic field models, which rely increasingly on sedimentary data as models are extended back in time, are to reproduce higher frequency (centennial or even decadal) field variations. Nilsson et al. (2014) showed that the independent ages of sedimentary paleosecular variation (PSV) data used to constrain the global pfm9k.1a geomagnetic field model are systematically older than model predictions. This is particularly apparent over the past 3000 years where the model is more heavily constrained by archeomagnetic data, which suggests that post-depositional lockin might be a widespread phenomena in lake/marine sediments all over the world.

So far, investigations of pDRM effects in ancient sediments have required high-precision chronologies to be able to distinguish between lock-in delay and dating uncertainties. Mellström et al. (2015) were able to obtain such high-precision chronologies (uncertainties of ±20 years, 95% confidence) spanning the period between 3000 and 2000 years BP using the <sup>14</sup>C wiggle-match dating technique in two Swedish lake sediment records (Gyltigesjön and Kälksjön). By comparing PSV data from these records with archeomagnetic field model predictions that are not affected by post-depositional smoothing and delay, they investigated and evaluated the suitability of a range of lock-in filter functions proposed in the literature. Application of this approach is, however, limited by the high costs involved (15 and 19 radiocarbon dates used, respectively) and special conditions of having floating varve chronologies that provide roughly annual resolution between each <sup>14</sup>C date.

Here we present a new Bayesian method to simultaneously model lock-in delay and construct an age-depth model based on paleomagnetic data and archeomagnetic field model predictions. As well as being the first paleomagnetic method that addresses problems related to pDRM acquisition in age-depth modeling, the method provides key insights into pDRM acquisition. Most importantly, it is able to recover high frequency field variations that are lacking in Holocene field models (Nilsson et al., 2014). To illustrate these applications of our new method, we use the two Swedish lake sediment records, Gyltigesjön (Snowball et al., 2013) and Kälksjön (Stanton et al., 2010, 2011), as case studies.

#### METHODS

Our new Bayesian method models simultaneously lock-in and constructs an age-depth model based on paleomagnetic data. Details of the method are provided in the Supplementary Information. The age-depth model is built on the widely used "Bacon" dating software of Blaauw and Christen (2011) to combine both radiocarbon and paleomagnetic measurements. It is based on controlling sediment accumulation rates using an autoregressive gamma process. Radiocarbon age determinations R<sup>j</sup> taken along the sediment core at depth z<sup>j</sup> provide tie points to which the age-depth model can be fitted.

#### Paleomagnetic Data

Paleomagnetic measurements, inclination I<sup>i</sup> and declination D<sup>i</sup> , taken at depth z<sup>i</sup> can be used indirectly to constrain the age-depth model through correlation to a reference curve with known ages. Here we chose to include only paleomagnetic field directions and not relative paleointensity estimates, which are usually associated with larger and more poorly defined uncertainties.

#### Sedimentary Paleomagnetic Data Uncertainties

Uncertainties in sedimentary paleomagnetic directional data can be difficult to assess. The experimental error, e.g., the maximum angular deviation (MAD, see Kirschvink, 1980), typically only accounts for a fraction of the total error budget. Other potential sources of error include: (i) compaction leading to some degree of inclination flattening (Blow and Hamilton, 1978; Anson and Kodama, 1987), (ii) chemical and physical overprinting (Snowball and Thompson, 1990), (iii) compression and extension of sediments caused by coring equipment (Bowles, 2007) and sub sampling into u-channels or boxes (Gravenor et al., 1984). In addition, potential problems with core orientation and unwanted and unknown core rotation, which affects declination data, will also add to the error budget. Traditionally, the total error is estimated by smoothing data and calculating the dispersion of paleomagnetic directions over the selected range of depths/ages, preferentially from several independently oriented cores. This method has an obvious disadvantage for our analyses that it introduces additional data smoothing, which is not related to lock-in processes (Mellström et al., 2015). The error estimated using this approach could also underestimate data uncertainties severely if the record is only based on one core, which is often the case for Holocene paleomagnetic data.

To avoid binning or smoothing of data we use data at sample level. To account for experimental error, we convert the MAD into an α<sup>95</sup> angle, following Khokhlov and Hulot (2016), which in turn can be expressed as an α<sup>63</sup> angle, the Fisher (1953) distribution equivalent of a standard error, using:

$$
\alpha\_{63} = \frac{81}{140} \alpha\_{95}.\tag{1}
$$

To account for uncertainties related to coring and subsampling etc., we introduce an additional orientation error, which is added in quadrature to the α63. Based on Stanton et al. (2011), who measured the tilt of core penetration based on the angle of sediment laminations, we use an orientation error of 2◦ . This additional error effectively limits the measurement precision to ±2 ◦ but does not account for potential systematic errors, e.g., induced by core rotation. To minimize the influence of such systematic errors, both inclination and declination data are treated as relative values and adjusted so that the mean inclination (declination) matches the mean reference inclination (declination) determined over the same depths. Inclination and declination data are affected by different sources of uncertainty and are therefore treated independently. The inclination error is derived directly from the α<sup>63</sup> angle (σ<sup>I</sup> = α63) and the declination error, which depends on inclination, is determined as follows:

$$
\sigma\_{\rm D} = \frac{\alpha\_{63}}{\cos I}.\tag{2}
$$

#### Geomagnetic Field Reference Curve

To constrain an age-depth model, paleomagnetic data are correlated to a geomagnetic reference curve with "known" ages. Such reference time series, B(t) , are readily obtained from geomagnetic field model predictions for the site coordinates. Geomagnetic field model predictions frequently also come with some form of uncertainty estimates, σB(t), usually based on bootstrap sampling of the paleo/archeomagnetic data (Korte et al., 2009).

A range of both regional and global geomagnetic field models are available and the choice of which model to use largely depends on the age range under investigation. To account for potential lock-in effects in sediments, the choice of model is restricted to those based exclusively on archeomagnetic data to ensure that they are not affected by lock-in processes. This includes measurements on archeological artifacts and igneous rocks but not on lake or marine sediments. There are many archeomagnetic models from which to choose (e.g., Korte et al., 2009; Licht et al., 2013; Pavón-Carrasco et al., 2014) and as shown by Mellström et al. (2015), the choice of model has some implications for the results. Data limitations restrict the effective range of these models, and consequently application of the dating method, to the past 3–4 millennia.

The number of paleomagnetic data points used in our analyses must remain constant for all investigated age-depth models. This means that the ages associated with a paleomagnetic measurement are restricted to the age range of the geomagnetic field model used. We, therefore, introduce a paleomagnetic depth limit zlim, below which paleomagnetic measurements are discarded.

#### Outliers

The total error budget for a given paleomagnetic measurement, as outlined above, will account reasonably for errors arising from sampling and measurements as well as the model uncertainty, but some data may still have much greater dispersion. To account for the possibility of such outlying data, we follow a suggestion of Sivia (1996) who derived an expression for the likelihood of normally distributed data with uncertain error estimates (**Figure 1**). Take measurement y with an expectation µ and estimated uncertainty ω . We assume that the true uncertainty (ς) is greater or equal to our estimate (ω) and assign the prior distribution <sup>ω</sup> ς 2 , for ς ≥ ω. The unknown parameter, ς , is then marginalized to obtain the following likelihood:

$$L\left(\wp|\mu,\omega\right) = \frac{1}{\omega\sqrt{2\pi}}\frac{1-e^{-\frac{\chi^2}{2}}}{\chi^2},\tag{3}$$

where χ 2 = (y−µ) 2 ω2 . Although this form is only strictly valid for normally distributed (rather than spherically distributed) data, in the scheme used here declination and inclination are treated separately, and are assumed to be approximately normally distributed.

#### Lock-In Modeling

The lock-in process is modeled following Roberts and Winklhofer (2004). First the geomagnetic field input timeseries, B(t), is converted to a depth-series, B(z), using a given age-depth model (see Supplementary Information). A pDRM lock-in filter function, (z ′ ) (**Figure 2**, see below for details), describes the relative contribution of layer to z ′ the total pDRM acquired at each depth (z) down to a depth λ below the sediment/water interface, such that:

$$P\left(z\right) = \int\_0^\lambda B\left(z - z'\right) F\left(z'\right) dz',\tag{4}$$

where P (z)is the pDRM obtained by convolving the geomagnetic field depth-series. The primed coordinates refer to the depth during deposition and the unprimed coordinates represent the actual depth of the recovered sediments (compaction during burial is not taken into account). The z ′ variable takes values from 0, at the start of deposition of a sediment horizon, to λ as this horizon is progressively buried. The presence of a bioturbated or mixed layer is not taken into account in the model because this is not relevant to our case studies, where both records consist of undisturbed varved sediments (Mellström et al., 2015).

#### Filter Functions

The shape of the lock-in filter function remains a major uncertainty in pDRM modeling, because the process is still not well understood. Several lock-in functions have been proposed in the literature: exponential (Løvlie, 1976; Denham and Chave, 1982; Hyodo, 1984; Kent and Schneider, 1995; Meynadier and Valet, 1996; Roberts and Winklhofer, 2004), constant (Bleil and Dobeneck, 1999), linear (Meynadier and Valet, 1996; Roberts and Winklhofer, 2004), cubic (Roberts and Winklhofer, 2004), and Gaussian (Suganuma et al., 2011). These are shown in **Figure 2A** as the percentage of pDRM locked-in with depth, i.e., the cumulative functions C<sup>F</sup> z ′ . Apart from the Gaussian function, all of these lock-in functions rely on the assumption that the pDRM results from progressive consolidation and dewatering of sediments with gradual expulsion of interstitial water (Irving and Major, 1964; Kent, 1973; Hamano, 1980; Otofuji and Sasajima, 1981).

For the purpose of our analysis we constructed a simple parameterized description of the lock-in function, modified from Meynadier and Valet (1996):

$$F\left(z'\right) = \frac{\beta + 1}{\lambda^{\beta + 1}} \left(\lambda - z'\right)^{\beta},\tag{5}$$

with the corresponding cumulative function:

$$C\_F\left(z'\right) = 1 - \frac{\left(\lambda - z'\right)^{\beta + 1}}{\lambda^{\beta + 1}}\tag{6}$$

By varying one parameter, the lock-in shape factor β , we are able to reproduce several previously published functions, for example: constant (β = 0), linear (β = 1), and cubic (β = 3) (**Figure 2B**). Large β values produce lock-in functions with long tails, approaching an exponential function, where pDRM acquisition is mainly confined to the uppermost sediments. Negative β values (−1 < β < 0) instead lead to pDRM acquisition that increases with depth and then suddenly ceases at a depth λ . This latter behavior is physically not well justified but we retained such lockin functions in our analyses, acknowledging that we still do not understand the process.

#### Numerical Modeling

For numerical modeling we transform the depth domain to a regular grid B(z) => B(j). The pDRM output is calculated by convolving the geomagnetic input signal with a discretized lock-in function:

$$P\left(j\right) = \sum\_{k=0}^{L} \,^L B\left(j-k\right) F\_d(k),\tag{7}$$

where L is the lock-in depth (λ) in the discrete domain, which will depend on the resolution of the model (Fres). The discretized lock-in function F<sup>d</sup> k is determined from the cumulative filter function C<sup>F</sup> z ′ divided into L + 1 equal parts and calculating the difference between each part. By construction, F<sup>d</sup> (0) + F<sup>d</sup> (1) + . . . + F<sup>d</sup> (L) always sums to 1 and F<sup>d</sup> therefore performs acceptably for convolutions at small L approaching the grid resolution. For large L relative to the grid size, which is desirable, the benefit of using Fd(k) instead of F(k) is negligible. Convolution of the geomagnetic model prediction is performed in Cartesian coordinates Bx, By, B<sup>z</sup> and is then converted to obtain the convolved inclination and declination predictions(P<sup>I</sup> , PD).

The convolution with our different lock-in filters essentially acts as a running average, so the precision of the geomagnetic field prediction will generally increase with increasing λ . This

change in σ<sup>P</sup> with λ is not trivial to calculate but depends on (i) the shape of the lock-in shape function (β), (ii) geomagnetic field variations over the convolved range, and (iii) auto correlation of the geomagnetic field depth-series. We, therefore, chose a conservative estimate of no change in the magnitude of the error and simply convolved geomagnetic field model errors (σBI, σBD) according to (8). For the uppermost sediments, lockin is incomplete and it is necessary to normalize the convolved model error, which would otherwise go to zero at j = 0:

$$
\sigma\_P \text{ (j)} = \frac{\sum\_{k=0}^j B \left(j - k\right) F\_d(k)}{C\_F \text{ (j)}} \tag{8}
$$

for j < L.

#### The Lock-In Parameter Space

The lock-in parameter space comprises information on both the lock-in depth and the lock-in function. Since the full lockin depth λ depends strongly on β we chose to use the half lock-in depth l (Hyodo, 1984) instead, which corresponds to the depth at which 50% of the pDRM is acquired. The main difference between our lock-in filters (β ≥ 0) is the length of the tail, therefore, the effects on the convoluted data when using a "constant" (β = 0) or a "cubic" (β = 3) filter with the same l will be relatively similar. The actual lock-in depth λ is related to l according to:

$$\lambda = \frac{l}{1 - 2^{-\frac{1}{\beta + 1}}}.\tag{9}$$

The parameterized lock-in function covers a wide spectrum of possible shapes of lock-in filter functions. However, we are mainly interested in the range found in the literature, which corresponds roughly to β = 0 ("constant") to β = 6 (an approximation of the "exponential" filter cut-off at 99.9%). To explore this range more effectively, we introduce a transformation 0 < b < 1 (**Figure 2B**), from which β can be derived according to:

$$
\beta = \log\left(\frac{1}{b^2}\right) - 1.\tag{10}
$$

As shown in section Exploring Lock-in Parameters, a uniform prior distribution of b was chosen to obtain an elliptical posterior distribution of l and b (e.g., Sivia, 1996).

#### t-Walk MCMC Sampler

Following Blaauw and Christen (2011) we use a self-adjusting Markov Chain Monte Carlo (MCMC) algorithm, the t-walk (Christen and Fox, 2010), to explore the posterior distribution of the age-depth model as well as the lock-in depth and shape of the lock-in filter.

#### APPLICATIONS OF THE METHOD

First, we demonstrate the potential of our method as a dating tool in the absence of independent radiocarbon dates. The degree of smoothing exhibited by the paleomagnetic data is used as the sole means to constrain lock-in parameters. We then explore the more usual situation encountered where a few radiocarbon dates per record are available to see how well we can identify lock-in parameters and compare the results to the findings of Mellström et al. (2015). Finally, we demonstrate how the identified lock-in parameters could be used to recover high-frequency geomagnetic variations. All data used for these examples, both paleomagnetic and chronologic data, were retrieved from the GEOMAGIA50.v3 database (Brown et al., 2015).

#### Gyltigesjön Paleomagnetic Age-Depth Model

The PSV record from lake Gyltigesjön, located in southwestern Sweden, consists of paleomagnetic measurements on 546 discrete samples collected over a 8-m-long sediment sequence from six partly overlapping cores (Snowball et al., 2013). The sequence has been dated using <sup>137</sup>Cs activity and 22 radiocarbon determinations, 15 of which were used for the wiggle-match dating (Mellström et al., 2013).

As noted in the original publication, core rotation can be a problem in the top meter of the cores. After further inspection of the data, we suspect this may extend down to 2 m of each core and, therefore, exclude the declination data from these parts in our analysis. In addition, we also exclude inclination data from the top 1 m in core GP4 (see Snowball et al., 2013). Excluded data are shown in **Figures 3B,C** as open symbols (without error bars).

Based on the findings of Mellström et al. (2015), we use the ARCH3k\_cst.1 model (Korte et al., 2009), constrained by a selection of high quality archeomagnetic data, as our input geomagnetic reference curve. The ARCH3k\_cst.1 model does not have any uncertainty estimates, so we assume a constant inclination error σBI = 2 ◦ and obtain the declination error σBD from (6) (see Mellström et al., 2015). The model is extrapolated to the year the cores were retrieved (2010 AD) and beyond, assuming no change in the field, which is necessary when considering that lock-in may still be on-going in the uppermost sediments. We also use the full, but less supported, age-range of the model, extending back to 2000 BC, to maximize the number of paleomagnetic measurements included in the analysis (zlim = 550 cm).

The parameters used for the age-depth model are summarized in **Table 1.** For more details, see Supplementary Information and Blaauw and Christen (2011). Sediment accumulation rates in the upper part of the core investigated here are expected to vary between c. 3 and 10 years/cm (Snowball et al., 2013), so we set our prior for the accumulation rate to a gamma distribution with shape parameter 1.5, recommended by Blaauw and Christen (2011), and mean of 5 years/cm. We do not expect drastic accumulation rate changes and, therefore, set the prior for accumulation rate variability as a beta distribution with memory strength 20 and mean 0.9 ("high memory"). We divide the sequence into 25 slices, which results in (corresponding to zlim/25). Based on the results of Mellström et al. (2015), we set lmax = 50 cm and the resolution for the discrete convolution Fres = 1 cm.

The intention with our new method is to combine all available chronological information in one age-depth model. However, for comparison purposes we present two separate agedepth models and compare them to the independent wigglematch dated chronology. In **Figure 3A** we show the posterior distribution of the age-depth model based on the paleomagnetic data only. For reference we show the 95% confidence limits (dotted black lines) of an age-depth model based on all (six) available macrofossil radiocarbon determinations for the investigated depths (Snowball et al., 2013), but without using paleomagnetic data. Both age-depth models agree well within their respective uncertainties and also overlap reasonably well with the wiggle-matched chronology between 3000 and 2200 years BP (solid cyan line). Importantly, there are no apparent systematic offsets between the paleomagnetic- and radiocarbonbased chronologies.

To visualize the model-data comparison we plot paleomagnetic data and the convolved geomagnetic field model prediction (**Figures 3B,C**) on a "preferred" age-depth model, determined as the weighted mean age for each depth, using the "preferred" l and b (determined from the posterior density in **Figure 3F**) to derive λ and β for the model convolution. The comparison between the original model prediction (red lines) and the convolved model (blue lines) shows the smoothing and temporal offset (c. 200 years) induced by lock-in. If not accounted for, this would lead to systematically young ages in the paleomagnetic-based chronology (**Figure 3A**). Since we do not include any hard temporal constraints apart from the top of the sequence, the lock-in depth must be constrained mainly by the degree of smoothing in the paleomagnetic data compared to the geomagnetic field model prediction.

The posterior for accumulation rate is similar to the prior distribution (**Figure 3D**), but the posterior for the memory (**Figure 3E**) indicates that the data require more accumulation rate variability than suggested by the prior distribution. The posterior distribution of the lock-in parameters (**Figure 3F**) suggests a half lock-in depth l of 30.7 ± 6.5 cm (mean and standard deviation) and a lock-in shape β > 0, i.e., negative values of β are rejected.

#### Kälksjön Paleomagnetic Age-Depth Model

Lake Kälksjön is located in west central Sweden. The PSV record consists of paleomagnetic measurements on 940 discrete samples from six partially overlapping piston cores taken from a 7-m-long sequence (Stanton et al., 2010, 2011; Mellström et al., 2015). The upper 2 m of the record (relevant for this study) have been dated by a combination of varve counting, <sup>137</sup>Cs activity, Pb pollution dating, 22 bulk radiocarbon dates, 19 of which were used for the wiggle-match dating, and four macrofossil radiocarbon determinations (Mellström et al., 2015).

Similar to Gyltigesjön, we exclude declination data from the uppermost 75 cm of each core from our analyses, due to suspected core rotation. We use the same geomagnetic field input model and set zlim = 200 cm. The prior for the accumulation rate is set to a gamma distribution with shape 1.5 and mean 20. As above, we divide the record into 25 slices, 1c = 8 cm (corresponding to zlim/25), but due to the smaller value of we allow for more accumulation rate variability with memory strength 10 and mean 0.8. Based on the results of Mellström et al. (2015), we set lmax = 25 cm and the resolution for the discrete convolution to Fres = 0.5 cm (**Table 1**).

We produce a paleomagnetic-based age-depth curve (ignoring all radiocarbon dates) and one radiocarbon-based

ARCH3k\_cst.1 geomagnetic field model prediction (solid red) filtered with the preferred lock-in function (solid blue) with 1σ uncertainty envelope. Rejected data are shown as open black symbols. In the lower panels: The prior (red) and posterior (gray) distributions of (D) accumulation rate and (E) memory. (F) The posterior distribution of the half lock-in depth (*l*) and lock-in shape(β), where the blue circle marks the "preferred" parameters, determined by the grid cell with highest probability.

age-depth curve, using all (four) macrofossil determinations, and compare these to the independent wiggle-matched chronology (**Figure 4A**). As was the case for Gyltigesjön, all three chronologies are consistent within their errors and do not contain any evidence for systematic errors. We note that the paleomagnetic-based age-depth models for both Gyltigesjön and Kälksjön indicate a distinct accumulation rate change around 2600 years BP, which suggests inconsistencies between paleomagnetic data and geomagnetic field model predictions and/or unsuitable lock-in filter functions. The posterior distribution of the lock-in parameter space (**Figure 4F**) suggests a half lock-in depth l of 16.6 ± 3.3 cm (mean and standard deviation), with a weakly constrained lock-in shape, again rejecting negative β values.

Also shown in **Figure 4A** (green line) is the "corrected" varve chronology from Stanton et al. (2010). Based on Pb pollution data and PSV correlations, Stanton et al. (2010) found that the original varve chronology was missing ∼270 varves in the younger part of the record. The varve chronology was, therefore, corrected by evenly distributing the missing varves over the part spanning the last 1000 years. However, the wiggle-match dating by Mellström et al. (2015) suggests that an additional 213 varves may be missing further down the core. In **Figure 5**, we plot the Kälksjön <sup>206</sup>Pb/207Pb data on our paleomagnetic-based

#### TABLE 1 | Age-depth model parameters.


FIGURE 4 | (A) The age-depth posterior distribution for Kälksjön based on paleomagnetic data (yellow/red shading) with the 95% confidence envelope (dashed black) compared to the wiggle-match chronology (solid cyan), the adjusted varve chronology (solid green) and the 95% confidence envelope of the radiocarbon-based chronology (dotted black). Calibrated probability density curves of available macrofossil radiocarbon dates are shown for reference as green shaded areas. The blue line represents the preferred age-depth model based on the weighted mean age for each depth. (B) Inclination and (C) declination data (filled black circles) with 1σ error compared to the ARCHk\_cst.1 geomagnetic field model prediction (solid red) filtered with the preferred lock-in function (solid blue) with one sigma uncertainty envelope. Rejected data are shown as open black symbols. In the lower panels: The prior (red) and posterior (gray) distributions of (D) accumulation rate and (E) memory. (F) The posterior distribution of the half lock-in depth (*l*) and lock-in shape (β), where the blue circle marks the "preferred" parameters, determined by the grid cell with highest probability.

chronology and compare them with three independently dated records of Pb pollution from different parts of Sweden (Brännvall et al., 1999; Bindler et al., 2009), one of which (Koltjärn) was used originally by Stanton et al. (2010). Our paleomagnetic-based chronology places the Greek-Roman Pb peak (corresponding to a decline in <sup>206</sup>Pb/207Pb ratios around 2000 years BP) roughly 100 years too early but well within the chronological error. Otherwise the data are in good agreement, which demonstrates that the paleomagnetic-based chronology is consistent with both the wiggle-match chronology and Pb pollution data.

also plotted.

similar independently dated records from Kassjön, Koltjärn, and Kalven are

### Exploring Lock-In Parameters

To further constrain the lock-in parameters, we produced a third age-depth model for Gyltigesjön and Kälksjön, respectively, including both paleomagnetic data and all macrofossil radiocarbon determinations (see Supplementary Information). Additional temporal constraints, provided by radiocarbon determinations, do not lead to any significant changes to the paleomagnetic-based age-depth curves but provide more information on the lock-in parameters (**Figures 6A,B**). To compare our results to those of Mellström et al. (2015), we repeated their analysis using our new parameterized lock-in function (**Figures 6C,D**) with the same paleomagnetic data (binned every 7 cm and 4 cm, respectively) and Gaussian error distribution as in the original study (for a more detailed description, see Supplementary Information). Even though the two analyses were conducted over different age ranges, with variable or fixed wiggle-matched timescales and with discrete or binned data, the results still look similar. In all four cases, the posterior distribution of the lock-in parameters is centered around β = 1 (i.e., a "linear" lock-in function) with l ≈ 30 cm and l ≈ 15 cm, respectively.

#### Recovering High-Frequency Field Variations

To demonstrate how the lock-in parameters, identified with our method, can be used to recover high-frequency field variations we use them to deconvolve the entire Gyltigesjön PSV record. Deconvolution is sensitive to noise, so to attempt to isolate a robust signal we fit the inclination and declination data with cubic smoothing splines weighted by data uncertainties. The spline fits were then transformed to Cartesian coordinates (assuming a unit sphere) and deconvolved by division in the frequency domain using the lock-in filter with highest probability from **Figure 6C** (β = 1.3, l = 28). The results are shown in **Figure 7** on the original published age-depth model, including the wiggle-matched dating (Snowball et al., 2013) and compared to ARCH3k\_cst.1 and pfm9k.1b (Nilsson et al., 2014) model predictions. The recent part of the deconvolved record (<3500 years) is essentially constrained to fit ARCH3k\_cst.1, although with a few clear differences such as the enhanced eastward declination swing around 2650 years BP (known as feature "f "; Turner and Thompson, 1981) and the shallower inclination around 2250 years BP. The older and independent part of the "restored" PSV record, which is the main gain of this approach, indicates that it is possible to extract a similar range of field variability in sedimentary data as in archeomagnetic data.

### DISCUSSION

#### The Lock-In Function

The method presented here enables the study of variations in the lock-in process in many lake or marine settings where expensive wiggle-match dating or other high-precision dating methods are either not available or applicable. Our analyses reproduce the results of Mellström et al. (2015), which supports the conclusion that the "linear" lock-in function is the most appropriate of available lock-in functions for the Gyltigesjön and Kälksjön records. In general, however, we conclude that the method does not distinguish much between variations in lock-in shape, but is mainly sensitive to the half lock-in depth, which is deeper in Gyltigesjön than in Kälksjön. The physical cause of pDRM in these sediments remains unclear but may be related to the unusually high organic content, in particular for Gyltigesjön (Mellström et al., 2015).

The parameterized lock-in function presented in this paper is able to reproduce most lock-in functions proposed in the literature. Although we have made no attempt to explain the underlying physical process, we propose that using this or other similar functions is better than ignoring lock-in altogether when considering paleomagnetic data for geomagnetic field modeling or for dating. Our new method provides a platform that could be used to experiment with other, more physically motivated (e.g., Egli and Zhao, 2015), variations of filter functions, as long as they can be defined with a limited number of parameters.

#### Age-Depth Model Applications

Although the method presented in this paper is primarily intended to be used to explore pDRM effects in sediments, it can also be used as a complement to conventional radiocarbon dating for sedimentary archives, particularly in areas where suitable materials for dating are scarce. The uncertainty of age-depth models based on paleomagnetic data will, however, only be as good as the regional geomagnetic field reference curve. When considering lock-in effects, the quality of the geomagnetic field model predictions depends on the availability of archeomagnetic data, which effectively restricts the method to sites in the northern hemisphere and the past 3–4 millennia. However, we note that the method could potentially also

be used with a wider range of geomagnetic field models based on both archeomagnetic and sedimentary data, such as pfm9k.1b, with the caveat that temporal offsets caused by lock-in processes may not be appropriately adjusted for. This would increase the potential application of the dating method in time, beyond the archeomagnetic models, and to the southern hemisphere.

#### Implications for Geomagnetic Field Reconstructions

As noted above, both paleomagnetic age-depth models for Gyltigesjön (**Figure 3**) and Kälksjön (**Figure 4**) contain a change in accumulation rate around 2600 years BP, which is not supported by the wiggle-match timescales (Mellström et al., 2015). It appears that the age-depth models require high accumulation rates prior to 2600 years to minimize lock-in smoothing of the geomagnetic field model prediction over this time interval and to better capture the large amplitude eastward declination swing (feature "f ") predicted by the data. This implies that the geomagnetic field model prediction might be too smooth at around this time, i.e., that the declination swing was more pronounced than the model predicts, which is supported by comparison to the deconvolved Gyltigesjön record (**Figure 7**). Likewise, the shallow inclination predicted by the deconvolved record at around 2250 years BP could also be a real feature of the field that is not captured by the archeomagnetic model.

These observations highlight the potential limitations of archeomagnetic field models, which lack data from high latitude locations such as the lake sites, and why it is important to include sedimentary data in global geomagnetic field models. By including assumptions about lock-in processes we have shown that it is possible to recover high-frequency variations from sedimentary paleomagnetic data and remove potential systematic age uncertainties (Nilsson et al., 2014). In **Figure 7**, we demonstrate this by deconvolving the Gyltigesjön record. However, with such an approach it will always be difficult to isolate the real signal and avoid amplifying noise. We propose that a better method would be to incorporate the full agedepth method presented here into the modeling, and to fit the sedimentary data with a lock-in filtered model prediction (Nilsson and Suttie, 2016). With a correct understanding of the

smoothing effect of the lock-in process, e.g., using updated lockin filter functions, it should be possible to recover the higher frequency geomagnetic field components that are currently captured by archeomagnetic data throughout the Holocene.

#### CONCLUSIONS

We have presented a novel Bayesian method to simultaneously model lock-in processes and construct an age-depth model based on paleomagnetic data and archeomagnetic field model predictions. The method provides, for the first time, a way to constrain the effects of pDRM acquisition in sediments in the absence of high-precision dating methods. This is the first paleomagnetic method that addresses problems related to lockin delay in age-depth modeling and has potential to be used to complement conventional radiocarbon dating for sedimentary archives. Finally, we have demonstrated how the identified lock-in parameters can be used to recover high-frequency geomagnetic field variations from sedimentary paleomagnetic data by deconvolving the Gyltigesjön Lake PSV record which can lead to an increase in the fidelity of geomagnetic field models.

The deconvolved Gyltigesjön PSV data shown in **Figure 7** are provided in the EarthRef.org Digital Archive (ERDA, http://www.earthref.org) by searching for Gyltigesjön or the title of this study.

#### AUTHOR CONTRIBUTIONS

AN wrote main manuscript with contributions from both NS and MH. AN and NS did most of the work on developing the methodology.

#### ACKNOWLEDGMENTS

This work was funded by the Swedish Research Council (2014- 4125), the Crafoord Foundation (20150843), and the Natural Environment Research Council, UK (NE/I013873/1). We thank Andrés Christen and Maarten Blaauw for providing helpful advice regarding the Bacon implementation. Andrew P. Roberts and Shuhui Cai are thanked for their thorough reviews, which substantially improved the manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2018.00039/full#supplementary-material

#### REFERENCES


pollution history and statistical correlation. Quat. Geochron. 5, 611–624. doi: 10.1016/j.quageo.2010.03.004


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Nilsson, Suttie and Hill. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.