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

# Archaeointensity record of weak field recurrence in Japan: new data from Late Yayoi and Kofun ceramic artefacts

E. Tema,<fnref n="1" />,<fnref n="2" /> Y. Santos,<fnref n="3" /> R. Trindade,<fnref n="4" /> G.A. Hartmann,<fnref n="5" />,<fnref n="6" /> T. Hatakeyama,<fnref n="6" /> F. Terra-Nova,<fnref n="7" /> N. Matsumoto,<fnref n="8" /> J. Mitsumoto<fnref n="9" /> and M. Gulmini<fnref n="3" />

Accepted 2022 December 10. Received 2022 December 10; in original form 2022 October 10

## SUMMARY

We present new absolute archaeointensity data from six archaeological sites situated in the Okayama Prefecture, Japan. The materials studied are well-dated fragments from pottery, ceramic coffins and haniwa artefacts. Their ages range from 160 AD to 675 AD, covering the Late Yayoi and Kofun periods. Rock magnetic experiments suggest the presence of magnetite and/or Ti-magnetite as the main carrier of the remanence, with a possible minor contribution of higher coercivity minerals. After thermal demagnetization experiments, the most magnetically stable samples were selected for archaeointensity analysis performed following the double-heating method proposed by Thellier and modified by Coe. Partial thermoremanent magnetization (pTRM) checks and pTRM tail-checks were performed for monitoring possible chemical alterations during heating. All measurements were corrected for both anisotropy and cooling-rate effects. Successful archaeointensity determinations, following rigorous selection criteria, were obtained for samples from all the investigated archaeological sites. Compared with literature data from Japan, the new high-quality data show significantly lower intensity values. They also reveal possible fast secular variation changes during the Late Yayoi period and very weak geomagnetic intensity field around 630 AD. Such values offer evidence of a possible recurrence of weak intensity field in East Asia, suggesting an ancient recurrence of the West Pacific Anomaly. The new data might change the archaeomagnetic field models interpretations in the area, even though more data are still necessary to better understand the secular variation in Japan and the temporal evolution of the geomagnetic field's behaviour in East Asia.

**Key words:** Japan; Archaeomagnetism; Palaeointensity; Palaeomagnetic secular variation.

## 1 INTRODUCTION

The Earth's magnetic field of internal origin is continuously changing in time, undergoing both long- and short-term variations. Such variations are of great interest since they represent the main source of observational information about the past behaviour of the geomagnetic field and about the outer's core flow magnitude and pattern. Several studies based on archaeomagnetic data showed that the geomagnetic field has been characterized by untypical episodes

of very fast- and short-lived local changes in both direction and intensity (e.g. Shaar et al. 2016; Cai et al. 2017; Osete et al. 2020), which are not observed in the present field (Alken et al. 2021). The origin, duration and geographical occurrence of such extreme geomagnetic field variations are subject of open discussion (Livermore et al. 2014; Davies & Constable 2018; Korte & Constable 2018; 2020; 2021) and their behaviour, feasibility and interpretation in terms of core dynamics, remain very challenging. To better understand and constrain the ancient geomagnetic field behavior, the

<footnote n="1">¹Dipartimento di Scienze della Terra, Università degli Studi di Torino, via Valperga Caluso 35, 10125 Torino, Italy. E-mail: evdokia.tema@unito.it</footnote>

<footnote n="2">²ALP-CIMAN Alpine Palaeomagnetic Laboratory, via Luigi Massa 6, 12016 Peveragno, Italy</footnote>

<footnote n="3">³Dipartimento di Chimica, Università degli Studi di Torino, via Pietro Giuria 7, 10125 Torino, Italy</footnote>

<footnote n="4">⁴Departamento de Geofísica, Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo, Rua do Matão, 1226, 05508-090, São Paulo, Brazil</footnote>

<footnote n="5">⁵Instituto de Geociências, Universidade Estadual de Campinas, Rua Carlos Gomes, 250, 13083-855, Campinas, Brazil</footnote>

<footnote n="6">⁶Faculty of Humanities and Social Sciences, Okayama University, Kita-ku, Okayama 700-8530, Okayama, Japan</footnote>

<footnote n="7">⁷Laboratoire de Planétologie et Géosciences, CNRS UMR 6112, Nantes Université, Université d'Angers, Le Mans Université, 2 rue de la Houssinière, F-44000 Nantes, France</footnote>

<footnote n="8">⁸Research Institute for the Dynamics of Civilizations, Okayama University, Kita-ku, Okayama 700-8530, Okayama, Japan</footnote>

<footnote n="9">⁹Graduate School of Humanities and Social Sciences, Okayama University, Kita-ku, Okayama 700-8530, Okayama, Japan</footnote>
study of very well-dated archaeological clays and volcanic eruptions is necessary. Such archaeomagnetic records offer important information and are essential for the construction of global geomagnetic field models and for core dynamics assessments.

Nowadays, several geomagnetic models are available but their reliability and capacity to detect past geomagnetic field changes and to reliable reconstruct the fine characteristics of the geomagnetic field behaviour, highly depend on the quality and the temporal and spatial distribution of the reference data (Constable et al. 2016; Campuzano et al. 2019; Pavón-Carrasco et al. 2021). Unfortunately, obtaining a uniform geographical distribution of the archaeomagnetic data is a very difficult (if not impossible) task. Presently, several areas are poorly covered by data either due to the limited availability of suitable materials (e.g. inhabited areas such as deserts, isolated islands and high mountains and/or difficult accessible areas such as war and conflict zones) or to the lack of scientific and/or economical resources for palaeomagnetic analyses (Di Chiara et al. 2020). The quality of the reference data is also an important issue. Even though many quality criteria have been proposed to select the most reliable palaeomagnetic records, the quality of the data, both in terms of dating and experimental determination, is still one of the main challenges faced by the palaeomagnetic community. Such problem is particularly important when dealing with palaeointensity records which are obtained through complex laboratory procedures. Palaeointensity data are thus often characterized by low success level in respect to the directional determinations, significantly influencing the reliability of the geomagnetic field models and their ability to detect short and abrupt intensity variations (Pavón-Carrasco et al. 2014; Livermore et al. 2021).

Some of the very first archaeomagnetic studies in Japan were published by Watanabe (1958, 1959) who reported directional data from a very large number of baked clay archaeological structures and proposed the first secular variation curve for chronological dating in Japan. Following, Nagata et al. (1963) published data from basaltic lavas and baked-clay potteries covering the last 5000 yr. Since then, several other studies enriched the directional and archaeointensity records in Japan, based on both archaeological material and volcanic rocks (e.g. Sasajima 1965; Kitazawa 1970; Sakai & Hirooka 1986; Yu 2012; Yamamoto et al. 2015). Among them, some of the most known and widespread archaeodirection studies were published by Hirooka (1971, 1983), who measured the archaeodirection of over 220 well-dated baked kilns, calculating a directional secular variation curve which is even now used for dating in Japan. However, even though the directional experimental protocols have only slightly changed over time, continuous research on archaeointensity has introduced substantial advances, demonstrating the need of strict and rigorous protocols to obtain reliable absolute intensity data. In Japan, the majority of the available up to now archaeointensity data come from studies carried out more than 20 yr ago, and thus important advances on the archaeointensity determination protocols, such as cooling rate (CR) and anisotropy corrections, were not considered. According to the GEOMAGIA50.v3.4 database (Brown et al. 2021), since 2000 only two studies from archaeological material (Kitahara et al. 2018; 2021) and four studies from volcanic rocks (Yoshihara et al. 2003; Mochizuki et al. 2004; Yamamoto & Hoshi 2008; Yu 2012) have been published. Kitahara et al. (2021) recently presented new archaeointensity results from archaeological relics in Osaka Prefecture, obtained with the Tsunakawa-Shaw method (Yamamoto et al. 2015). The comparison of these new data with previous data from Japan showed a significant discrepancy with the much higher archaeointensities

