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 and Planetary Science Letters 533 (2020) 116047

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

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 M.L. Osete a,b,∗ , A. Molina-Cardín a,b , S.A. Campuzano c , G. Aguilella-Arzo d , A. Barrachina-Ibañez d , F. Falomir-Granell d , A. Oliver Foix d , M. Gómez-Paccard b , F. Martín-Hernández a,b , A. Palencia-Ortas a , F.J. Pavón-Carrasco a,b , M. Rivero-Montero b a

Departamento de Física de la Tierra, Astronomía y Astrofísica I (Geofísica y Meteorología), Facultad de Ciencias Físicas, Universidad Complutense de Madrid, Plaza de las Ciencias 1, 28040 Madrid, Spain b Instituto de Geociencias (UCM-CSIC), Facultad de Ciencias Físicas, Universidad Complutense de Madrid, Avda Complutense s/n, 28040 Madrid, Spain c Istituto Nazionale di Geofisica e Vulcanologia, Sezione Roma 2, Via di Vigna Murata 605, 00143 Roma, Italy d Servicio de Investigaciones Arqueológicas y Prehistórica de la Diputación de Castellón, Museu Belles Arts, Av. Germans Bou, 28, 12003 Castelló, Spain

a r t i c l e

i n f o

Article history: Received 2 August 2019 Received in revised form 16 December 2019 Accepted 18 December 2019 Available online 14 January 2020 Editor: M. Ishii Keywords: archaeomagnetism archaeointensity secular variation Iberian Peninsula Levantine Iron Age Anomaly (LIAA) geomagnetic modelling

a b s t r a c t Variations of geomagnetic field in the Iberian Peninsula prior to Late Iron Age times are poorly constrained. Here we report 14 directional and 10 palaeointensity results from an archaeomagnetic study carried out on 17 combustion structures recovered from six archaeological sites in eastern Spain. The studied materials have been dated by archaeological evidences and supported by radiocarbon dates (8th-5th centuries BC). Rock magnetic experiments indicate that the characteristic remanent magnetization is carried by a low coercivity magnetic phase with Curie temperatures of 500-575 ◦ C, most likely titanomagnetite/maghemite with low titanium content. Archaeointensity determinations were carried out by using the classical Thellier-Thellier experiment including pTRM-checks and magnetic anisotropy corrections. A new full vector Iberian Paleosecular Variation Curve for the Iron Age is presented. High fluctuation rates on both directions and intensities are observed during the Early Iron times that seems to be related with the Levantine Iron Age Anomaly (LIAA), the most prominent anomaly of the geomagnetic field of the last three millennia. Two intensity maxima were observed at Iberian coordinates, the oldest around 750 BC (associated with easterly declinations of around 23◦ ) and the second 275 yrs later (475 BC) with northerly directions. The related virtual axial dipole moment was up to 14 · 1022 Am2 for the oldest materials (750 BC) and reaching 16 · 1022 Am2 for the materials corresponding to the end of the Early Iron Age. In order to investigate the origin of the unusually high fluctuations of the palaeofield we have developed a new global geomagnetic field reconstruction, the SHAWQ-IronAge model, which is based on a critical revision of the global archeomagnetic and volcanic dataset. The new model provides an improved description of the evolution of the LIAA, which is related to a normal flux patch at the core-mantle boundary (CMB) below Arabian Peninsula clearly observed at around 950 BC. This flux patch expanded towards the north-west, while decreasing in intensity, reaching Iberia at around 750 BC. Around 600-500 BC, it underwent a revival below the European continent after that it seems to vanish in situ. © 2019 Elsevier B.V. All rights reserved.

1. Introduction

*

Corresponding author. E-mail addresses: mlosete@fis.ucm.es (M.L. Osete), [email protected] (A. Molina-Cardín), [email protected] (S.A. Campuzano), [email protected] (G. Aguilella-Arzo), [email protected] (A. Barrachina-Ibañez), [email protected] (F. Falomir-Granell), [email protected] (A. Oliver Foix), [email protected] (M. Gómez-Paccard), [email protected] (F. Martín-Hernández), [email protected] (A. Palencia-Ortas), [email protected] (F.J. Pavón-Carrasco), [email protected] (M. Rivero-Montero). https://doi.org/10.1016/j.epsl.2019.116047 0012-821X/© 2019 Elsevier B.V. All rights reserved.

Archaeomagnetic data provide important constraints on the local and global geomagnetic field variation during the historical past, contributing to the investigation of the geodynamo theory. Rapid and unusually high intensities during the Iron Age were observed in the Levantine region (Ben-Yosef et al., 2009; Shaar et al., 2011; 2016). These rapid variations, not seen in the present-day

2

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

and historical field, are of particular interest to constrain the full range of core dynamics. Ben-Yosef et al. (2009) and Shaar et al. (2011) named these geomagnetic anomalies characterized by short-lived, extremely high field maximum as “geomagnetic spikes”. Subsequent studies from Turkey (Ertepinar et al., 2012) and Georgia (Shaar et al., 2013) confirmed these high regional intensities around 900 BC. Shaar et al. (2016) proposed the term of “Levantine Iron Age Anomaly” (LIAA) to describe these features charactering the field surrounding the Levant region for a period of 350 yrs (∼1050 BC to ∼700 BC). During this period of generally high field, at least two geomagnetic “spikes” took place reaching Virtual Axial Dipole Moment (VADMs) values of 16–19 · 1022 Am2 (Ben-Yosef et al., 2017). Associated with these high intensity maxima directional deviations of more than 22◦ from axial dipole field direction were also observed (Shaar et al., 2016; 2018). Recent studies corroborate the occurrence of similar anomalies during the first Millennium BC from other areas around the world (see Korte and Constable, 2018 for a summary). High palaeointensities were observed at the Canary Islands (De Groot et al., 2015; Kissel et al., 2015), Azores (Di Chiara et al., 2014), central Europe (Hervé et al., 2017), Bulgaria (Kovacheva et al., 2014) and recently from Iberia (Molina-Cardín et al., 2018), that occurred at around 600 BC, that is, a few centuries after the LIAA. Older palaeointensity maxima have also been found in Korea (Hong et al., 2013) and eastern China (Cai et al., 2017). Livermore et al. (2014) indicated that the reported rates of secular variation from the Levant are not compatible with commonly accepted core-flow dynamics. Davies and Constable (2017) suggested that the source of the LIAA could be a flux patch at the Core Mantle Boundary (CMB) first growing in place and then migrating north and westward. However, recently Korte and Constable (2018) using spherical harmonic modelling, suggested that the spikes do not migrate but they are originated from intense normal flux patches growing and decaying mostly in situ, combined with stronger and more variable dipole moment than shown by previous global geomagnetic field models. A recent study on Mediterranean sediments Béguin et al. (2019), suggests that the LIAA moved westwards at a rate of 15–30◦ in 1000 yr, while decaying in field strength from 15 to 11 · 1022 Am2 in the Levant region. They observed low intensities from one core from Alborán (south Spain) for the last 3 millennia suggesting a western limit for the LIAA at least until 250 BC. In contrast, a recent full-vector Palaeosecular Variation Curve (PSVC) for Iberia spanning the last three millennia based on high quality archaeomagnetic data (Molina-Cardín et al., 2018) observed high intensities at around 600 BC related to eastern declinations of around 20◦ . They also compiled several PSVC from the Middle East to Iberia (relocated to Madrid latitude) obtaining East-West Hovmöller diagrams (Hovmöller, 1949) for each paleomagnetic element to investigate the evolution of the palaeofield, at the Earth’s surface. They suggested that the intensity maximum observed in Iberia at around 600 BC seemed to be related to the LIAA feature moving westwards. Iberia occupies a unique situation to investigate the possible westward migration of the LIAA. The Iberian PSVC (Molina-Cardín et al., 2018) is well constrained for the Islamic period (720-1492 AD) and for the Late Iron Age, ibero-roman and roman periods (300 BC-400 AD). However, the Early Iron and Late Bronze (1300550 BC) segment is based on few quality data and needs some improvement. So far, few directional data prior to the Late Iron Age have been investigated from Iberia and only few palaeointensity data meet the currently standard quality criteria. In this study, we present a new set of archeomagnetic data (14 directions and 10 palaeointensities) coming from 17 combustion structures of 6 archaeological settlements located in the Castellón

