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

The large and well-studied archaeological record of Israel offers a unique opportunity for collecting high resolution archaeomagnetic data from the past several millennia. Here, we initiate the ﬁrst catalog of archaeomagnetic directions from Israel, with data covering the past four millennia. The catalog consists of 76 directions, of which 47 fulﬁll quality selection criteria with Fisher precision parameter ( k ) ≥ 60, 95% cone of conﬁdence ( α 95 ) < 6 ◦ and number of specimens per site ( n ) ≥ 8. The new catalog complements our published paleointensity data from the Levant and enables testing the hypothesis of a regional geomagnetic anomaly in the Levant during the Iron Age proposed by Shaar et al. (2016, 2017). Most of the archaeomagnetic directions show < 15 ◦ angular deviations from an axial dipole ﬁeld. However, we observe in the tenth and ninth century BCE short intervals with ﬁeld directions that are 19 ◦ -22 ◦ different from an axial dipole ﬁeld and inclinations that are 20 ◦ -22 ◦ steeper than an axial dipole ﬁeld. The beginning of the ﬁrst millennium BCE is also characterized with fast secular variation rates. The new catalog provides additional support to the Levantine Iron Age Anomaly hypothesis.


INTRODUCTION
Despite decades of intense paleomagnetic research, many details of geomagnetic secular variations have still remained elusive. It is well accepted that secular variations average out globally to an axial dipole field over long geological timescales. Yet, many aspects concerning the spatial and temporal characteristics of secular variations remain unclear, especially when dealing with periods preceding direct human observational data. For example, while regional deviations of field direction from an axial dipole field are widely recognized, neither the degree limits nor the lifetime of these deviations are fully known. Today, the largest deviation from an axial dipole field, between 20 • and 30 • , occupies a confined area in the southern Atlantic associated with a low field intensity anomaly termed "South Atlantic Anomaly" (SAA) (Thebault et al., 2015). The question whether the SAA is a typical secular variations characteristic or, instead, a unique geomagnetic phenomenon is yet to be tested. Equivalently, it is not fully understood if rates and amplitudes of secular variations measured during the past few centuries (Jackson et al., 2000) also represent the characteristic behavior of the geomagnetic field in earlier periods. To fill these gaps in knowledge there is a growing need for reliable and precise paleomagnetic datasets in sub-millennial temporal resolution from periods preceding direct measurements of the geomagnetic field.
Archaeomagnetic data from in-situ archaeological structures, such as ovens, furnaces, kilns, and burnt buildings provide an excellent opportunity to capture the direction of the ancient field. When these structures cooled from high temperatures they acquired thermoremanent magnetization (TRM) parallel to the ambient field, thus preserving an instantaneous recording of the ancient field. In many cases, the age of the TRM can be precisely dated using radiocarbon, historical constraints, archaeological correlation, indicative pottery, coins, or a combination of these methods. In this perspective, the long, continuous, well-studied archaeological record of Israel offers a unique opportunity for archaeomagnetic research.
The global role of archaeomagnetic data from Israel is illustrated in Figure 1A, which shows a map of the published archaeomagnetic directional data available in the GEOMAGIA50 database (Korhonen et al., 2008;Brown et al., 2015) from the past four millennia. Our study area is located in an important geographic area that extends the densely scattered data from Europe to the southeast. To date, only several archaeomagnetic directions from Israel were published in journal articles (Aitken and Hawley, 1967;Segal et al., 2003;Shaar et al., 2016;Shahack-Gross et al., 2018). However, there are considerable unpublished data that are available only as unpublished theses (Segal, 2003;Hassul, 2015). The purpose of this work is to gather and compile all the available data from these sources and provide the first catalog of archaeomagnetic directions from Israel. To this end, we have collected the archaeological and the chronological information from the above sources, added new data from 15 additional sites, gathered all the raw paleomagnetic measurement data (if it exists) and translated them to a community standard MagIC format . The combined data were then re-analyzed using identical standards and selection criteria. The resulting catalog includes new secular variations data spanning the past four millennia.