published before 1986. Such older data were mainly obtained with the conventional Thellier method and its variants, but they missed experimental corrections currently used. This underlines the need of obtaining new high-quality archaeointensity data for Japan in order to reconstruct a reliable geomagnetic field intensity curve in the past and contribute to the improvement of the geomagnetic field models prediction in the area.

In this perspective, we present here new archaeointensity data from six archaeological sites, situated in the Okayama Prefecture, Japan. The material studied come from well-dated pottery sherds, haniwa figures (ceremonial earthenware made to be placed on the Kofun burial mounds) and coffin ceramic fragments, with ages ranging from 160 CE to 675 CE. The new data, together with previously published archaeomagnetic records and global geomagnetic field models are used to better understand the Earth's magnetic field variation in Japan during the last 2000 yr. Evidence of sharp intensity variations during the Late Yayoi period and implications on weak field recurrence in East Asia are also discussed.

## 2 MATERIALS AND METHODS

Samples come from the archaeological sites of Tatezaka (TATEZ), Tatetsuki (TATET), Tenguyama (TEN), Nima Ohtsuka (NIMA), Sada Higashizuka (SADA1) and Sada Nishizuka (SADA2 and SADA4), all of them situated at the Okayama Prefecture, Japan (Fig. 1). From the Sada Nishizuka site, samples from two ceramic coffins were collected, coffin no 2 (SADA2) and coffin no 4 (SADA4). From each site (or collection), 7–10 samples were selected belonging to fragments from pottery, baked clay coffins and haniwa artefacts. All samples are very well-dated based on archaeological evidence such as stratigraphy, pottery style and coffin production (Table 1). Their ages vary from the Late Yayoi to the Kofun periods, covering almost six centuries from 160 CE to 675 CE. This period is of great cultural relevance in Japan, characterized by increasing social stratification. Such social differentiation into hierarchically superposed classes is reflected in the monumental burial mounds constructed for persons of the upper social layers. Information about the archaeological sites, the studied samples and their dating is given in the Supporting Information and is summarized in Table 1.

For each independent fragment and depending on its dimensions, 3 to 7 small cubes of about 1 cm edge were prepared and used for rock magnetic and archaeointensity experiments. Measurements were carried out at the palaeomagnetic laboratories of ALP-CIMAN Alpine Palaeomagnetic Laboratory (Peveragno, Italy), Okayama University of Science (Okayama, Japan) and University of São Paulo (São Paulo, Brazil).

The magnetic mineralogy of the samples was investigated through isothermal remanent magnetization (IRM) acquisition curves and thermal demagnetization of a composite three-axis IRM component (Lowrie 1990). The IRM curves were stepwise obtained up to 5 T with a Quantum Design MPMS XL-5 SQUID Magnetometer, in the Research Instruments Center, Okayama University of Science. Lowrie experiments were obtained by imparting a 1.2 T direct field along the z-axes of the samples, a 0.6 T field along the y-axes and a 0.16 T along the x-axes with an ASC Scientific pulse magnetizer; then, samples were thermally stepwise demagnetized up to 540 °C with a TD-48 SC (ASC Scientific) furnace and their remanent magnetization was measured with a JR6 Spinner magnetometer (Agico).
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Map of Japan with the location of the Okayama prefecture; (b) photo of the samples collected from the Tenguyama (TEN) site and (c) details of two fragments from the Tatetsuki site (TATET).</figcaption></figure>

<table><caption>**Table 1.** Information about the studied archaeological sites and their dating. Columns: Archaeological site; Code; Latitude (°); Longitude (°); Number of independent fragments; Type of ceramic samples; Archaeological age based on stratigraphic evidence and pottery typology.</caption><thead><tr><th>Archaeological Site</th><th>Code</th><th>Lat (°)</th><th>Long (°)</th><th>N</th><th>Type</th><th>Archaeological age</th></tr></thead><tbody><tr><td>Tatetsuki</td><td>TATET</td><td>34.66</td><td>133.82</td><td>8</td><td>Pottery</td><td>160–190 AD</td></tr><tr><td>Tatezaka</td><td>TATEZ</td><td>34.66</td><td>133.67</td><td>7</td><td>Pottery</td><td>160–190 AD</td></tr><tr><td>Tenguyama</td><td>TEN</td><td>34.62</td><td>133.72</td><td>10</td><td>Haniwa</td><td>470–500 AD</td></tr><tr><td>Nima Ohtsuka</td><td>NIMA</td><td>34.62</td><td>133.71</td><td>8</td><td>Haniwa</td><td>550–600 AD</td></tr><tr><td>Sada Higashizuka</td><td>SADA1</td><td>34.94</td><td>133.63</td><td>9</td><td>Coffin</td><td>630 AD</td></tr><tr><td>Sada Nishizuka, coffin no 2</td><td>SADA2</td><td>34.94</td><td>133.63</td><td>7</td><td>Coffin</td><td>630–650 AD</td></tr><tr><td>Sada Nishizuka, coffin no 4</td><td>SADA4</td><td>34.94</td><td>133.63</td><td>7</td><td>Coffin</td><td>650–675 AD</td></tr></tbody></table>

The sample's magnetic stability was analyzed by both stepwise thermal (THD) and alternating field (AFD) demagnetization. Two fragments per site (a total of 14 specimens) were AF demagnetized up to 180 mT with an automated AF demagnetizer and spinner magnetometer (Natsuhara-Giken Dspin). Then, one specimen from all fragments (a total of 54 specimens) was stepwise thermally demagnetized with a TD-48 SC (ASC Scientific) furnace while the magnetic remanence was measured with a JR6 Spinner magnetometer (Agico). After each demagnetization step, the bulk magnetic susceptibility was measured with a KLY3 Kappabridge (Agico) to monitor possible magnetic mineralogy changes during heating.

Archaeointensity experiments were performed using the classical Thellier & Thellier (1959, TT) double-heating palaeointensity method, as modified by Coe (1967), with first zero field (Z) and then in-field (I) heating-cooling steps (TT-ZI protocol). At least three specimens per fragment were used for archaeointensity determinations, while one more fresh specimen was used for the CR experiments. All archaeointensity analysis were performed at the Palaeomagnetic Laboratory of the University of São Paulo (USP-Mag). Experiments were carried out in stepwise increasing temperature steps between 100 and 575 °C, using an ASC Scientific single-chamber oven with an inducing field coil. A laboratory field of 60 µT was applied parallel to the z-axis of the specimen. Heating-cooling cycles were performed in air during 30 min each for both zero-field and in-field steps. Thermochemical alteration was monitored with additional partial thermoremanent magnetization (pTRM) checks carried out every two experimental temperature steps (Coe et al. 1978). Five additional checks were performed at 200, 300, 350,

400 and 500 °C (pTRM-tail checks; Riisager & Riisager 2001) to monitor multidomain effects. The anisotropy of thermoremanent magnetization (ATRM) was investigated at both 350 and 500 °C through pTRM acquisition along six positions (X, -X, Y, -Y, Z, -Z), following the procedure described by Hartmann et al. (2011) and Poletti et al. (2016). The anisotropy correction factor was calculated according to Veitch et al. (1984). The CR effect was also investigated using an additional fresh specimen. For the CR correction factor calculation, specimens were heated and cooled in an applied field of 60 µT for three times: once following a rapid cooling of around 30 min, then following a slow cooling of around 12 hr and finally repeating one more rapid cooling (of 30 min). The slow cooling step was performed introducing thermal insulators into the oven and turning off the cooling fan. Further details on the experimental procedures can be found in Hartmann et al. (2019).

# 3 RESULTS

## 3.1 Magnetic mineralogy

IRM acquisition curves obtained for one sample per site indicate the dominance of a low-coercivity magnetic mineral, often accompanied by a minor high-coercivity magnetic component (Fig. 2). The low-coercivity component is saturated at applied fields of around 0.2–0.4 T while application of a direct field up to 5 T clearly shows that in some samples, mainly from the NIMA and TATET sites, a high-coercivity component is also present (Figs 2c and d).
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Stepwise IRM acquisition curves up to 5 T for representative samples from (a) TEN; (b) SADA1; (c) NIMA and (d) TATET sites.</figcaption></figure>

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