province (Eastern Spain). Investigated materials were dated by radiocarbon and archaeological evidences between the 8th and the 5th centuries BC. These data along with the previous archaeomagnetic quality data allow us to investigate the secular variation rates during the Early Iron Age in Iberia and to test the westward displacement of the LIAA. 2. Archaeological context and sampling A total of 129 oriented hand samples were obtained from 17 combustion structures from Eastern Spain (see Fig. 1). Six archaeological sites were investigated from the Castellón province: Puig de la Nau (PUIG); Puig de la Misericordia (PUM); Monte Calvario (CAL); Los Morrones de Arenoso (MRR); El Tossal de la Vila (TOS) and El Mortorum (MRT). In this area the Early Iron Age began at the end of the 8th century BC and ended around the middle of the 6th century BC, when the Iberian Culture began. The Old Iberian period developed during the entire second half of the 6th century BC and the first third of the 5th century BC. Detailed information on the archaeological sites is provided in the supplementary material. Sampling was mostly performed by means of gypsum bandages and plaster of Paris. Flat plaster surfaces adhering on selected sampling points were prepared and oriented by a sun and/or a magnetic compass. In few cases the structures sampled presented a flat and non-porous surface that hindered the adherence of the plaster; in these cases, we used the upper flat face as a reference plane. In the laboratory, hand samples were consolidated by impregnation in sodium silicate (waterglass) and were cut in 8 cm3 cubic specimens. In all cases samples were taken from the floor of the combustion structures. The structures sampled were small circular or oval shape cooking hearths (with diameter or semi-major axis ranging from 30 cm up to 120 cm). 3. Methods Palaeomagnetic measurements were conducted at the palaeomagnetism laboratory of Complutense University of Madrid (Spain). Rock magnetic studies include thermomagnetic curves and hysteresis analyses that were performed with a KLY3 kappabridge, a CS3 temperature control unit (Agico Ltd.) and a coercivity meter (Jasonov et al., 1998) with a maximum applied field of 500 mT. Additional thermomagnetic curves were performed by a Magnetic Measurements variable field translation balance (MMVFTB). Natural remanent magnetization (NRM) was measured on a modified Minispin spinner magnetometer (Molspin) and on a 2G Enterprises SRM755R superconducting rock magnetometer. Stepwise thermal and alternating field demagnetizations were conducted by MMTD80 or MMTDSC (Magnetic Measurements) thermal demagnetizers and by the AGICO AF demagnetizer LDA-5 (up to 200 mT) or by using the demagnetizer incorporated to the 2G Enterprises SRM755R superconducting rock magnetometer. After pilot experiments (carried out on two specimens per site), specimens were subjected to stepwise alternating field (AF) demagnetization in 10-15 steps up to the maximum peak AF of 100–200 mT and by thermal demagnetization (TH) from 50-100 ◦ C to temperatures up to 550-600 ◦ C. One step was included (at 350-450 ◦ C) in all thermally demagnetized specimens to check the possible influence of the anisotropy of the thermoremanent magnetization (ATRM) in the estimations of palaeomagnetic directions (PalenciaOrtas et al., 2017). Low field bulk susceptibility was measured after each step of thermal demagnetization using a KLY3 Kappabridge (Agico) or a MS3 susceptibility meter (Bartington) in order to detect magneto-mineralogical alterations during the experiments.

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

3

Fig. 1. Locations of the studied archaeological sites (a) and representative studied structures (b-d).

In addition to directional studies, archaeointensity determinations were carried out following the classical Thellier palaeointensity method including pTRM-checks (Thellier and Thellier, 1959). This method is based on the comparison of the NRM lost and the partial thermoremanent magnetization (pTRM) gained in a known laboratory field (50 μT), applied along the Z axis, at progressive higher temperature steps. The ATRM tensor was calculated for each specimen from the acquisition of a TRM in six different directions at a heating step between 350 ◦ C and 475 ◦ C, depending on thermal demagnetization behaviour of the specimens, following the method described by Veitch et al. (1984) and Chauvin et al. (2000). 4. Results 4.1. Rock magnetic results Initial NRM intensities and magnetic susceptibilities χ varied between 2.9 · 10−4 and 7.4 A/m and between 4.7 · 10−6 and 2.1 · 10−2 (SI) respectively. The lowest NRM values were observed in specimens from PUIG2, MRT1 and MRT2 structures, probably indicating that these specimens did not reach high temperatures. NRM vs. χ is plotted in Fig. 1S (supplementary material) together with Koenigsberger ratios (Qn ). High NRM and χ values, typical for well-baked argillaceous materials, are observed in most structures except in the previous mentioned sites. Especially well grouped and values of Qn around 10 are obtained on specimens from MRR2, PUM and CAL1 sites. Magnetic parameters from other structures are more scattered reflecting, probably, the position of the sample

in relation with the source of heat. Values of Qn larger than one support a thermoremanent origin of the NRM acquired during the use of the hearths and kilns. The highest values of Qn (up to 35) were obtained in specimens from TOS4 site. Magnetic hysteresis cycles are very similar in all cases, they display closed hysteresis loops that saturate at about 150 mT and with coercivities about 6-12 mT (Fig. 2a). Exceptionally, the structure MRT1 exhibited a non-saturated cycle at 500 mT. This hearth was badly preserved and displayed characteristics indicating that the material sampled did not reach high temperatures: very low intensities and non-consistent NRM directions. This site was excluded (see section below). Hysteresis parameters according to Dunlop (2002) would indicate a PSD range for magnetite (Fig. 2c). Corresponding thermomagnetic curves (Fig. 2d) show a dominant phase with a Curie point between 500 and 575 ◦ C. This magnetic phase is probably a Ti-poor titanomagnetite/titanomaghemite. In many cases another magnetic phase is also observed associated with a slight bulge at low temperatures (ranging from 100 ◦ C to 300 ◦ C) in the heating segment of the thermomagnetic curve. This phase is unstable and is not observed in the cooling curve. It is not clear what mineralogy could be responsible for this low coercivity and unstable magnetic phase. Thermomagnetic curves present different degrees of irreversibility, which makes it necessary to control the existence of chemical transformations during the thermal demagnetization, especially, during palaeointensity experiments.