Sites and Locations
The Israeli archaeomagnetic catalog is assembled from a collection of several sources: Two unpublished Masters theses: (A) Segal (2003) that also includes two sites published in Segal et al. (2003) and (B) Hassul (2015); Two published articles: (C) Shaar et al. (2016), and (D) Shahack-Gross et al. (2018); and (E) New data from 15 structures labeled hereafter "this study". Here, we also revise and augment the directions previously reported in Shaar et al. (2016) with new measurements. Therefore, the paleomagnetic interpretations reported here are slightly different from those previously published and replace the previous interpretations. In the catalog we follow the standard paleomagnetic hierarchy nomenclature and define "location" as a collection of structures from the same area (i.e., an archaeological site), "site" as a single structure (i.e., cooling unit), "sample" as an oriented piece from a given site, and "specimen" as the part from the sample that was measured.
Sites can be classified under one of the three categories shown in Figure 2. Cooking ovens (tabuns, Figure 2A) are rounded structures, typically about 1 m in diameter, frequently found in domestic settings. Although the ovens used fire as a heating source, burnt remains that could be used for radiocarbon dating usually did not survive. Hence, the age of the ovens are typically dated by the age of living stratum in which it was found. Burnt walls ( Figure 2B) are in-situ remains of large mudbrick structures, which were burnt as a whole during historical destruction events (e.g., Shahack-Gross et al., 2018). In the case of a large conflagration there may be a large amount of burnt organic material that can be directly dated and crosscorrelated with known historical military campaigns, leading to high precision dating, sometimes with an uncertainty of several years (e.g., Tel Megiddo; Finkelstein and Piasetzky, 2009), Tel Hazor (Sandhaus, 2013;Zuckerman, 2013), Tel 'Eton (Faust, 2008), Bethsaida (Arav, 2014), and Lachish (Ussishkin, 1990)]. Furnaces and kilns ( Figure 2C) are large industrial structures that were used to manufacture ceramics. The kilns can be dated from the type of the ceramics and the typology of other finds, such as coins.
The Supplementary Material provides a short summary of the archaeological context and the location of each of the structures. We note that some sites were collected from old and presently inactive excavations, where a revised inspection of the exact archaeological contexts and the corresponding ages is not always possible. Also, in some previously studied sites there might be an ambiguity regarding their precise age because the dates of the corresponding archaeological strata have been refined throughout the years as more archaeo-chronological data have been accumulated. Therefore, we clearly stress the need for a detailed review of the archaeological contexts and the ages in this catalog. This effort requires a thorough and detailed archaeological study that is beyond the scope of this research, but we make such a future investigation possible with the information given in the Supplementary Material.

Sampling and Lab Procedures
Oriented samples were obtained either by drilling standard paleomagnetic cores, 1" in diameter, with a portable electrical drill or by hand samples. In the case of hand samples, orientations of flat surfaces were measured and marked in the field before detaching the oriented samples from the structure (site). In some cases, the material was hardened in the field with epoxy before being measured and removed. Specimens were prepared from the hand samples by sawing small cubes from the sample and gluing them inside non-magnetic paleomagnetic plastic sampling boxes. Samples were oriented in-situ with a Brunton compass prior to detachment and collection for both the cores (using Pomeroy or ASC orientation device) and the hand samples. A declination correction was added to the azimuth measurements  (Korhonen et al., 2008;Brown et al., 2015). The enlarged map shows locations of archaeomagnetic sites that passed the acceptance criteria in Table 2 and listed in Table 3. Symbols and colors are as in Figure 3, where closed symbols show data from this study and open symbols show data from Segal (2003). Locations of sedimentary cores data shown in Figure 4 are shown with asterisks with the same color code as in Figure 4.
for all sites, except those of Segal (2003) from which we do not have the original measurement data. The latter results in a possible declination offset of the Segal (2003) dataset by 2-3 degrees.
The sampling, the formation, and the condition of the structures (sites) have a considerable effect on the uncertainty of the paleomagnetic directions. Man-made archaeological structures can collapse, break apart, incline or tilt, especially when the archaeological layer had been buried under a heavy overburden before being excavated. Therefore, extra care should be taken during paleomagnetic sampling of archaeological objects. To enable comparison between sites in the catalog, we assign to each site a "site formation quality index" (Qi) with fourscale grading indices ( Table 1). The highest score (Qi = 1) is granted when the entire periphery of an oven or kiln was sampled or when several bricks from at least two walls in a burnt structure were sampled. A medium score (Qi = 2) is given when only a segment of a wall, oven, or a furnace was sampled. A low score (Qi = 3) is given to samples with high orientation uncertainty, and Qi=4 indicates poorly oriented samples.