The thermal demagnetization of a composite three-axes IRM further confirms the presence of a low-coercivity mineral as the main magnetization carrier. In all cases, the curves show the dominance of the magnetically soft fraction (<0.16 T) while the medium and hard components are very minor (Fig. 3). All samples are completely demagnetized at 540 °C, suggesting the presence of magnetite or/and Ti-magnetite as the main magnetic mineral. In the samples from TATET and NIMA sites the presence of a hard magnetic component is also observed. This high-coercivity component is characterized by low Curie temperature of around 300 °C, and could be identified as the High Coercivity, magnetically Stable, Low Temperature (HCSLT) magnetic phase previously recognized in several archaeomagnetic samples as a product of heating of (iron-rich) clays (McIntosh et al. 2007; Hartmann et al. 2011; 2011; Kosterov et al. 2021). This phase could partially contribute to the remanent magnetization together with the dominant magnetite/Ti-magnetite magnetic signal.

## 3.2 Magnetic stability

The magnetic stability of the samples was controlled by both AFD and THD experiments. The results, reported as equal area projections, Zijderveld diagrams and intensity decay plots generally show a stable, single component of magnetization that passes through the origin of the Zijderveld plots (Fig. S1a, Supporting Information).

All samples are completely demagnetized at 540 °C, further confirming the presence of magnetite and/or Ti-magnetite, almost excluding the contribution of hematite as high-coercivity component. Magnetic susceptibility monitoring after each thermal demagnetization step generally shows very small variations, suggesting that most of the samples are thermally stable and suitable for archaeointensity experiments. AFD plots also show very good magnetic stability. Some samples are not completely AF demagnetized at peak field of 180 mT, confirming the presence of the HCSLT magnetic phase, which however generally shows very good magnetic stability. Some samples show a clear secondary magnetic component or disturbed Zijderveld diagrams and were therefore excluded from any further consideration (Fig. S1b, Supporting Information). Moreover, all samples from the SADA2 coffin show significant magnetic susceptibility increase at temperatures around 300–360 °C suggesting insufficient heating in the past and poor thermal stability. The THD and AFD experiments results were used as a pre-selection criterion, rejecting all samples with magnetically unstable behaviour, multi-magnetization components and magnetochemical changes during heating.

## 3.3 Archaeointensity results

A total of 145 specimens from at least four independent ceramic fragments from each studied samples group were measured and
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Thermal demagnetization of orthogonal three-axes IRMs for representative fragments: (a) NIMA-6; (b) TATET-5; (c) TATEZ-2; (d) SADA1-3; (e) SADA4-1 and (f) TEN-2.</figcaption></figure>

analysed: TATET (8 fragments, 29 specimens), TATEZ (9 fragments, 35 specimens), TEN (5 fragments, 14 specimens), NIMA (5 fragments, 24 specimens), SADA1 (7 fragments, 28 specimens) and SADA4 (4 fragments, 15 specimens). Results were elaborated with the Thellier Tool software (Leonhardt *et al.* 2004) and the quality of the archaeointensity determination at specimen level was estimated following rigorous acceptance criteria. Such criteria are similar to those proposed by Hartmann *et al.* (2019) and are based on a quality selection at both specimen and fragment levels.

At specimen level, the intensity value was computed involving the same temperature interval used for isolating the main magnetic component in the orthogonal diagram and at least four temperature steps ($N \ge 4$) [comprising at least 40 per cent of NRM ($f \ge 0.4$)] were used for line slope fitting. Moreover, standard error should be less than 15 per cent ($\beta \le 0.15$), quality index ($q$) above 5 ($q \ge 5$), anchored and unanchored maximum angular deviation (MAD) below 10°, angular difference between anchored and free-floating best-fitting directions below 5° ($\alpha \le 5°$), maximum difference of normalized pTRM checks ($\delta CK$) smaller than 10 per cent ($\delta CK \le 10$ per cent), the normalized ratio of alteration-corrected to non-corrected estimates ($\delta pal$) smaller than 15 per cent ($\delta pal \le 15$ per cent), and maximum difference of normalized pTRM-tail ($\delta TR$) checks below 10 per cent ($\delta TR \le 10$ per cent). At fragment level, a minimum of two successful specimens was necessary for computing an intensity value with an error below 10 per cent after anisotropy correction.

From the 145 studied specimens, 99 (coming from 30 independent fragments) passed the selection criteria reported above, resulting to a success rate of 68 per cent (Fig. 4). Most of the rejected specimens failed the quality selection due to magnetic alteration during the archaeointensity experiments or to unsuccessful magnetic anisotropy determinations. Results at fragment and site level are summarized in Table 2 and analytically presented in Supporting Information (Table S1).

# 4 DISCUSSION

## 4.1 New data and secular variation path in Japan during the last three millennia

The archaeointensity signal is recorded during the last firing of the archaeological baked clay artefacts and thus represents a unique instantaneous picture of the geomagnetic field in the past. When reliable evidence exists to confirm that the different fragments collected from an archaeological site were contemporaneously fired (e.g. baked clay pieces from a kiln), a mean archaeointensity value per site can be calculated, representing a consistent estimation of the geomagnetic field's intensity at the time of their last firing. Nevertheless, it is often difficult (if not impossible) to certify that a set of ceramic fragments found in an archaeological site were fired exactly at the same time and place. For this reason, the calculation of a mean value per site based on independent ceramic fragments should be cautiously performed, as it can hamper the identification of fast geomagnetic field variations (Livermore et al. 2021). In fact, pottery fragments collected from well-dated and short-chronologically constrained archaeological contexts can offer valuable information about rapid changes in the geomagnetic field intensity in a certain geographical area (Shaar et al. 2016; Hervé et al. 2017).

The new high-quality archaeointensity data at fragment level obtained in this study are plotted in Fig. 5(a), together with the mean archaeointensity value calculated for each site. This plot clearly shows that for all samples, the intensity records obtained at fragment level are very well defined. Site mean values for the samples from TEN, NIMA and SADA1 are also well constrained confirming the high-quality of the new data. For site SADA4, the results are divided in two groups, with two fragments showing intensity values of around 40 µT and two fragments showing values of around 50 µT. The samples studied from this site are baked clay pieces from
<figure><img src="image_6.png" /></figure>

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

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

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

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

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