4

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

Fig. 2. Representative rock-magnetic results. hysteresis curves (a) (red corrected by the paramagnetic contribution, black uncorrected), (b) SIRM and back field; (c) Day plot of hysteresis ratios; (d) Susceptibility versus temperature Heating/cooling is shown in red/blue. (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)

4.2. Archaeomagnetic directions Most specimens presented a good magnetic behaviour during AF and TH demagnetization. After demagnetization of a small viscous component, which is usually fully removed at 100-150 ◦ C or 5-10 mT, a single component that demagnetizes towards the origin is well isolated (Fig. 3 a-f). The characteristic remanent magnetization (ChRM) is fully demagnetized at 80-100 mT fields. This confirms the dominant presence of magnetically soft minerals. Maximum unblocking temperatures ranged from 450 up to 575 ◦ C. Well-grouped ChRM directions were obtained in most sites. The ATRM correction was small and it does not alter statistically mean directions. Consequently, mean ChRM directions were computed together from specimens subjected to either AF, TH and Thellier treatments. Results are summarized in Table 1 and plotted in Fig. 2S (supplementary material). In four structures (MRT2, CAL2, MRR1 and TOS2) ChRM directions were more scattered (α95 > 5◦ ), probably related to the degree of alteration of these structures. In addition, some of them showed evidences of insufficient heating. These structures probably lost their outer layer (the layer that suffered the highest temperatures). In addition, a lower number of specimens could be obtained from MRT2, CAL2 and MRR1 sites due to the frailty of these structures. In spite of this, consistent directions with the other structures sampled at the same archaeological site were obtained. We maintained the mean directions obtained in these sites, although their statistical parameters were not as good as the previous ones. Few specimens were sampled at PUIG3 and MRT1 sites due to the bad preservation state of these structures. The specimens behaved very fragile during thermal demagnetization. In addition,

unstable directions during pilot AF and TH treatments revealed that these structures did not probably reach high temperatures (this is also supported by the low values of Qn∼1). These sites were rejected. Only in TOS4 site two different components were found. An intense low coercivity component, demagnetized at 25-40 mT (Fig. 3h), presented a flat and well grouped direction (N = 10, D = −39.7◦ , I = 9.8◦ , K = 112, α95 = 4.6◦ ). The origin of this component is uncertain. This nearly flat direction is unlikely to represent the Earth’s magnetic field at the site latitude. After demagnetization of this low coercivity component a northerly direction is revealed by AF cleaning (Fig. 3h). By thermal demagnetization the first component could be isolated up to 300-450 ◦ C (N = 9, D = −39.2◦ ; I = 12.6◦ , K = 161, α95 = 4.0◦ ) but the high temperature component was not properly isolated. An apparently linear segment going throughout the origin was observed in some of the orthogonal plots (e.g. Fig. 3g) at specimen level. However, when computing the mean direction of the higher unblocking temperature we obtained a non-fisherian distribution, revealing that the high temperature component was not properly isolated by thermal treatment. Initial NRM intensities were high (in most cases ranging from 0.8 · 10−1 up to 5 A/m) as well as Qn ratios (Qn>10). These high values are related to the contribution of this low coercivity/low temperature phase. Horizontal magnetic fields are naturally produced by lightning strikes (e.g. Verrier and Rochette, 1999; Wasilewski and Kletetschka, 1999; Tauxe et al., 2003) but they are usually associated to very high NRM intensities and high Qn ratios. In our case, initial NRM intensities and Qn ratios from TOS4 site are higher than from other structures, but the obtained values are not anoma-

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

5

Table 1 Summary of palaeomagnetic information for the studied structures. Structure: label of the combustion structure; Site: Name of the studied archaeological site; Lat, Lon, coordinates; Age (yrs), assigned age; Dating method (C14, radiocarbon; Arch, archeological constraints); Nd, number of specimens retained for mean direction, D, mean declination; I, mean inclination; k, precision parameter; α95 , semiangle of 95%. confidence limit; Ni number of specimens retained for mean archaeointensity calculation; F, mean intensity and standard deviation (after ATRM correction). VADM, Virtual axial dipole moment of each structure. *: Data excluded (see text for details). Structure

Site

Lat (◦ N)

Lon (◦ E)

Age (BC)

Dating method

Nd

D

I

k

α95

Ni

F (μT)

VADM (1022 Am2 )

PUIG-1 PUIG-2 PUIG-3 PUIG-4 PUIG-5 PUIG-6 MRT2 PUM CAL1 CAL2 MRR-1 MRR-2 TOS1 TOS2 TOS3 TOS4* MRT1

Puig de la Nao Puig de la Nao Puig de la Nao Puig de la Nao Puig de la Nao Puig de la Nao El Mortorum P. Misericordia El Calvario El Calvario Los Morrones Los Morrones Tossal de la Vila Tossal de la Vila Tossal de la Vila Tossal de la Vila El Mortorum

40.42 40.42 40.42 40.42 40.42 40.42 40.15 40.47 40.03 40.03 40.19 40.19 40.27 40.27 40.27 40.27 40.15

0.43 0.43 0.43 0.43 0.43 0.43 0.05 0.47 −0.55 −0.55 −0.54 −0.54 −0.02 −0.02 −0.02 −0.02 0.05

412.5±37.5 412.5±37.5 412.5±37.5 412.5±37.5 412.5±37.5 412.5±37.5 550 ± 25 650± 50 625± 25 625± 25 650± 50 650± 50 775± 25 775± 25 775± 25 775± 25 1250 ±150

C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch C14+arch

14 14 4 23 13 14 7 18 19 7 7 23 20 6 34 13 –

−2.6 −1.3

58.2 57.7 – 60.8 61.6 61.4 63.1 60.9 58.7 63.9 56.4 62.9 68.1 57.9 55.1 14.6 –

426 336 – 216 391 457 50.3 340 369 107 123 327 179 101 45.5 182 –

1.9 2.2 – 2.1 2.1 1.9 8.6 1.9 1.7 5.8 5.5 1.7 2.4 6.7 3.7 3.1 –

4 3 – 10 6 5 – 5 6 – – 5 4 – 4 4 –

79.5±5.1 72.4±5.8 – 97.1±9.1 85.9±10.1 80.0±8.0 – 64.3±5.4 62.1±3.5 – – 55.4±6.5 71.5±3.8 – 77.9±1.0 118.0±18 –

13.9 12.7 – 17.0 15.0 14 – 11.2 10.9 – – 9.7 12.5 – 13.6 20.6 –

–

−2.9 −3.4 3.6 −4.3 7.8 18.8 18.3 21.8 11.5 11.3 17.3 10.5 −40.2 –

Fig. 3. Representative orthogonal vector projections of the remanent magnetization of studied samples during stepwise thermal (a, c, e, g) and alternating field (b, d, f, h) demagnetization. TOS4 site exhibited two magnetization components more clearly identified by alternating field demagnetization (h) than by thermal cleaning (g).