Paleomagnetic Measurements and Data Analysis
The sample set of Segal (2003) was measured in the paleomagnetic lab in the Geophysical Institute of Israel using a 2G cryogenic magnetometer, and in the paleomagnetic laboratory at the University of Rennes, France, using a spinner magnetometer and a Leti cryogenic magnetometer. Part of the sample set of Hassul (2015) was measured in the paleomagnetic laboratory at the Helmholtz Centre Potsdam, GFZ Germany using a 2G cryogenic magnetometer. The majority of the data were measured in the paleomagnetic laboratory at the Institute of Earth Sciences, the Hebrew University of Jerusalem using a 2G cryogenic magnetometer. Specimens from all sites underwent progressive demagnetizations with Alternating Field (AF).
All the raw measurement data, except the dataset of Segal (2003) were translated into the community standard MagIC format  and merged into a single measurement file. All the data with the exception of Segal (2003) were re-analyzed, including previously published data, using the Demag GUI program, which is part of the PmagPy software package . Paleomagnetic directions of specimens and site means were calculated using the principal component analysis technique (Kirschvink, 1980) and Fisher statistics (Fisher, 1953). The interpretations follow a fairly strict set of selection criteria listed in Table 2, accepting only specimens with MAD (Kirschvink, 1980) ≤ 5, DANG (Tauxe and Staudigel, 2004) ≤ 5, and sites with n (number of specimen per site) ≥ 8, k (Fisher, 1953) ≥60, and α 95 ≤ 6. The measurement data and the interpretations (except Segal's dataset) are available in the MagIC database (https://www2.earthref.org/ MagIC).

RESULTS AND DISCUSSION
The catalog consists of 76 paleomagnetic sites collected from 33 different locations (archaeological sites) in Israel. Forty seven sites passing the acceptance criteria in Table 2 are listed in Table 3, and their locations are shown in Figure 1B. The declinations and the inclinations of the accepted sites are shown in Figures 3A,B, where the inclination error bar is the α 95 value and the declination error bar is calculated using the equation: △ D = sin −1 sin α 95 cosI , where D is the declination error and I is the inclination. Sites that do not pass the criteria in Table 2 are listed in Table 4. In the catalog we distinguish between sites that have all their measurement data available in the MagIC database and can be downloaded and re-interpreted by others using different criteria (labeled "this study" in Tables 3, 4 and marked with filled symbols in Figures 1, 3) and sites of which we have only the site's mean parameters (labeled " Segal, 2003" in Tables 3, 4 and marked with open symbols in Figures 1, 3). From Figures 3A,B it can be seen that some periods have several coeval sites that show overlapping directions and demonstrate internal consistency and cross correlation between archaeological locations. These include: the fourteenth century BCE with two sites (Tel Megiddo and Tel Rehov); the thirteenth century BCE with two sites (Tel es-Safi and Tel Hazor); the beginning of the eighth century BCE with two sites (Tel Azekah, Tel Hazor,) the end of the eighth century with two sites (Bethsaida, Tel 'Eton), the third century BCE with two sites (Bnei Brak), and the second century BCE with two sites (Tel Shimron). An exceptional period with inconsistent directions is the tenth to ninth century BCE, which includes 11 sites from Tel Megiddo and Tel es-Safi, with non-overlapping directions. We interpret the latter as a time interval with fast changes in the geomagnetic field, and discuss this result in section Non-axial Dipole Field During the Iron Age Anomaly below.

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