<figure><img src="image_12.png" /><figcaption>**Figure 4.** Representative Arai plots with the corresponding Zijderveld diagrams for (a)–(f) successful samples and (g)–(j) for rejected samples from the different studied sites. In the Arai plots, the empty circles represent NRM remaining against pTRM gained and triangles show pTRM checks performed every two temperature steps.</figcaption></figure>
<table><caption>**Table 2.** New archaeointensity results from Japan. Columns: Frag ID; *n*: number of specimens per fragment; $F_m$ and $\sigma F_m$: mean and error of uncorrected intensity for each fragment (in $\mu$T); $F_{ATRM}$: archaeointensity corrected for the effect of the thermoremanent anisotropy; $F_{m\_ATRM}$ and $\sigma F_{m\_ATRM}$: mean and error of anisotropy corrected intensity for each fragment (in $\mu$T); $F_{m\_ATRM+CR}$ and $\sigma F_{m\_ATRM+CR}$: mean and error of both anisotropy and CR corrected intensity for each fragment (in $\mu$T); $F$ Mean (site) and $\sigma F$ Mean (site): mean archaeointensity value and associated error after all corrections, computed considering at least three fragments per site (in $\mu$T).</caption><thead><tr><th>Frag ID</th><th>n</th><th>F<sub>m</sub></th><th>σF<sub>m</sub></th><th>F<sub>ATRM</sub></th><th>F<sub>m-ATRM</sub></th><th>σF<sub>m-ATRM</sub></th><th>F<sub>m-ATRM+CR</sub></th><th>σF<sub>m-ATRM+CR</sub></th><th>F Mean (site)</th><th>σf Mean (site)</th></tr></thead><tbody><tr><td colspan="11"><em>Tatetsuki (TATET): 160-190 AD [fragments: 5/8; specimens: 14/29]</em></td></tr><tr><td>TATET-2</td><td>2</td><td>50.86</td><td>3.56</td><td>31.90</td><td>32.77</td><td>1.22</td><td>30.64</td><td>1.22</td><td>39.78</td><td>13.03</td></tr><tr><td>TATET-3</td><td>3</td><td>62.09</td><td>0.54</td><td>66.40</td><td>63.10</td><td>4.38</td><td>60.71</td><td>4.38</td><td></td><td></td></tr><tr><td>TATET-6</td><td>5</td><td>37.83</td><td>1.81</td><td>34.07</td><td>32.73</td><td>1.96</td><td>29.46</td><td>1.96</td><td></td><td></td></tr><tr><td>TATET-9</td><td>2</td><td>39.97</td><td>2.82</td><td>37.59</td><td>38.73</td><td>1.61</td><td>34.08</td><td>1.61</td><td></td><td></td></tr><tr><td>TATET-10</td><td>2</td><td>67.15</td><td>0.35</td><td>47.85</td><td>46.33</td><td>2.14</td><td>44.02</td><td>2.14</td><td></td><td></td></tr><tr><td colspan="11"><em>Tatezaka (TATEZ): 160-190 AD [fragments: 5/9; specimens: 19/35]</em></td></tr><tr><td>TATEZ-2</td><td>3</td><td>55.24</td><td>0.14</td><td>46.47</td><td>46.95</td><td>1.02</td><td>42.99</td><td>1.02</td><td>47.99</td><td>9.77</td></tr><tr><td>TATEZ-3</td><td>4</td><td>58.55</td><td>4.17</td><td>34.59</td><td>36.91</td><td>1.59</td><td>35.07</td><td>1.59</td><td></td><td></td></tr><tr><td>TATEZ-7</td><td>4</td><td>66.10</td><td>9.11</td><td>60.95</td><td>58.51</td><td>2.36</td><td>54.97</td><td>2.36</td><td></td><td></td></tr><tr><td>TATEZ-8</td><td>4</td><td>74.87</td><td>0.58</td><td>64.15</td><td>63.64</td><td>1.18</td><td>59.82</td><td>1.18</td><td></td><td></td></tr><tr><td>TATEZ-10</td><td>4</td><td>59.55</td><td>1.30</td><td>44.68</td><td>50.10</td><td>3.77</td><td>47.09</td><td>3.77</td><td></td><td></td></tr><tr><td colspan="11"><em>Tenguyama (TEN): 470-500 AD [fragments: 5/5; specimens: 14/14]</em></td></tr><tr><td>TEN-4</td><td>3</td><td>62.61</td><td>1.50</td><td>54.46</td><td>52.17</td><td>2.49</td><td>49.06</td><td>2.49</td><td>43.45</td><td>3.43</td></tr><tr><td>TEN-6</td><td>3</td><td>51.52</td><td>1.24</td><td>43.75</td><td>42.75</td><td>0.86</td><td>40.33</td><td>0.86</td><td></td><td></td></tr><tr><td>TEN-7</td><td>3</td><td>58.04</td><td>1.48</td><td>47.20</td><td>45.28</td><td>1.69</td><td>42.31</td><td>1.69</td><td></td><td></td></tr><tr><td>TEN-8</td><td>2</td><td>56.51</td><td>2.23</td><td>45.59</td><td>46.61</td><td>1.44</td><td>44.10</td><td>1.44</td><td></td><td></td></tr><tr><td>TEN-9</td><td>3</td><td>50.05</td><td>1.58</td><td>44.68</td><td>44.88</td><td>1.14</td><td>41.42</td><td>1.14</td><td></td><td></td></tr><tr><td colspan="11"><em>Nima Ohtsuka (NIMA): 550-600 AD [fragments: 4/5; specimens: 16/24]</em></td></tr><tr><td>NIMA-2</td><td>3</td><td>48.20</td><td>2.10</td><td>51.83</td><td>47.78</td><td>3.79</td><td>44.31</td><td>3.79</td><td>44.40</td><td>1.60</td></tr><tr><td>NIMA-4</td><td>5</td><td>47.20</td><td>4.30</td><td>46.90</td><td>45.60</td><td>2.30</td><td>46.20</td><td>2.30</td><td></td><td></td></tr><tr><td>NIMA-7</td><td>3</td><td>55.48</td><td>2.01</td><td>42.91</td><td>43.24</td><td>0.45</td><td>42.29</td><td>0.45</td><td></td><td></td></tr><tr><td>NIMA-8</td><td>5</td><td>53.10</td><td>3.00</td><td>51.53</td><td>49.90</td><td>2.40</td><td>44.80</td><td>2.40</td><td></td><td></td></tr><tr><td colspan="11"><em>Sada Higashizuka (SADA1): 630 AD [fragments: 7/7; specimens: 26/28]</em></td></tr><tr><td>SADA1-1</td><td>4</td><td>39.42</td><td>3.89</td><td>37.89</td><td>36.33</td><td>1.80</td><td>33.78</td><td>1.80</td><td>31.97</td><td>2.21</td></tr><tr><td>SADA1-2</td><td>4</td><td>36.19</td><td>0.65</td><td>30.01</td><td>32.58</td><td>2.37</td><td>29.98</td><td>2.37</td><td></td><td></td></tr><tr><td>SADA1-3</td><td>3</td><td>37.76</td><td>3.37</td><td>33.96</td><td>34.20</td><td>1.43</td><td>32.14</td><td>1.43</td><td></td><td></td></tr><tr><td>SADA1-4</td><td>4</td><td>38.18</td><td>2.36</td><td>33.19</td><td>31.30</td><td>1.82</td><td>29.42</td><td>1.82</td><td></td><td></td></tr><tr><td>SADA1-5</td><td>4</td><td>41.60</td><td>3.52</td><td>38.89</td><td>35.52</td><td>2.61</td><td>33.39</td><td>2.61</td><td></td><td></td></tr><tr><td>SADA1-6</td><td>3</td><td>46.99</td><td>2.54</td><td>39.25</td><td>39.41</td><td>1.98</td><td>35.07</td><td>1.98</td><td></td><td></td></tr><tr><td>SADA1-9</td><td>4</td><td>35.32</td><td>3.84</td><td>34.81</td><td>31.23</td><td>3.05</td><td>29.98</td><td>3.05</td><td></td><td></td></tr><tr><td colspan="11"><em>Sada Nishizuka—Coffin no. 4 (SADA4): 650-675 AD [fragments: 4/4; specimens: 10/15]</em></td></tr><tr><td>SADA4-1</td><td>3</td><td>45.67</td><td>1.77</td><td>37.26</td><td>41.50</td><td>3.94</td><td>37.77</td><td>3.94</td><td>45.60</td><td>7.83</td></tr><tr><td>SADA4-2</td><td>2</td><td>59.44</td><td>4.37</td><td>63.16</td><td>60.60</td><td>3.62</td><td>53.33</td><td>3.62</td><td></td><td></td></tr><tr><td>SADA4-4</td><td>3</td><td>53.71</td><td>4.55</td><td>58.17</td><td>53.96</td><td>3.70</td><td>51.26</td><td>3.70</td><td></td><td></td></tr><tr><td>SADA4-7</td><td>2</td><td>48.46</td><td>0.98</td><td>45.38</td><td>46.04</td><td>0.93</td><td>40.05</td><td>0.93</td><td></td><td></td></tr></tbody></table>

the same ceramic coffin and it would be thus expected to give similar archaeointensity values, representative of the Earth's magnetic field intensity at the time of the coffin's production. A difference of ~10 $\mu$T is unexpectedly high (taking into consideration the very rigorous experimental protocol applied in this study), but it is however comparable with the often highly scattered archaeointensity values at site level reported in the literature. Data from SADA1 coffin are well concentrated around the mean value and suggest the presence of an intensity low at the beginning of the 7th century AD. However, what is particularly interesting in the new results is the large variability of the intensity values observed for the contemporaneous sites of TATET and TATEZ (Fig. 5a). Fragments from these two sites show very similar behaviour, with some of them being characterized by high intensity values of around 60 $\mu$T and others by low values of around 30–35 $\mu$T. These sites are both very well dated at 160–190 AD, representing a well constrained and short-time period. Careful typological examination of the fragments collected from both sites excludes possible contamination from different stratigraphic levels, confirming that all samples were produced in the same period