lously high. During the archaeological excavation it was observed that the last use of this structure was not as a hearth or a kiln but a working place. We cannot propose a physical mechanism (different than a distant lightning strike) capable of producing a remagnetization with the observed characteristics. As we do not fully understand the origin of the remanence this site was not considered. In summary a total amount of 14 new archaeomagnetic directions have been obtained for the Early Iron Age in Iberia. Ten struc-

tures provided well grouped directions (α95 ≤ 3.7◦ ) and four structures yielded lower confidence parameters (α95 ≥ 5◦ ) but mean directions were consistent with similar structures from the same archaeological site. 4.3. Archaeointensities The 42% of investigated specimens provided good archaeointensity results. This low rate of success in palaeointensity determina-

6

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

Fig. 4. Representative Arai plots from studied sites. Filled circles represent the chosen interval for palaeointensity calculations. A small Zijderveld plot of each selected specimen is also represented showing the NRM structure.

tion is a common problem in old Iberian hearths, often related with the fragility of the samples, the construction techniques, the state of preservation of the thin heated slab, the mineralogy and/or the burial conditions. In addition, the temperatures reached in some of the domestic kilns and hearths were not very high and during the heating treatment in the laboratory alterations occurred preventing an adequate determination of the palaeointensity. Representative Arai plots are plotted in Fig. 4. 48 of the 114 investigated specimens were selected for palaeointensity determina-

tion (Table 1S). They showed positive pTRM-checks and palaeointensity was calculated with at least 6 points from the Arai plot (in most cases from more than 9 points). The NRM fraction used for archaeointensity determination for the majority of specimens ranges between 0.6 and 0.9 corresponding with quality factors q values between 10.0 and 39.8 which satisfied the commonly applied selection criteria (e.g. Gómez-Paccard et al., 2016). There are only 2 specimens that exhibited lower q and f values, but they behave similarly to the rest of the specimens from the same struc-

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

ture and provided consistent results. A mean archaeointensity was also computed excluding these specimens, but result does not differ from the previous calculations considering the whole group and were maintained. Intensity results at specimen level are presented in Table 1S before and after the ATRM correction. It can be seen that the ATRM correction was small, less than 3% in most of the structures. High palaeointensity determinations (72.4-97.1 μT) were obtained from the 5 investigated combustion structures from PUIG site (Fig. 4 a-g). They correspond to values of the past geomagnetic field ranging from 12.7 · 1022 Am2 up to 17 · 1022 Am2 in terms of virtual axial dipole moment (VADM). Maximum values are about the double of the present day VADM. Acceptable palaeointensity results were obtained from 5 of the 7 specimens investigated from PUM structure. Linear segments in Arai plots were obtained from 100-250 ◦ C Thellier step up to 450525 ◦ C (Fig. 4f). Mean palaeointensity was 64.3 ± 5.4 μT (Table 1S). The CAL1 structure showed a very good magnetic behaviour. Consistent palaeointensities were obtained from all analyzed specimens. After demagnetization of a low temperature component linear segments were obtained from 200-300 ◦ C Thellier step up to 525 ◦ C (Fig. 4g). Easterly declinations and moderate palaeointesities (62.1 ± 3.5 μT) were derived from this structure. In contrast to CAL1, the CAL2 hearth was not well preserved. Pilot specimens thermally demagnetized indicated that this structure was not a good candidate for palaeointensity studies due to the fragility of the samples which were easily destroyed by heating. Palaeointensity determinations were attempted on 6 specimens from MRR2 combustion structure; 5 of them provided useful results. Linear segments in the Arai plots are observed (Fig. 4i) from 200-250 ◦ C. This structure provided the lowest intensity value obtained in this study (55.4 ±6.5 μT). Curved Arai plots associated with more than one component were observed in pilot palaeointensity specimens from MRR1 site (Fig. 4h). Thermally demagnetized samples indicated that temperatures reached by the recovered materials were low and with a high gradient within the structure. In addition, some samples were altered during thermal demagnetization. This site was therefore rejected for palaeointensity estimations. Palaeointensity determinations were attempted on 13 specimens from TOS1 domestic kiln, only 4 were retained. Rejected specimens presented curved Arai plots. In contrast, selected specimens presented a good linear trend (Fig. 4j) and acceptable statistical parameters. The past intensity derived from selected specimens was 71.5 ± 3.8 μT. From TOS2 structure few hand samples could be recovered due to the fragility of the hearth. Studied specimens did not provide quality intensity results. Curved Arai plots or evidences of alterations produced by heating were detected. Linear NRM-TRM diagrams have been observed up to 425– 550 ◦ C in 10 specimens (16 investigated) from the domestic TOS3 hearth (Fig. 4k). Mean intensity computed from these specimens was high (88.0 ± 9.4 μT). Although not all samples meet the required quality criteria (4 specimens had a q less than 10 and exhibited high palaeointensities and 2 presented similar intensitiy values, Table 2S). Excluding these specimens, a very well grouped palaeointensity was obtained (77.9 ± 1.0 μT) that is very consistent with TOS1 structure. We kept this last value. The mean palaeointensity value derived from rejected specimens was 94.6 ± 5.0 μT (information at specimen level is provided in Table 2S). As previously explained, a different behaviour was observed in TOS4 combustion structure. Two directional components were observed. Most of them presented curved Arai plots or had evidences of alteration. From the Thellier experiment, we only could estimate the palaeointensity of the low temperature component (after correction by the higher temperature phase) in few specimens that

7

provided a linear Arai plot (Fig. 4l). Very high palaeointensities were obtained. We have estimated a mean palaeointensisity of 118.0 ±18.1 μT from 4 specimens (showed in Table 2S, supplementary material). But, as mentioned, if the first component is an IRM induced by lightning strikes, the obtained palaeointensities have no significance in terms of PSV (further analyses comparing the Thellier results on the NRM and on a TRM and IRM imparted at laboratory are provided in supplementary material, Fig. 3S). At this point, we would like to remark the importance of studying in situ combustion structures. Because the materials researched in this study were oriented it was possible to doubt the origin of the first component in samples obtained from the TOS4 structure. With non-oriented structures, such as ceramic fragments, this component could have been interpreted as a partial thermoremanence (pTRM) and the derived palaeointensities could have been considered. This fact highlights the importance of studying in situ hearths and kilns, when available. In summary, from the 114 specimens analyzed in this work, 48 have been retained and used for calculating 10 new palaeointensities for Early Iron Age from Eastern Iberia which are summarized in Table 1. The ATRM was computed at specimen level (Table 1S). Low ATRM corrections were found. Differences between the uncorrected and ATRM corrected values were insignificant (≤3%). 5. Discussion 5.1. The Iberian database for the Late Bronze and Iron Ages The new archaeomagnetic data obtained in this study were relocated to Madrid coordinates (40.4◦ N latitude, 3.7◦ W longitude) using the virtual geomagnetic pole method (Noel and Batt, 1990), and compared with previous palaeomagnetic information from Iberia and surroundings areas (within a radius of 900 km centred at Madrid) for the Late Bronze-Iron Ages (Fig. 5). As already mentioned, there were few archaeomagnetic directional studies carried out in Iberia for the Late Bronze/Early Iron Ages. To date, only four directional data coming from well-dated archeological structures have been published (Catanzariti et al., 2008; Palencia-Ortas et al., 2017; Molina-Cardín et al., 2018; orange dots Fig. 5a-b). Recently, four archeological hearths from the Early Iron Age complex of Sant Jaume (Eastern Spain) of have also been archeomagnetically dated (Gómez-Paccard et al., 2019) but since the age of this site is still under debate we decide to discard here the directions proposed. We provide 14 new directional data (red dots Fig. 5 a-b) form this time period (Late Bronze/Early Iron). Recent directional data from Iberia of the Late Iron Age (PalenciaOrtas et al., 2017) have also been included in Fig. 5 (solid orange dots). In contrast to directional studies, a significant number of archaeointensity investigations have been carried on Iberian remains from 1400 BC up to 0 AD (Fig. 5c), but, unfortunately, they do not fulfil the current quality protocols (e.g. those described in PavónCarrasco et al., 2014b). 82 archaeointensity data from Spain and Portugal were determined on pottery sherds based on single specimens (Burakov et al., 2006; Nachasova et al., 2007; Nachasova and Burakov, 2009) which are represented on Fig. 5c (open orange dots). Only 7 recently published palaeointensity data meet the current quality standards (Molina-Cardín et al., 2018; Osete et al., 2016) solid orange dots Fig. 4c-d). We provide 10 new paleointensity data (red dots). As the number of palaeomagnetic information from Iberia is still very scarce, we have also considered data from adjacent regions (from sites located within 900 km radius from Madrid) available from the current GEOMAGIA database (Brown et al., 2015) as well as from recent publications. Data were also relocated to Madrid coordinates (blue dots). We included only data with an age