Comparison With Sedimentary Data
There are several advantages of archaeomagnetism over sedimentary magnetic data. The archaeomagnetic TRM does not suffer from inclination shallowing, lock-in depth, and post-depositional effects associated with depositional remanent magnetization (DRM). Hence, paleomagnetic directions from archaeomagnetism can sometimes provide better precision and age control than sedimentary data. However, sedimentary magnetism provides continuous datasets spanning much larger time intervals than archaeomagnetism. Both types of records are available in Israel and we compare them in Figure 4 that shows data from four Holocene cores available in the GEOMAGIA50 database (Brown et al., 2015). Figures 4A-C show data from three piston cores raised from the Holocene Dead Sea (Frank et al., 2007a,b). These cores were obtained without azimuthal orientation and Frank et al. (2007a,b) corrected their declination profile by setting the mean declination value of the core top to zero. Their age model is based on a large number of radiocarbon measurements (Migowski et al., 2004). Given the uncertainty in the age models of the Dead Sea cores, there is fairly good agreement between the archaeomagnetic and the sedimentary data, especially for the past two millennia. The high inclination values observed in the ninth century BCE are not observed in the sedimentary data. This could be a result of inclination shallowing, post depositional magnetization, and smoothing of the sedimentary data. As the Holocene Dead Sea sediments are dominated by authigenic greigite (Ron et al., 2006;Frank et al., 2007b;Thomas et al., 2016;Ebert et al., 2018) a complicated magnetic acquisition mechanism is expected. Thus, differences between the archaeomagnetic and the sedimentary data are likely. Figure 4D shows data from a core taken in the Birket Ram maar lake in northern Israel, which was dated using only two radiocarbon ages (Frank et al., 2002). Here, the trends in the inclination and the declination profiles agree with the archaeomagnetic data, but the temporal resolution in Birket-Ram is much lower than at the Dead Sea. In summary, we observe a reasonable correlation between the sedimentary and the archaeomagnetic data highlighting the potential of combining these two types of records into a single joint master secular variation curve for the Levant. Yet, owing to the large uncertainties in both the sedimentary magnetic acquisition mechanism and the sedimentary age models this challenge requires a more detailed investigation. Shaar et al. (2016Shaar et al. ( , 2017 hypothesized a positive local geomagnetic anomaly in the Levant between 1050 BCE to 750 BCE which they termed "The Levantine Iron Age Anomaly" (LIAA). The paleointensity data from the Levant supporting the LIAA hypothesis ( Figure 5A) show high field values between the mid-Eleventh century BCE and the eighth century BCE and two geomagnetic spikes (virtual axial dipole moment, VADM > 160 ZAm 2 ) (Ben-Yosef et al., 2009Shaar et al., 2011Shaar et al., , 2016. Fourteen sites in Figure 3, from Tel Megiddo, Tel es-Safi, Tel Azekah, Tel Hazor, and Bethsaida cover the interval between 900 and 750 BCE. Of these sites, 12 show inclinations above 60 • , while two sites from Tel Megiddo show exceptionally high inclinations of 73 • (mgq05t2, Shaar et al., 2016) and 71 • (QTMB; Hassul, 2015; this study), considerably higher than the expected geocentric axial dipole (GAD) inclination in Jerusalem (51 • ). We note that one site in Segal (2003) dataset (Ceramic kiln from Kfar Menachem, see Supplementary Material) showed an even higher inclination of 81 • around this time. Yet, this site did not pass the selection criterion for the age as the age of the kiln was not supported by any direct dating method, only by a correlation with nearby archaeological sites. The declinations during the Levantine Iron Age Anomaly interval show large scatter in the ninth century with declinations ranging from −3 • to 23 • , and closely grouped values around 5 • in the eighth century. The angle between the archaeomagnetic directions and the axial dipole field ( Figure 3C) (Brown et al., 2015) following Frank et al. (2002Frank et al. ( , 2007a.

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

Archaeo-Chronological Applications
The forty seven archaeomagnetic data points reported here comprise a significant first step toward a master archaeomagnetic curve for the Levant, from which archaeomagnetic dating can be developed, using an approach similar to that described in Lanos (2004) and Pavon-Carrasco et al. (2014). Still, even in the absence Frontiers in Earth Science | www.frontiersin.org   with data of (Ben-Yosef et al., 2017)] are shown in black, and data from Syria are shown in gray (for references list see Shaar et al., 2016). (B) Inclination profile from this study, showing that the highest inclinations occur during the Levantine Iron Age anomaly, supporting the hypothesis of a local non-dipolar anomaly.
of a continuous secular variation curve, archaeo-chronological insights can be obtained by merely comparing archaeomagnetic data from structure whose age is not tightly constrained, with the available data. The burnt structure MGDF from Shahack-Gross et al. (2018) is an example of this approach. For dating the destruction that caused the fire at this archaeological level, Shahack-Gross et al. (2018) compared the paleomagnetic direction of MGDF with the available archaeomagnetic directions of periods with known destructions in Megiddo (Late Bronze III, Late Iron I, Iron IIA, and Iron IIB) and found that only one period fit to the archaeomagnetic data. Another promising potential of archaeomagnetism is the opportunity to use archaeomagnetic time-series of fast variation in the geomagnetic field, such as the Iron Age, to complement radiocarbon data, in particular for periods when the radiocarbon calibration curve is flat. For a more detailed review on addressing the Iron Age chronology problems with archaeomagnetism see Stillinger et al. (2016Stillinger et al. ( , 2018 and Herve and Lanos (2017).

SUMMARY AND CONCLUSIONS
We initiate here the first catalog of archaeomagnetic directions from Israel, with 76 sites that have been collected and analyzed since 2003. From this catalog, 47 sites pass a set of strict selection criteria and serve as a basis for future development of a master archaeomagnetic secular variation curve for the Levant. The measurement data of this catalog can be approached from MagIC database (https://www2.earthref. org/MagIC) and the archaeological information is provided in the Supplementary Material of this article. Some misfits with the global models prior to the second half of the first millennium BCE imply that more archaeomagnetic data are needed to refine our knowledge of fast secular variations during early periods in the Holocene. The most prominent feature in the new catalog is the high inclinations and high declination interval in the ninth century BCE. During this interval two sites, one published here, and another published in Shaar et al. suggesting an unusual field behavior. The large angular deviation occurred during a period with extremely high field intensity in the Levant, providing additional support to the Levantine Iron Age Anomaly hypothesis of Shaar et al. (2016Shaar et al. ( , 2017) of a regional high field anomaly at the beginning of the first millennium BCE.