(160–190 AD). The high quality of the new intensity data, accompanied by the very reliable and precise dating of these two sites may offer, for the first time, evidence of very fast geomagnetic field variations in Japan during the Late Yayoi period, even though more high-quality contemporaneous data are still necessary to better constrain such event.

Comparison of the new data with literature data for Japan (GEOMAGIA50.v3.4 database, Brown et al. 2021) shows that the new data are characterized by systematically lower intensity values with respect to the previously published data (Fig. 5b). Such important difference could be explained by the almost complete lack of CR and anisotropy corrections in the literature data. Nowadays, it is well known that the fast-cooling time used during the archaeointensity laboratory experiments overestimates the intensity results and thus the correction of the intensity data for this effect is very important (Genevey et al. 2008). Moreover, archaeological artefacts such as pottery, bricks and tiles may be highly anisotropic (Tema 2009), significantly influencing the archaeointensity determination that needs to be corrected (Chauvin et al. 2000). As previously
<figure><img src="image_13.png" /><figcaption>**Figure 5.** (a) New archaeointensity results at fragment level (dot) plotted together with the mean archaeointensity value calculated for each sample collection (star); (b) new data plotted together with the intensity values from Kitahara *et al.* (2021), literature data from Japan (GEOMAGIA50.v3.4 database) and data from South Korea (Hong *et al.* 2013), as well with the most recent global geomagnetic field models. All data are relocated to the latitude of Tokyo (35.65°N).</figcaption></figure>

discussed, most of the literature data from Japan mainly come from
studies carried out more than 20-30 yr ago, when such corrections
were not so well established yet, most probably resulting to higher-
than-expected intensity values. Re-evaluating such old data is an
extremely difficult, if not impossible task, as in most cases the orig-
inal data are not available, and it is not possible to have access
to the raw measurements for evaluating their quality. Even if we
could apply an arbitrary CR correction of 5 per cent decrease for
the data which were not originally corrected for the CR effect (fol-
lowing Genevey et al. 2008), it would be impossible to correct them

a posteriori for the anisotropy effect or to estimate their thermal
stability (and thus reliability) without the performance of pTRM
checks.

Nevertheless, comparison of our new data with the results ob-
tained from the most recent archaeointensity study in Japan pub-
lished by Kitahara *et al.* (2021) shows very good agreement
(Fig. 5b). Kitahara *et al.* (2021) studied archaeological baked clays
from the Suemura old kiln complex (located in the Senboku hills,
Osaka Prefecture), dated from 400 to 700 AD, by applying the
Tsunakawa–Shaw archaeointensity method and they also found
much lower intensity values in respect to the previously published Japanese data. Unfortunately, no contemporaneous data from nearby to Japan countries are available for comparison with our new data. However, data from South Korea (Hong *et al.* 2013), if compared with the literature data for Japan, generally show lower intensity values, too (Fig. 5b). The data from South Korea are corrected for the anisotropy effect through the investigation of the anisotropy of the anhysteretic remanent magnetization, but they still miss CR correction (Hong *et al.* 2013).

Comparison with geomagnetic field models shows that the new intensities are lower than the model's predictions for the Late Yayoi and Kofun periods (Fig. 5b). Among the models used for comparison, ARCH3k.1 (Korte *et al.* 2009), pfm9k.1a (Nilsson *et al.* 2014), HFM.OL.A1 (Constable *et al.* 2016) and SHAWQ2k (Campuzano *et al.* 2019), clearly follow a much higher intensity trend in respect to the Kitahara's *et al.* (2021) and our new data values. This is expected for ARCH3k.1 and SHAWQ2k models as they are only based on data from archaeological artefacts and volcanic rocks, being strongly influenced by the higher intensity values of the older reference data for Japan. The new data are in better agreement with the lower values given by the CALS3k.4 model (Korte & Constable 2011) which shows low intensities, being based on sedimentary records too. Such comparison underlines the importance of the new high-quality data presented in this study for the improvement of the geomagnetic field models' predictions in the area.

## 4.2 Sharp intensity variations in contemporaneous fragments

The new absolute intensity data from TATET and TATEZ sites show significant differences in intensities which could be interpreted as fast geomagnetic field variations when values at fragment level are considered. In fact, these two data sets show differences in intensities among contemporaneous fragments that vary from 29.46 ± 1.96 up to 60.71 ± 4.38 μT over an interval of around 30 yr, as constrained by chronological and archaeological evidence. In both sites, intensity rates suggest faster variations than those observed in the current field (Alken et al. 2021) and the upper bound value from purely advection scenarios of the field generated within the Earth's outer core (Livermore et al. 2014). The fast intensity changes depicted from the new data can reach values comparable to extreme intensity field variations, the so-called geomagnetic spikes (Ben-Yosef et al. 2009; Shaar et al. 2011; Livermore et al. 2014; Shaar et al. 2016; Ben-Yosef et al. 2017), even though with much milder maximum values than those seen in the Levantine Iron Age Anomaly (LIAA). However, when the site mean intensities of these two sites are considered, the intensity variation rates are much smoother pointing that only a strict fragment level interpretation of our data sustains very sharp intensity variations. Indeed, Livermore et al. (2021) showed that data sets interpreted as geomagnetic spikes can have much lower rates of intensity variation when data are considered only at the level of a group of fragments from the same archaeological site. They also showed that defining and detecting fast geomagnetic field variations strongly depends on the approach imposed to the selection and interpretation of the reference data (both in terms of experimental robustness and well-constrained dating). Even though our new data from TATET and TATEZ archaeological sites encourage the possibility of likewise Levantine sharp intensity magnetic field variations in Japan during the Late Yayoi period, which might be crucial for further considerations on the spatial and temporal occurrence of such geomagnetic

field phenomena, inconsistency in the interpretation of the results at fragment level cannot be excluded. Undoubtedly, more data from this period are still needed to assert the puzzling contemporaneous intensity values found in the TATET and TATEZ sites and to better constrain the evidence of such sharp geomagnetic field variations in Japan.

## 4.3 Evidence for recurrence of the West Pacific Anomaly

The new absolute intensity data from TEN, NIMA and SADA1 sites are very well constrained at both fragment and site mean levels, showing a clear drop in intensity from the late 5th to early 7th century AD, with a well-defined intensity low at 630 AD. Such weak intensity values, clearly lower than those depicted in the previous archaeomagnetic data available for the Late Kofun period in Japan (Brown et al. 2021), are further supported by the recent data from Kitahara et al. (2021), which also show weak intensity values for the same period, including the intensity drop at the beginning of the 7th century AD, even though less sharp (Fig. 5). The minimum intensity value obtained from the SADA1 site is lower than the intensity minima observed in South Korea at 700 AD, and 1600 AD (Hong et al. 2013, see Fig. 5b). To further investigate the occurrence of such weak geomagnetic field intensity field in Japan, we have used the SHAWQ2k global model (Campuzano et al. 2019), to create snapshots of the surface intensity field at 500, 600, 700, 1500, 1600 and 1700 AD, identifying local minima centres as extreme low intensity values (Fig. 6). All snapshots clearly show the global intensity minimum located at Africa, while two other minima are observed at the west of South America and at the east of East Asia. The East Asia minimum is not evident at the 500 and 600 AD snapshots while it only weakly appears at 700 AD (Fig. 6, left-hand side). On the contrary, such intensity minimum is much more evident in the 1500–1700 AD period (Fig. 6, right-hand side). This is in contrast with our new data which depict lower values in Japan during the 500–700 AD period than those seen in South Korea for the 1500–1700 AD (Hong et al. 2013). The 1600 AD minimum has been attributed to the influence of a weak surface field anomaly centered in the equatorial position at East Asia, the so-called West Pacific Anomaly (WPA, e.g. He et al. 2021). He et al. (2021) showed that this anomaly prompts equatorial aurorae at South Korea between 1590 and 1720 AD. Furthermore, archaeomagnetic results in Cambodia show an intensity dip between 1100 and 1300 AD, accompanied by fast directional change between 1200 and 1300 AD (Cai et al. 2021). Our new data from Middle to Late Kofun period, give evidence of an ancient occurrence of the WPA influence over East Asia, being also supported by the recent archaeomagnetic records from Japan and China (Cai et al. 2021; Kitahara et al. 2021). All these recent data are not included in the SHAWQ2k model's calculations, preventing it from depicting the possible recurrence of weak ancient surface geomagnetic field in East Asia suggested here.