8

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

Fig. 5. Declination (a), inclination (b) and intensity (c) data from Iberia and southern France (within a radius of 900 km centred at Madrid), for the Late Bronze-Iron Ages, relocated to Madrid coordinates. Red dots: data from this study, orange and blue, data from Iberia and southern France respectively. Solid/open circles represent high/low quality data. d) Intensity data with mean values computed from individual specimens (see text for explanation).

error ≤ 250 yrs. The considered data base is summarized in Table 3S (supplementary material). 114 archaeointensities and 54 directional data were initially considered. From a simple visualization of Fig. 5a-b, a consistent directional evolution pattern is observed. Easterly declinations characterize the 900-550 BC time span, followed by northerly or northwest directions during the end of the Early Iron Age and the Late Iron/Iberian period. In inclination one clear outlier is observed (Ruiz-Martinez et al., 2008) which is automatically excluded during the PSVC computation procedure (see next section). In contrast to the directional database, the scatter of the archaeointensity data is very puzzling. This is a common problem observed at a global scale. This problem has been addressed in different ways. Some authors impose quality criteria to reject data that do not meet these criteria (e.g. Hervé et al., 2013b) and others believe that with a sufficient number of data, the erroneous data will be eliminated by a statistical calculation process (e.g. Tema and Kondopoulou, 2011). Here we will address the problem following a revised classification criteria based on the previous proposal of Pavón-Carrasco et al. (2014b). Archeointensity data are classified in two categories: high or low quality. The high quality archaeomagnetic data correspond to categories A and B and C1 of Pavón-Carrasco et al. (2014b). That is: 1) palaeointensities must be obtained by the original or derived Thellier-Thellier method with pTRM checks; with TRM anisotropy correction if derived from potential high-anisotropic objects (such as potteries or tiles) and 2) mean archeointensities should be derived from at least three specimens. The lower quality group includes the rest of data. With these criteria 42 palaeointensitity data from the 1400 BC0 AD interval were considered of high quality (10 of them coming from this study) and 87 data were labelled as low quality data. Most of low quality data are palaeointensities derived from pot-

tery fragments from Spain and Portugal (Burakov et al., 2006; Nachasova and Burakov, 2009) but also a few of them (5) were obtained from French and Moroccan sites. The first problem of the low quality group is that the data were not corrected by the effect of the anisotropy of thermoremanent magnetization (ATRM) on archaeointensity estimations. This effect can be particularly important for the materials analyzed (ceramics) and reach, in some cases, very high correction factors of more than 80% of the determined field strength (e.g., Chauvin et al., 2000; Genevey et al., 2008; Osete et al., 2016). In well preserved Iberian (Portugal) hearth and kilns from the Early Iron Age, Molina-Cardín et al. (2018) observed differences between the uncorrected and ATRM corrected palaeointensity values that reached up to 11.4%. Higher values of the ATRM correction (up to 86%) were found in younger structures from the Iberian period in Iberia (Osete et al., 2016). The second problem of the low quality group is that the determination of many archaeointensities from Iberia were based on one single specimen. In order to homogenize the Iberian palaeointensity database and not give, artificially, more weigh to archaeointensities derived from one specimen we have recalculated a mean palaeointensity from those specimens studied by Burakov et al. (2006) and Nachasova and Burakov (2009) coming from the same site and having similar dates. Data derived from one single specimen were not considered. With these new calculations the total number of archaeointensity data is 68 (42 of high quality) for the period comprised between 1400 BC and 0 AD (Fig. 5d). 5.2. The full vector Iberian PSVC for Late Bronze/Iron Ages To analyze the present Iberian archaeomagnetic database and investigate more in detail the effect of including quality categories

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

9

Fig. 6. a) Full vector Iberian PSVC based on high quality data, together with the data used for its construction. It has been represented with a 2-sigma error band. Open circles represent the outliers rejected for curve calculation. b) Iberian palaeosecular variation curves obtained with different input datasets (dark green: only quality data, light green: all data from the database, green: all data with averages of data based on one specimen). c) Iberian PSVC compared with other curves and synthetic curves derived from regional and global models.

in determination of Iberian PSV curves we have computed 3 full vector PSVCs. Differences in the curves are related to the input

database used. In the first case (Fig. 6a), only quality data were considered. Fig. 6b shows the previous one (dark green), compared

10

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

with the curves obtained including all data with the new means computed (green) and with all data published in their original version (light green). The methodology followed to compute the time-continuous full vector PSV curves for Iberia was similar that the one used in Molina-Cardín et al. (2018). The three elements of the paleofield (i.e. declination D, inclination I and intensity F) were fitted all together to define a full-vector of the geomagnetic field at the given point (i.e., Madrid’s coordinates). The time evolution is represented by cubic b-splines with knot points every 50 yrs. In the fitting approach, a rejection criterion is applied to detect possible outliers. All data with residuals higher than 3 times the average root mean square error were rejected. The palaeofield during the Iron Age in Iberia is characterized by two intensity maxima (at 750 BC and at 475 BC) and a rapid directional variation. Differences between the three curves are related with the value reached in the first intensity maximum (at 750 BC). The curve based on the whole database (light green curve in Fig. 6b), which includes non-quality data obtained from ceramics, seems to overestimate the palaeofield. It could be due, probably, to the ATRM effect not corrected in these specimens. We prefer to be cautious and consider for the reference Iberian Iron Age PSVC the one based on quality data (Fig. 6a). Easterly declinations characterize the paleofield in Iberia during the Late Bronze/Early Iron Ages (900 BC) up to around 500-400 BC. With a maximum deviation of around 23◦ reached at ∼750 BC. Coeval with this directional deviation from an axial dipole the first maximum of palaeointensity reaching 73 μT is also evidenced. This palaeointensity maximum is followed by a decrease in intensity (62 μT) at ∼600 BC. A younger maximum was observed at ∼475 BC, reaching the highest palaeofield values of the Iron Ages (84 μT). This maximum is related with northerly declinations. Therefore, the two intensity maxima are clearly distinguishable from the archaeomagnetic point of view. They correspond to two different palaeofields. Declination moves rapidly from eastern declinations to the north/north-west at a mean rate of 8◦ /century. Westerly declinations of around −5◦ are observed at 350 BC. Inclinations are increasing during the first part of the Iron Age up to 500-475 BC, when a maximum value of ∼65◦ is reached coeval with the second intensity maximum. Rapid changes of the geomagnetic field strength are observed in this area during the Early Iron Ages. Sudden increases of 1016 μT/century followed by rapid decreases of about 4-9 μT/century occurred in NE Iberia during the 10-5th centuries BC. 5.3. Comparison with other western European PSVCs and palaeomagnetic reconstruction models The proposed reference curve is compared with other local and model-derived PSVCs. The previous Iberian full vector PSVC (Molina-Cardín et al., 2018), the French curve for directions (Hervé et al., 2013a), the central Europe paleointensity curve (Hervé et al., 2017), both relocated to Madrid coordinates, and the model predictions derived from regional and global models are represented in Fig. 6c. We used the predictions given by the SCHA.DIF.3k regional model for directions (Pavón-Carrasco et al., 2009) and the regional model B from Pavón-Carrasco et al. (2014b) for intensities and the global models HFM.OL1.A1 (Constable et al., 2016) and SCHA.DIF.14 (Pavón-Carrasco et al., 2014a). A general agreement is observed in declination with the exception of the HFM.OL1.A1 global reconstruction. This evidences that reconstructions including sedimentary data in the input database are too smoothed and they lose a part of the variability of the secular variation of the Earth’s magnetic field. However, when an archaeomagnetic and volcanic database is used in the gener-

ation of global or regional geomagnetic models (e.g. SHA.DIF.14k, SCHA.DIF.3k) the short wavelengths are better reconstructed. Most of predictions and data suggest that inclination values reached a maximum during the end of the Early Iron Age (600-500 BC). Major differences are observed in palaeointensity. The maximum observed around 600 BC in the last Iberian PSVC (MolinaCardín et al., 2018) is now subdivided into two maxima (at 750 BC and 475 BC) separated by a relative minimum at 600 BC. Our new curve is a refinement of the previous one, which was based on only two palaeointensity data to define the maximum. The VADM was 13.2 · 1022 Am2 for the first maximum (750 BC) and 14.9 · 1022 Am2 for the youngest maxima (475 BC). The regional reconstruction by Pavón-Carrasco et al. (2014b) shows two maxima during the Iron Ages as well; the first one at about 900 BC and the youngest at 650 BC. The central European palaeointensity curve (Hervé et al., 2017, in violet) exhibited also two maxima, but at different times: at 950 BC (the smallest one, with a VADM of ∼ 10.5 · 1022 Am2 ) and at 550 BC (the highest one, VADM: 14.0 · 1022 Am2 ). That is, a similar behaviour to the Iberian palaeointensity curve but slightly smoothed and displaced in time. A displacement of two centuries is observed for the oldest maximum and around 75 yrs for the youngest and most intense maximum. 5.4. Evolution of the Levantine Iron Age Anomaly (LIAA) revealed by the new paleomagnetic reconstruction SHAWQ-Iron Age model The observed high speed variations in directions and intensities seems to be related to the Levantine Iron Age Anomaly. As mentioned, Molina-Cardín et al. (2018) by performing a joint analysis of European quality data suggested that the LIAA could be envisaged as a local magnetic palaeofield feature drifting westward. In contrast, Korte and Constable (2018) from geomagnetic modelling considered that LIAA grew up and vanished in situ. Whether geomagnetic intensity spikes can move or are static or how big is the region affected by them is a current and important debate (e.g. Davies and Constable, 2017; Korte and Constable, 2018; Béguin et al., 2019). Here we have analyzed the complete data set for the Late Bronze-Iron Ages (1400 BC-100 AD) and constructed a global model weighted by the quality of the data, the SHAWQ-Iron Age model from 1300 BC to 0 AD (Figs. 7 and 8). The methodology used to construct this model was that described in Campuzano et al. (2019). It is based on the Spherical Harmonic Analysis (SHA) technique in space and the penalized cubic B-splines in time (Korte and Constable, 2005). Input data were weighted according to the quality criteria described above (5.1 section). The maximum degree of the harmonic expansion was 10 and the temporal knot points was fixed every 25 yrs. In order to estimate model uncertainties a bootstrap technique (similar to that described in Korte and Constable, 2008) was applied. More information on the details of the model are described in the supplementary material (Figs. 47S and Table 4S). According to the SHAWQ-Iron Age palaeoreconstruction an increase in the dipolar moment is observed at 950 BC, reaching values of about 10 · 1022 Am2 (Fig. 7S). A value that is higher than the present day field, but not enough to explain the observations of the VADM described in the Levant region. Therefore, the LIAA is clearly linked to a non-dipolar feature of the geomagnetic field. Fig. 7 shows that the observed high intensities at the Earth’s surface (the LIAA) can be associated to a flux patch of normal polarity (NFP) of the radial field (Br) at the CMB located beneath the Arabic Peninsula at 1100-900 BC (Fig. 7, left). The evolution of this NFP can be followed more in detail visualizing the Movie 1S (supplementary material). Due to the non-dipolar nature of the LIAA we have also analyzed the non-dipolar terms of the geomagnetic field to enhance its characteristics. Fig. 8 shows the 100 yrs aver-

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

11

Fig. 7. Evolution of the geomagnetic elements at the Earth’s surface together with the radial field at the CMB from the SHAWQ-Iron Age reconstruction. Left columns: snapshots of declination, inclination and intensity plotted each 200 yr. Right column: mean radial field (Br) at the CMB averaged every 200 yr using 200 yr window widths. A zoom on the European-African-Levant region is also shown. Note that the last average plotted corresponds to the whole studied period (1300 BC–0 AD).

aged evolution of the non-dipolar terms contributing to Br at the CMB (a movie is also given as the supplementary material, Movie 2S). At about 1100 BC an isolated NFP is located below the Levant region that grows up in situ reaching its highest value at around 950 BC. Then it expands towards the northwest while decreasing in intensity. This lobe of normal polarity affected the European region with a lower intensity around 800-700 BC. At around 650 BC its centre is located below Europe. A reinforcement of the anomaly below the European region is observed at 500-400 BC. Therefore, only a minor migration to the NW of the NFP centre is evidenced. An apparent migration is observed in Europe because the NFP expands to the NW while decreasing in amplitude. From 350 BC and in advance the NFP is vanishing in situ. This NFP was probably present before the LIAA was observed at the Earth’s surface. It seems that about 1200 BC it was already present at the CMB together with a reversed flux patch located close to the north. The presence of these two close lobes produced a cancellation, and therefore it was not observed as an anomaly on the Earth’s surface. More archaeomagnetic information is still needed to analyze more in detail the emergence of this anomaly. There is also a normal polarity patch below Indonesia. The two patch could be connected or not (more data from the region of India are still required).