The persistent location of recurrence of surface intensity minima (Terra-Nova et al. 2019; Nilsson et al. 2022), which might include the West Pacific as supported by our data, is related to core–mantle dynamics. Archaeomagnetic field models reveal that lateral mantle heterogeneities can localize core–mantle boundary (CMB) reversed flux patches (Terra-Nova et al. 2016). The expansion of reversed flux patches below South Atlantic during the historical period prompts a large region of weak surface field (e.g Jackson et al. 2000), known as South Atlantic anomaly (SAA), for which mantle driven anomaly has been proposed (Tarduno et al. 2015; Terra-Nova et al. 2019).
<figure><img src="image_14.png" /><figcaption>**Figure 6.** Surface intensity field as predicted by the SHAWQ2k global geomagnetic field model (Campuzano *et al.* 2019) at (a) 500 AD, (b) 600 AD, (c) 700 AD, (d) 1500 AD, (e) 1600 AD and (f) 1700 AD. White diamonds indicate the local minima centres.</figcaption></figure>

Based on archaeomagnetic data from Africa which reveal sharp directional changes and an intensity drop, Tarduno *et al.* (2015) proposed that the radial length scale of the field might be reduced by the presence of the Africa large low-shear-velocity province (LLSVP), thus prompting the emergence of reversed flux patches. Using numerical dynamo simulations that account for lateral lower mantle heterogeneities, Terra-Nova *et al.* (2019) showed that the recurrent position of surface minima can be associated to localized fluid upwelling at the top of the core maintained by regions where less heat extraction occurs across the CMB, at LLSVPs. Using core flow models, He *et al.* (2021) showed that the wax and wane on intensity requires localized horizontally diverging structures, suggesting that upwelling flows are key to intensity drops and downwelling to the restrengthening of local field. Similar results were obtained by Trindade *et al.* (2018). They used synthetic CMB magnetic field kinematic scenarios to show that a reversed flux patch, expanding and drifting west and southward below the measure point, could represent an optimal scenario for sharp local intensity changes. Additionally, the centres of the ancient WPA and the present-day SAA are located relatively close, at the edges of the Pacific and Africa LLSVPs, respectively (Terra-Nova *et al.* 2019; He *et al.* 2021). The WPA might, thus, be a mantle-driven anomaly similar to the SAA.

Fig. 7 shows the time dependence of field intensity minima in longitude identified for a set of archaeomagnetic field models (Korte *et al.* 2009; Nilsson *et al.* 2014; Constable *et al.* 2016; Campuzano *et al.* 2019), and the longitudinal stack (the sum of all values for each latitude, illustrated with the blue line in Fig. 7) of seismic shear wave

velocity anomalies from the tomography model of Masters *et al.* (2000) at the lower mantle, truncated at spherical harmonic degree 6. We assume that the seismic shear wave velocity anomalies reflect the lateral heterogeneities of the heat extraction from the outer core by the mantle (e.g. Gubbins *et al.* 2007). At locations where there is more heat extraction, the fluid is colder than sinks, while below regions of less heat extraction fluid upwelling occurs (e.g. Olson & Christensen 2002). Upwelling plumes carry toroidal field lines upwards accumulating magnetic flux at the top of the core which is then expelled by magnetic diffusion and creates pairs of normal and reversed flux patches on the CMB (Bloxham 1986). Several minima located at the edges of the LLSVP can be identified (Fig. 7), as well as at its centre (e.g. around longitude 0°). It is worth noting that surface intensity minima between 2000–6000 BC tend to be less scatter, being recurrent around the position of the SAA and WPA, thus possibly reflecting the location of upwelling plumes at the top of the core. The drifting patterns of these anomalies depend on the accuracy of their identification (Amit *et al.* 2021). Here, surface minima drift both west and eastward. Edges of LLSVP might act as anchors where surface minima anomalies present slower drifting velocities and therefore are seen for longer time intervals. This could explain why the SAA present-day westward drift is getting slower (Terra-Nova *et al.* 2017) as it passes through the edges of the African LLSVP. Evidence for recurrent weak surface field anomalies in the Holocene due to the emergence of reversed flux patches, similar to the present-day SAA are supported not only by the intensity field drop but also by sharp changes in the field's direction as observed by archaeomagnetic studies in South America (Trindade *et al.* 2018),
<figure><img src="image_15.png" /><figcaption>**Figure 7.** Identification of intensity minima longitudes for up to the last 10 000 yr for some of the most used global geomagnetic field models. Blue line denotes the non-dimensional sectorial average of the shear seismic wave anomaly for the lower mantle calculated based on the tomographic model of Masters *et al.* (2000), truncated at spherical harmonic degree 6. Dashed lines indicate the position of surface weak intensity field anomalies: red for the present-day SAA (Terra-Nova *et al.* 2017) and black for an ancient position of the WPA (He *et al.* 2021).</figcaption></figure>

Africa (Tarduno *et al.* 2015; Hare *et al.* 2018) and West Pacific (Cai *et al.* 2021).

All archaeomagnetic field models fail to reproduce the intensity drop seen in our data around the beginning of the 7th century AD (Fig. 5b), most probably due to the models’ spatial and temporal resolution limitation caused both by the poor data coverage as well as by the often low archaeomagnetic data quality (Brown *et al.* 2021). Actually, only the CALS3k.4 model (Korte & Constable 2011) shows intensity values close or lower than 30 µT, similar to our new data. The most recently published models show clearly stronger intensity field in Japan during the 500–700 AD period, being clearly influenced by the older Japanese archaeointensity data. Fig. 8 shows snapshots of the radial field at CMB below Japan at 500, 600 and 700 AD as calculated from the archaeomagnetic field models CALS3k.4 (Korte & Constable 2011), HFM.10k.1A (Constable *et al.* 2016) and SHAWQ2k (Campuzano *et al.* 2019). Intense normal flux patches are located below Japan for both HFM.10k.1A and SHAWQ2k models. In the HFM.10k.1A model, an extensive patch dominates the field while in the SHAWQ2k model patches drift, merge and concentrate below Japan (Fig. 8). Oppositely, the CALS3k.4 model resolves weak normal radial field below Japan. These differences between the predictions of the radial field below Japan based on the different models used, clearly evidence that more high-quality data are still needed and improvements on the treatment of the data’s uncertainties are still necessary to reliably constrain the radial field below Japan and thus investigate the source of the intensity decrease observed in our data. We hope that the incorporation of our new data in the calculations of the future archaeomagnetic field models could contribute to improve our current knowledge of the geomagnetic field variation in East Asia during the first millennium AD and help in the prediction of the WPA temporal evolution.

## 5 CONCLUSIONS

This study reports 30 new archaeointensity results from six archaeological sites located at the Okayama Prefecture, Japan, belonging to the Late Yayoi and Kofun periods. The new data, which are corrected for both anisotropy and CR effects and are carefully selected, clearly show that most of the previous archaeointensity determinations from Japan are overestimated, most probably due to the lack of necessary corrections. The new data show a decrease in the geomagnetic field intensity from 5th to 7th centuries AD, with a well constraint intensity low around 630 AD. Significant differences in intensity values obtained from well chronologically constrained, contemporaneous ceramic fragments from the Late Yayoi period (around 200 AD) may give evidence of high intensity variation rates, reaching values similar to those calculated during the LIAA, but with much lower absolute intensity values. The new high-quality data, together with data from South Korea and China, support a possible ancient recurrence of the WPA. Our assessment using archaeomagnetic field models and topographic shear wave seismic anomalies at the lowermost mantle supports mantle driven WPA. The new data offer a significant contribution to improve the global geomagnetic field models predictions in Japan, even if more data are still necessary to better understand the Earth’s magnetic field evolution in East Asia.