6. Conclusions 10 new mean archaeointensities and 14 directional data were obtained from 17 archaeological structures (6 sites) sites from Castellon, Northeast Iberia. The new data were dated between 800 and 400 BC and complete the poor high-quality Iberian database for this period. A new full vector Iberian secular variation curve has been obtained showing a double oscillation of the geomagnetic field strength with two maxima, one during the 8th century BC reaching 72 μT and another at circa 500 BC of 82 μT. The two maxima are separated by a minimum of 62 μT near 600 BC. This evolution resulted in fast secular variation rates up to 16 μT/century. Coeval with these fluctuations in intensity, directions were deviated more than 20◦ from the axial dipole direction. Our results confirm that strong fluctuations of the geomagnetic field occurred over Europe during the Early Iron Ages. A new global geomagnetic field reconstruction has been performed, the SHAWQ-Iron Age model, which is based on a critical revision of the global archaeomagnetic and volcanic dataset. The new model provides an improved picture on the evolution of the Levantine Iron Age Anomaly (LIAA). The LIAA seems to be linked to a singular and important nondipolar anomaly at the CMB. A lobe of normal polarity was observed in the CMB below the Arabian Peninsula around 1000 BC, it extended to the northwest, while decreasing in intensity. This anomaly affected the European region around 750-500 BC, in-

12

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

Appendix A. Supplementary material Supplementary material related to this article can be found online at https://doi.org/10.1016/j.epsl.2019.116047. References

Fig. 8. Evolution of non-dipole radial field at the CMB computed for the 1300 BC – 0 AD time period averaged every 200 yrs using the SHAWQ-Iron Age model.

creased its intensity in this area around 500-400 BC, and then vanished in the Levant region at around 0 AD. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements Financial support was given by the Spanish Ministry of Education and Science (CGL2017-87015-P, and CGL2017-92285-EXP projects and FPU14/02422 grant) and from the Ministry of Economy and Industry and Competitivity (Contract of APO). SAC thanks the LIMADOU Scienza project (ASI, Italy). We would like to thank the excellent work during laboratory procedures of Aida Adsuar Vaquero and Enós Delso Calcerrada as well as the help of José C. Perez-Fuentes and Javier Carmona during the field work. All data or their sources (GEOMAGIA50.v3, http://geomagia.gfzpotsdam.de/index.php; HISTMAG, http://www.conrad-observatory. at/zamg/index.php/data-en/histmag-database) are indicated within the main text or the supporting information.

Béguin, A., Filippidi, A., de Lange, G.J., de Groot, L.V., 2019. The evolution of the Levantine Iron Age geomagnetic Anomaly captured in Mediterranean sediments. Earth Planet. Sci. Lett. 511, 55–66. https://doi.org/10.1016/j.epsl.2019.01.021. 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 (3–4), 529–539. https://doi.org/10.1016/j.epsl.2009.09.001. 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. 114 (9), 2160–2165. https://doi.org/10.1073/pnas. 1615797114. Brown, M.C., Donadini, F., Korte, M., Nilsson, A., Korhonen, K., Lodge, A., Lengyel, S.N., Constable, C.G., 2015. GEOMAGIA50. v3: 1. General structure and modifications to the archeological and volcanic database. Earth Planets Space 67 (1), 83. https://doi.org/10.1186/s40623-015-0232-0. Burakov, K.S., Nachasova, I.E., Mata, C., 2006. Geomagnetic field intensity in the first millennium BC from data on ceramics of the Los Villares archaeological monument (Spain). Izv. Phys. Solid Earth 42 (11), 942–950. 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 (1), 39–44. https://doi.org/10.1073/pnas.1616976114. 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 paleomagnetic reconstruction SHAWQ2k. Earth Planet. Sci. Lett. 512, 17–26. Catanzariti, G., McIntosh, G., Soares, A.M.M., Díaz-Martínez, E., Kresten, P., Osete, M.L., 2008. Archaeomagnetic dating of a vitrified wall at the Late Bronze Age settlement of Misericordia (Serpa, Portugal). J. Archaeol. Sci. 35 (5), 1399–1407. https://doi.org/10.1016/j.jas.2007.10.004. 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 (1–2), 111–136. 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., Constable, C., 2017. Geomagnetic spikes on the core-mantle boundary. Nat. Commun. 8, 15593. https://doi.org/10.1038/ncomms15593. De Groot, L.V., Béguin, A., Kosters, M.E., van Rijsingen, E.M., Struijk, E.L., Biggin, A.J., Hurst, E.A., Langereis, C.G., Dekkers, M.J., 2015. High paleointensities for the Canary Islands constrain the Levant geomagnetic high. Earth Planet. Sci. Lett. 419, 154–167. https://doi.org/10.1016/j.epsl.2015.03.020. Di Chiara, A., Tauxe, L., Speranza, F., 2014. Paleointensity determination from São Miguel (Azores Archipelago) over the last 3 ka. Phys. Earth Planet. Inter. 234, 1–13. Dunlop, D.J., 2002. Theory and application of the Day plot (Mrs/Ms versus Hcr/Hc) 1. Theoretical curves and tests using titanomagnetite data. J. Geophys. Res., Solid Earth 107 (B3). EPM-4. Ertepinar, P., Langereis, C.G., Biggin, A.J., Frangipane, M., Matney, T., Ökse, T., Engin, A., 2012. Archaeomagnetic study of five mounds from Upper Mesopotamia between 2500 and 700 BCE: further evidence for an extremely strong geomagnetic field ca. 3000 years ago. Earth Planet. Sci. Lett. 357, 84–98. 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). Gómez-Paccard, M., Osete, M.L., Chauvin, A., Pavón-Carrasco, F.J., Pérez-Asensio, M., Jiménez, P., Lanos, P., 2016. New constraints on the most significant paleointensity change in Western Europe over the last two millennia. A non-dipolar origin? Earth Planet. Sci. Lett. 454, 55–64. Gómez-Paccard, M., Rivero-Montero, M., Chauvin, A., i Rubert, D.G., Palencia-Ortas, A., 2019. Revisiting the chronology of the Early Iron Age in the north-eastern Iberian Peninsula. Archaeol. Anthropol. Sci. 1, 1–13. https://doi.org/10.1007/ s12520-019-00812-9. Hervé, G., Chauvin, A., Lanos, P., 2013a. Geomagnetic field variations in Western Europe from 1500 BC to 200 AD. Part I: directional secular variation curve. Phys. Earth Planet. Inter. 218, 1–13. https://doi.org/10.1016/j.pepi.2013.02.002. Hervé, G., Chauvin, A., Lanos, P., 2013b. Geomagnetic field variations in Western Europe from 1500 BC to 200 AD. Part II: new intensity secular variation curve. Phys. Earth Planet. Inter. 218, 51–65. Hervé, G., Faßbinder, J., Gilder, S.A., Metzner-Nebelsick, C., Gallet, Y., Genevey, A., Wittenborn, F., 2017. Fast geomagnetic field intensity variations between 1400 and 400 BCE: new archaeointensity data from Germany. Phys. Earth Planet. Inter. 270, 143–156.

M.L. Osete et al. / Earth and Planetary Science Letters 533 (2020) 116047

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. Hovmöller, E., 1949. The trough-and-ridge diagram. Tellus 1 (2), 62–66. https://doi. org/10.3402/tellusa.v1i2.8498. Jasonov, P.G., Nourgaliev, D.K., Burov, B.V., Heller, F., 1998. A modernized coercivity spectrometer. Geol. Carpath. 49 (3), 224–226. Kissel, C., Laj, C., Rodriguez-Gonzalez, A., Perez-Torrado, F., Carracedo, J.C., Wandres, C., 2015. Holocene geomagnetic field intensity variations: contribution from the low latitude Canary Islands site. Earth Planet. Sci. Lett. 430, 178–190. https:// doi.org/10.1016/j.epsl.2015.08.005. Korte, M., Constable, C.G., 2005. The geomagnetic dipole moment over the last 7000 years—new results from a global model. Earth Planet. Sci. Lett. 236 (1–2), 348–358. Korte, M., Constable, C.G., 2008. Spatial and temporal resolution of millennial scale geomagnetic field models. Adv. Space Res. 41 (1), 57–69. Korte, M., Constable, C.G., 2018. Archeomagnetic intensity spikes: global or regional geomagnetic field features? Front. Earth Sci. 6, 17. https://doi.org/10.3389/feart. 2018.00017. Kovacheva, M., Kostadinova-Avramova, M., Jordanova, N., Lanos, Ph., Boyadzhiev, Y., 2014. Extended and revised archaeomagnetic database and secular variation curves from Bulgaria for the last eight millennia. Phys. Earth Planet. Inter. 23, 79–94. Livermore, P.W., Fournier, A., Gallet, Y., 2014. Core-flow constraints on extreme archeomagnetic intensity changes. Earth Planet. Sci. Lett. 387, 145–156. https:// doi.org/10.1016/j.epsl.2013.11.020. Molina-Cardín, A., Campuzano, S.A., Osete, M.L., Rivero-Montero, M., Pavón-Carrasco, F.J., Palencia-Ortas, A., Martín-Hernández, F., Gómez-Paccard, M., Chauvin, A., Guerrero-Suárez, S., Pérez-Fuentes, J.C., McIntosh, G., Catanzariti, G., Sastre Blanco, J.C., Larrazabal, J., Fernández Martínez, V.M., Álvarez Sanchís, J.R., Rodríguez-Hernández, J., Martín Viso, I., Garcia i Rubert, D., Pérez-Fuentes, J.C., 2018. Updated Iberian archeomagnetic catalogue: new full vector paleosecular variation curve for the last three millennia. Geochem. Geophys. Geosyst. 19 (10), 3637–3656. https://doi.org/10.1029/2018GC007781. Nachasova, I.E., Burakov, K.S., Lorrio, A.J., 2007. Archaeomagnetic study of ceramics from the El Molon archaeological monument (Spain). Izv. Phys. Solid Earth 43 (10), 830–835. Nachasova, I.E., Burakov, K.S., 2009. Variation of the intensity of the Earth’s magnetic field in Portugal in the 1st Millennium BC. Izv. Phys. Solid Earth 45 (7), 595–603. Noel, M., Batt, C.M., 1990. A method for correcting geographically separated remanence directions for the purpose of archaeomagnetic dating. Geophys. J. Int. 102 (3), 753–756. Osete, M.L., Chauvin, A., Catanzariti, G., Jimeno, A., Campuzano, S.A., BenitoBatanero, J.P., Tabernero-Galán, C., Roperch, P., 2016. New archaeomagnetic data recovered from the study of celtiberic remains from central Spain (Numantia and Ciadueña, 3rd-1st centuries BC). Implications on the fidelity of the Iberian paleointensity database. Phys. Earth Planet. Inter. 260, 74–86.

13

Palencia-Ortas, A., Osete, M.L., Campuzano, S.A., McIntosh, G., Larrazabal, J., Sastre, J., Rodriguez-Aranda, J., 2017. New archaeomagnetic directions from Portugal and evolution of the geomagnetic field in Iberia from Late Bronze Age to Roman Times. Phys. Earth Planet. Inter. 270, 183–194. Pavón-Carrasco, F.J., Osete, M.L., Torta, J.M., Gaya-Piqué, L.R., 2009. A regional archeomagnetic model for Europe for the last 3000 years, SCHA.DIF.3K: applications to archeomagnetic dating. Geochem. Geophys. Geosyst. 10 (3). https:// doi.org/10.1029/2008GC002244. Pavón-Carrasco, F.J., Osete, M.L., Torta, J.M., De Santis, A., 2014a. A geomagnetic field model for the Holocene based on archaeomagnetic and lava flow data. Earth Planet. Sci. Lett. 388, 98–109. Pavón-Carrasco, F.J., Gómez-Paccard, M., Hervé, G., Osete, M.L., Chauvin, A., 2014b. Intensity of the geomagnetic field in Europe for the last 3 ka: influence of data quality on geomagnetic field modeling. Geochem. Geophys. Geosyst. 15 (6), 2515–2530. https://doi.org/10.1002/2014GC005311. Ruiz-Martinez, V.C., Pavón-Carrasco, F.J., Catanzariti, G., 2008. First archaeomagnetic data from northern Iberia. Phys. Chem. Earth Parts A/B/C 33 (6–7), 566–577. 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 (1–2), 297–306. https://doi.org/ 10.1016/j.epsl.2010.11.013. Shaar, R., Tauxe, L., Gogichaishvili, A., Rathert, M.C., Devidze, M., Licheli, V., 2013. Absolute geomagnetic field intensity in Georgia during the past 6 millennia. Latinmag Lett. 3, 1–4. 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. Shaar, R., Hassul, E., Raphael, K., Ebert, Y., Segal, Y., Eden, I., Vaknin, Y., Marco, S., Nowaczyk, N.R., Chauvin, A., Agnon, A., 2018. The first catalog of archaeomagnetic directions from Israel with 4000 years of geomagnetic secular variations. Front. Earth Sci. 6, 164. https://doi.org/10.3389/feart.2018.00164. Tauxe, L., Constable, C., Johnson, C.L., Koppers, A.A., Miller, W.R., Staudigel, H., 2003. Paleomagnetism of the southwestern USA recorded by 0–5 Ma igneous rocks. Geochem. Geophys. Geosyst. 4 (4). https://doi.org/10.1029/2002gc000343. Thellier, E., Thellier, O., 1959. Sur l’intensité du champ magnétique terrestre dans le passé historique et géologique. Ann. Geophys. 15, 285–376. Tema, E., Kondopoulou, D., 2011. Secular variation of the Earth’s magnetic field in the Balkan region during the last eight millennia based on archaeomagnetic data. Geophys. J. Int. 186 (2), 603–614. Veitch, R.J., Hedley, G., 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 (3), 359–373. Verrier, V., Rochette, P., 1999. Estimating peak currents at ground lightning impacts using remanent magnetization. Geophys. Res. Lett. 29, 1867. Wasilewski, P., Kletetschka, G., 1999. Lodestone: natures only permanent magnet. What it is and how it gets charged. Geophys. Res. Lett. 26 (15), 2275–2278.