## ACKNOWLEDGMENTS

This study has been partially financed by the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie RISE Action ‘BEYOND ARCHAEOLOGY: An advanced approach linking East to West through Science, Field
<figure><img src="image_16.png" /><figcaption>**Figure 8.** Radial field at the CMB calculated from the archaeomagnetic field models (a)–(c) SHAWQ2k (Campuzano *et al.* 2019); (d)–(f) HFM.OL1.A1 (Constable *et al.* 2016) and (g)–(i) CALS3k.4 (Korte & Constable 2011) at 500, 600 and 700 AD. White diamonds indicate the centre of intense patches.</figcaption></figure>

archaeology and Interactive Museum Experiences', project grant agreement no. 823826. The European Commission's support for the production of this publication does not constitute an endorsement of the contents, which reflect the views only of the authors, and the Commission cannot be held responsible for any use which may be made of the information contained therein. GH acknowledges the National Scientific and Technological Development Council-CNPq (grants nos 425728/2018-8 and 312737/2020-3). TH is supported by JSPS Kakenhi grant (no. 20H00028). FTN was supported by the Centre national d'etudes spatiales (CNES). The authors acknowledge the editor Andrew Biggin for the attention given to the manuscript. Ron Shaar and an anonymous reviewer are sincerely acknowledged for their constructive reviews.

## DATA AVAILABILITY

The data produced in this study can be downloaded from https://earthref.org/ERDA/2547 and are available under request.

## CONFLICT OF INTEREST

The authors declare that they have no conflict of interests.

## REFERENCES

* Alken, P. *et al.* 2021. International Geomagnetic Reference Field: the thirteenth generation, *Earth Planets Space*, **73**, 49, doi:10.1186/s40623-020-01288-x.

* Amit, H., Terra-Nova, F., Lezin, M. & Trindade, R. I., 2021. Non-monotonic growth and motion of the South Atlantic Anomaly, *Earth Planets Space*, **73**, doi: 10.1186/s40623-021-01356-w.

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

* Ben-Yosef, E., Millman, M., Shaar, R., Tauxe, L. & Lipschits, O., 2017. Six centuries of geomagnetic intensity variations recorded by royal Judean stamped jar handles, *Proc. Natl. Acad. Sci. U.S.A.*, **114**, 2160–2165.

* Bloxham, J., 1986. The expulsion of magnetic flux from the Earth's core, *Geophys. J. R. astr. Soc.*, **87**, 669–678.
Brown, M. C., Hervé, G., Korte, M. & Genevey, A., 2021. Global archaeomagnetic data: The state of the art and future challenges, *Phys. Earth planet. Inter.*, **318**, 106766, doi:10.1016/j.pepi.2021.106766.

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

Cai, S., Doctor, R., Tauxe, L., Hendrickson, M., Hua, Q., Leroy, S. & Phon, K., 2021. Archaeomagnetic results from Cambodia in Southeast Asia: evidence for possible low-latitude flux expulsion, *Proc. Nat. Acad. of Sci.*, **118**(11), e2022490118, doi:10.1073/pnas.2022490118.

Campuzano, S.A., Gómez-Paccard, M., Pavón-Carrasco, F.J. & Osete, M.L., 2019. Emergence and evolution of the South Atlantic Anomaly revealed by the new paleo-magnetic reconstruction SHAWQ2k, *Earth planet. Sci. Lett.*, **512**, 17–26.

Chauvin, A., Garcia, Y., Lanos, P. & Laubenheimer, F., 2000. Paleointensity of the geomagnetic field recovered on archaeomagnetic sites from France, *Phys. Earth planet. Inter.*, **120**, 111–136.

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

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

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

Davies, C.J. & Constable, C.G., 2018. Searching for geomagnetic spikes in numerical dynamo simulations, *Earth planet. Sci. Lett.*, **504**, 72–83.

Davies, C.J. & Constable, C.G., 2020. Rapid geomagnetic changes inferred from Earth observations and numerical simulations, *Nat. Commun.*, **11**, 3371, doi:10.1038/s41467-020-16888-0.

Di Chiara, A., Herrero-Bervera, E. & Tema, E. 2020. Geomagnetic field variations in the past: an introduction, *Geol. Soc., Lond., Spec. Publ.*, **497**(1), 1, doi:10.1144/SP497-2020-78.

Genevey, A., Gallet, Y., Constable, C.G., Korte, M. & Hulot, G., 2008. ArcheoInt: an upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment, *Geochem. Geophys. Geosyst.*, **9**(4), doi:10.1029/2007GC001881.

Gubbins, D., Willis, A.P. & Sreenivasan, B., 2007. Correlation of Earth's magnetic field with lower mantle thermal and seismic structure, *Phys. Earth planet. Inter.*, **162**, 256–260.

Hare, V.J., Tarduno, J.A., Huffman, T., Watkeys, M., Thebe, P.C., Manyanga, M., Bono, R. K. & Cottrell, R. D., 2018. New archeomagnetic directional records from Iron Age southern Africa (ca. 425–1550 CE) and implications for the South Atlantic Anomaly, *Geophys. Res. Lett.*, **45**, 1361–1369.

Hartmann, G., Poletti, W., Trindade, R., Ferreira, L. & Sanches, P., 2019. New archeointensity data from South Brazil and the influence of the South Atlantic Anomaly in South America, *Earth planet. Sci. Lett.*, **512**, 124–133.

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

He, F. *et al.* 2021. Equatorial auroral records reveal dynamics of the paleo-West Pacific geomagnetic anomaly, *Proc. Nat. Acad. of Sci.*, **118**(20), e2026080118, doi:10.1073/pnas.2026080118.

Hervé, G. *et al.* 2017. Fast geomagnetic field intensity variations between 1400 and 400 BCE: new archaeointensity data from Germany, *Phys. Earth planet. Inter.*, **270**, 143–156.

Hirooka, K., 1971. Archaeomagnetic study for the past 2,000 years in Southwest Japan. *Memoirs of the Faculty of Science, Kyoto University, Ser. Geol. Mineral.*, **38**(2), 167–207.

Hirooka, K., 1983. Results from Japan, in "Geomagnetism of Baked Clays and Recent Sediments", Creer, K.M., Tucholka, P. & Barton, C.E.(Eds.), Elsevier, Amsterdam, pp. 150–157.

Hong, H., Yu, Y., Lee, C.H., Kim, R.H., Park, J., Doh, S.J., Kim, W. & Sung, H., 2013. Globally strong geomagnetic field intensity circa 3000 years ago, *Earth planet. Sci. Lett.*, **383**, 142–152.

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

Kitahara, Y., Nishiyama, D., Ohno, M., Yamamoto, Y., Kuwahara, Y. & Hatakeyama, T., 2021. Construction of new archaeointensity reference curve for East Asia from 200 CE to 1100 CE, *Phys. Earth planet. Inter.*, **310**, 106596, doi:10.1016/j.pepi.2020.106596.

Kitahara, Y., Yamamoto, Y., Ohno, M., Kuwahara, Y., Kameda, S. & Hatakeyama, T., 2018. Archeointensity estimates of a tenth-century kiln: first application of the Tsunakawa–Shaw paleointensity method to archeological relics, *Earth Planets Space.*, **70**, doi: 10.1186/s40623-018-0841-5.

Kitazawa, K., 1970. Intensity of the geomagnetic field in Japan for the past 10,000 years, *J. geophys. Res.*, **75**, 7403–7411.

Korte, M. & Constable, C.G., 2018. Archeomagnetic intensity spikes: global or regional geomagnetic field features?, *Front. Earth Sci.*, **6**, 17.

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

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

Kosterov, A. *et al.* 2021. High-coercivity magnetic minerals in archaeological baked clay and bricks, *Geophys. J. Int.*, **224**, 1256–1271.

Leonhardt, R., Heunemann, C. & Krasa, D., 2004. Analyzing absolute paleointensity determinations: acceptance criteria and the software Thellier-Tool4.0., *Geochem., Geophys., Geosyst.*, **5**, doi:10.1029/2004GC000807.

Livermore, P., Gallet, Y. & Fournier, A., 2021. Archeomagnetic intensity variations during the era of geomagnetic spikes in the Levant, *Phys. Earth planet. Inter.*, **312**, 106657. doi:10.1016/j.pepi.2021.106657.

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

Lowrie, W., 1990. Identification of ferromagnetic minerals in a rock by coercivity and unblocking temperature properties, *Geophys. Res. Lett.*, **17**, 159–162.

Masters, G., Laske, G., Bolton, H. & Dziewonski, A., 2000. The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle: implications for chemical and thermal structure, in: Karato S, Forte A, Liebermann R, Masters G, Stixrude L, eds. Earth's Deep Interior: Mineral Physics and Tomography from the Atomic to the Global Scale. AGU Monography, 117, pp. 63–87.

McIntosh, G., Kovacheva, M., Catanzariti, G., Donadini, F. & Osete, M.L., 2011. High coercivity remanence in baked clay materials used in archaeomagnetism, *Geochem. Geophys. Geosyst.*, **12**(2), Q02003, doi:10.1029/2010GC003310.

McIntosh, G., Kovacheva, M., Catanzariti, G., Osete, M.L. & Casas, L., 2007. Widespread occurrence of a novel high coercivity, thermally stable, low unblocking temperature magnetic phase in heated archeological material, *Geophys. Res. Lett.*, **34**, L21302, doi:10.1029/2007GL031168.

Mochizuki, N., Tsunakawa, H., Oishi, Y., Wakai, S., Wakabayashi, K.I. & Yamamoto, Y., 2004. Palaeointensity study of the Oshima 1986 lava in Japan: implications for the reliability of the Thellier and LTD-DHT Shaw methods, *Phys. Earth planet. Inter.*, **146**, 395–416.

Nagata, T., Arai, Y. & Momose, K., 1963. Secular variation of the geomagnetic total force during the last 5000 years, *J. geophys. Res.*, **68**, 5277–5281.

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

Nilsson, A., Suttie, N., Stoner, J. S. & Muscheler, R., 2022. Recurrent ancient geomagnetic field anomalies shed light on future evolution of the South Atlantic Anomaly, *Proc. Natl. Acad. Sci.*, **119**(24), e2200749119, doi:10.1073/pnas.2200749119.

Olson, P. & Christensen, U. R., 2002. The time-averaged magnetic field in numerical dynamos with non-uniform boundary heat flow, *Geophys. J. Int.*, **151**, 809–823.
Osete, M.L. et al. 2020. Two archaeomagnetic intensity maxima and rapid directional variation rates during the Early Iron Age observed at Iberian coordinates. Implications on the evolution of the Levantine Iron Age Anomaly, *Earth planet. Sci. Lett.*, **533**, 116047, doi:10.1016/j.epsl.2019.116047.

Pavón-Carrasco, F. J., Campuzano, S.A., Rivero-Montero, M., Molina-Cardín, A., Gómez-Paccard, M. & Osete, M. L., 2021. SCHA.DIF.4k: 4,000 years of paleomagnetic reconstruction for Europe and its application for dating, *J. geophys. Res.: Solid Earth*, **126**, e2020JB021237, doi:10.1029/2020JB021237.

Pavo'n-Carrasco, F. J., Go'mez-Paccard, M., Herve', G., Osete, M. L. & Chauvin, A., 2014. Intensity of the geomagnetic field in Europe for the last 3 ka: influence of data quality on geomagnetic field modeling, *Geochem. Geophys. Geosyst.*, **15**, 2515–2530.

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

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

Sakai, H. & Hirooka, K., 1986. Archaeointensity determinations from Western Japan, *J. Geomagn. Geoelectr.*, **38**, 1323–1329.

Sasajima, S., 1965. Geomagnetic secular variation revealed in the baked earths in West Japan (part 2) change of the field intensity, *J. Geomagn. Geoelectr.*, **17**, 413–416.

Shaar, R., Ben-Yosef, E., Ron, H., Tauxe, L., Agnon, A. & Kessel, R., 2011. Geomagnetic field intensity: how high can it get? How fast can it change? Constraints from iron-age copper-slag, *Earth planet. Sci. Lett.*, **301**, 297–306.

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

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

Tema, 2009. Estimate of the magnetic anisotropy effect on the archaeomagnetic inclination of ancient bricks, *Phys. Earth planet. Inter.*, **176**, 213–223, doi:10.1016/j.pepi.2009.05.007.

Terra-Nova, F., Amit, H. & Choblet, G., 2019. Preferred locations of weak surface field in numerical dynamos with heterogeneous core-mantle boundary heat flux: consequences for the South Atlantic Anomaly, *Geophys. J. Int.*, **217**, 1179–1199.

Terra-Nova, F., Amit, H., Hartmann, G.A., Trindade, R.I.F. & Pinheiro, K.J., 2017. Relating the South Atlantic Anomaly and geomagnetic flux patches, *Phys. Earth planet. Inter.*, **266**, 39–53.

Terra-Nova, F., Amit, H., Hartmann, G. A. & Trindade, R. I. F., 2016. Using archaeomagnetic field models to constrain the physics of the core:

robustness and preferred locations of reversed flux patches, *Geophys. J. Int.*, **206**, 1890–1913.

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

Trindade, R.I.F. et al. 2018. Speleothem record of geomagnetic South Atlantic Anomaly recurrence, *Proc. Natl. Acad. Sci.*, **115**(52), 13198-13203.

Veitch, R., Hedley, I. & Wagner, J.J., 1984. An investigation of the intensity of the geomagnetic field during Roman times using magnetically anisotropic bricks and tiles, *Arch. Sci.*, **37**, 359–373.

Watanabe, N., 1958. Secular variation in the direction of geomagnetism as the standard scale for geomagnetochronology in Japan, *Nature*, **182**, 383–384

Watanabe, N., 1959. The direction of remanent magnetism of baked earth and its application to chronology for anthropology and archaeology in Japan; an introduction to geomagnetochronology, *J. Fac. Sci. Univ. Tokyo*, **2**, 1–188.

Yamamoto, Y. & Hoshi, H., 2008. Paleomagnetic and rock magnetic studies of the Sakurajima 1914 and 1946 andesitic lavas from Japan: a comparison of the LTD-DHT Shaw and Thellier paleointensity methods, *Phys. Earth planet. Inter.*, **167**, 118–143.

Yamamoto, Y., Torii, M. & Natsuhara, N., 2015. Archeointensity study on baked clay samples taken from the reconstructed ancient kiln: implication for validity of the Tsunakawa-Shaw paleointensity method, *Earth Planet Space*, **67**, 63, doi:10.1186/s40623-015-0229-8.

Yoshihara, A, Kondo, A, Ohno, M & Hamano, Y., 2003. Secular variation of the geomagnetic field intensity during the past 2000 years in Japan, *Earth planet. Sci. Lett.*, **210**, 219–231.

Yu, Y., 2012. High-fidelity paleointensity determination from historic volcanoes in Japan, *J. geophys. Res.*, **117**, doi:10.1029/2012JB009368.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

Fig. S1. Representative thermal demagnetization results. (a) Representative samples that showed stable magnetic behaviour and were selected for archaeointensity experiments; (b) samples characterized by important magnetic susceptibility enhancement, secondary magnetic components and/or nonlinear Zijderveld diagrams. These samples were excluded from further archaeointensity investigations. Left: intensity decay curves; middle: Zijderveld diagrams and right: monitoring of bulk magnetic susceptibility at room temperature after each heating/cooling step.

Table S1. Archeointensity results.

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