First Full Vector Archeomagnetic Data From Northern Mexico

Several regional secular variation curves of the geomagnetic field have been proposed for Mexico over the last millennia. Despite a fairly large number of archeomagnetic data, these curves remain imprecise because of an uneven quality and geographic distribution, with a lack of data in Northern Mexico. Nine pottery kilns were sampled in Casas de Fuego, an archeological site in the Chihuahua state. These kilns belong to the Casas Grandes culture and were in use between 1250 and 1450 CE. Rock magnetic experiments indicate that the main magnetic phase is a Ti‐poor titanomagnetite in the SD range. Mean characteristic directions per kiln were estimated by alternating field and/or thermal demagnetization, and archeointensities with the Thellier‐Thellier protocol, after correction for anisotropy and cooling effects. The nine directions and eight intensities agree with data from USA and Mexico but do not support the peak in inclination modeled by the Mexican secular variation curve around 1200–1300 CE. For the last millennium, the Western North American curve is consistent with the SHAWQ2k global model, and better reflects the secular variation in Northern Mexico. For intensity, neither the SHAWQ2k model nor the regional intensity curve can depict the rapid secular variation that likely occurs circa 1500 CE. The Casas de Fuego results are the first full vector determinations obtained in Northern Mexico.

Only a few data are available outside these areas: lava flows from the Cascade Mountains (e.g., Hagstrum & Champion, 2002) and archeological baked clays from the Mayan area in Southern Mexico and Guatemala (e.g., Alva-Valdivia et al., 2010;Fanjat et al., 2013;Morales et al., 2009;Wolfman, 1990). No data have been presently published in large areas, such as Northern Mexico. It is interesting to note that this distribution is partly related to the location of the experimental laboratories and does not reflect the archeomagnetic potential of these areas. Also, only a few dozens of data, generally acquired on volcanic flows are full vector (e.g., Mahgoub, Juárez-Arriaga, Böhnel, Siebe, & Pavón-Carrasco, 2019). Another weakness of the Northern and Central American data set is the inhomogeneous data quality, which clearly impacts the recovery of the secular variation in these regions (e.g., Hervé, Perrin, Alva-Valdivia, Tchibinda Madingou, et al., 2019;Jones et al., 2020). A large proportion of intensity data do not fulfill modern quality criteria, such as the necessity to appropriately correct for TRM anisotropy and cooling rate effects, at least for archeological baked clays (Hervé , Perrin, Alva-Valdivia, Tchibinda Madingou, et al., 2019;Mahgoub, Juárez-Arriaga, Böhnel, Manzanilla, & Cyphers, 2019).
As a new step to improve the Central American archeomagnetic database, nine pottery kilns were sampled in the archeological site of Casas de Fuego in Chihuahua state, Northern Mexico. This site is part of the Mogollon cultural area, one of the most important in Northern Mexico and Southwestern USA.

Archeological Context and Sampling
Northern Mexico and Southwestern USA have a high potential for archeomagnetism through their rich archeological heritage. During the pre-Columbian period, these regions belonged to the so-called Oasisamerica (Kirchhoff, 1954). They were divided into three main cultural areas in the first half of the second millennium CE: Ancestral Puebloan (Four Corners area, before named Anasazi), Hohokam (SW Arizona), and Mogollon (SE Arizona, SW New Mexico, West Texas, and northern parts of Chihuahua and Sonora states in Mexico). These three cultural complexes have similar characteristics, for example large buildings in adobe divided into many rooms. The main archeological site of the Casas Grandes culture, belonging to the Mogollon area, is Paquimé in NW Chihuahua. This famous site, classified to the UNESCO World Heritage List, was occupied from circa 700 to 1500 CE with a florescence period at Middle period between 1150 and 1450 CE (Pailes, 2017). This period was characterized by intense economic activity with long-distance trade with other Oasisamerica and Central Mexico regions. Fine painted potteries were produced, such as the Figure 1. Location of Casas de Fuego area (red diamond) and available archeomagnetic data (directional data are in blue, and intensity data in pink). Circles (triangles) are for data coming from sites using archeological (volcanic) material.

of 13
Mimbres and Ramos traditions (Alva-Valdivia et al., 2021). The population also significantly increased, as evidenced by the high number of smaller contemporaneous archeological sites close to Paquimé.
Casas de Fuego is located ∼18 km south of Paquimé (Cruz-Antillon & Maxwell, 2017) in the Arroyo Seco valley. The archeological area was discovered in 1999 and excavated in 2000-2002(Cruz-Antillon, 2003. Several settlements are dispersed over a few km 2 in two main areas on both sides of the valley, Arroyo Seco A and B. Our study focuses on Arroyo Seco B, with two archeological sites 1.5 km away (CF1 and CF2,Figure 2a). Buildings in adobe are rather well preserved with size varying from 1 to 100 m 2 . The larger complex was excavated on archeological site CF1 by R. Cruz-Antillón (Figures 2b and 2c), with square (2 × 2 m) or rectangular (2 × 4 m) rooms. In comparison to other known archeological sites in the Paquimé area, the Arroyo Seco B buildings have the peculiarity of presenting some rooms that were clearly exposed to fire with well burnt walls and ashy floors (Figures 2c and 2d). Such rooms are interpreted as kilns for pottery baking (Cruz-Guzman & Nava Maldonado, 2008), as supported by the presence of overheated pottery sherds related to a failed baking process. Archeological artifacts, dendrochronology and radiocarbon converge to an occupation of the Casas de Fuego area between 1250 and 1450 CE, that is, at the end of Casas Grandes Middle period.
The archeomagnetic sampling was performed on the best-preserved parts of the burnt walls of CF1 and CF2 (Figure 2c). At CF1 archeological site, five rooms (CF1-1 to CF1-5) were sampled in the main building ( Figure 2b) and another room (CF1-6) in a close small building. At CF2 archeological site, two rooms (CF2-7 and CF2-8) were sampled. Finally, we studied another burnt room (six blocks) in a smaller archeological site called El Papalote (PAP). This site is located in the same valley about 8.5 km north of CF1. No direct dating could be performed at PAP but the similarity of the archeological context let us assume that it was used at the same period as CF1-CF2. A total of 48 paleomagnetic block samples, between 4 and 8 per room, were sampled ( Figure S1) and oriented with magnetic and sun compasses on horizontally leveled plaster caps. After induration with sodium silicate in the laboratory, all block samples were cut in 4-15 cubic specimens, 2 cm side.  Figure S1 for labels), burnt (unburnt) walls in orange (brown); and (d) picture of room CF1-1 with plaster caps. Maps modified from Cruz-Antillon (2003).

Rock Magnetism
The magnetic carriers and their thermal stability were estimated by different rock magnetic experiments. Susceptibility versus temperature (K/T) curves were done with an AGICO Multi-frequency Kappabridge (MFK1-CS3). Two successive heating-cooling cycles, up to 350 and 620°C, were performed on powder specimens of ∼200 mg. Curie temperatures (Tc) were estimated at the point where the curvature of the concave part of the heating curve is maximal. First-order reversal curves (FORC) experiments were measured with a regular grid and a 200 ms averaging time on selected samples. Data were processed by the FORCinel routine (Harrison & Feinberg, 2008) with VARIFORC smoothing parameters (Egli, 2013).
Most kiln walls present a change in color from reddish when the wall was in contact with the fire to blackish inside the wall (Figure 3). To test the influence of the color gradient on the magnetic mineralogy, detailed rock magnetic analysis was done in four different zones from red to black (zones a-d, Figure 3). In general, heating and cooling K/T curves are fairly reversible, evidencing no mineralogical transformation, except for an example (Figure 3b) that shows the effect of oxidation with lower values of the cooling branch. The main magnetic phase is a Ti-poor titanomagnetite with Curie temperature ranging from 520 to 550°C. FORC diagrams present mainly a SD behavior with uniaxial, randomly oriented grains and no interaction indicated by a central ridge over B u direction and in one case a negative distribution ( Figure 3d).

Directional Analysis
All magnetic moments were measured with a 2G Enterprises cryogenic magnetometer. Alternating field (AF) demagnetization was done with an AF device in line with the cryogenic magnetometer, from 4 to 90 mT (3-5 mT steps). Thermal demagnetization was performed in an MMTD80A oven with temperatures ranging from 100 to 585°C (20-50°C steps). Sixteen pilot sister-specimens were demagnetized with both methods. As they give similar results, the rest of the specimens were AF demagnetized. Results were interpreted with the Puffin Plot software (Lurcock & Wilson, 2012), and the mean directions per kiln were calculated by Fisher statistic with the PmagPy software (Tauxe et al., 2016).
Two-third of the specimens present a single component of magnetization (Figures 4a and 4c), while the others also have a small secondary component (Figures 4b and 4d) partly maybe of viscous origin. Characteristic remanent directions are listed per specimen (Table S1) and plotted per kiln ( Figure S2). The mean directions per kiln, calculated at the block sample and at the specimen levels (Table 1), are less than 1° apart and the distribution parameters are similar. Considering the homogeneity between the two calculations, the specimen level was preferred because the larger number of determinations allows a more reliable calculation of the 95% confidence circles (α 95 ). Most mean directions are well-grouped ( Figure 5), except CF2-7 that has the larger α 95 and for which the kappa decreased at the block sample level. For these reasons, CF2-7 should be considered with caution.

Archeointensity Experiments
Archeointensities were determined using the Thellier-Thellier method (Thellier & Thellier, 1959), where specimens are heated twice at each temperature in a laboratory field applied in opposite directions (+z and −z). Every two steps, a partial thermoremanent magnetization (pTRM) check is made with an additional heating along +z to monitor the absence of mineralogical alteration during the experience. Temperature steps, ranging from 100 to 560°C, were chosen according to the unblocking temperatures estimated during thermal demagnetizations. All thermal experiments were performed in an ASC oven and a laboratory field of 40-50 μT (see Table S2) was applied during both heating and cooling. Results were displayed with Thelli-erTool4.0 software (Leonhardt et al., 2004), and the archeointensities were deduced from the ratios between the NRM lost and the pTRM gained at each temperature step.
Seventy-six specimens were measured and representative NRM/TRM plots and orthogonal diagrams are presented in Figure 6. Roughly one-third of our specimens, mainly from the reddish part of the walls, present a linear NRM/TRM plot with a good reproducibility of pTRM checks (within 5%) and only one component of magnetization on the orthogonal diagram (Figure 6a). On these most reliable results, the archeointensities   were estimated between 100-150°C and 540-550°C. Another third of the specimens presents a slight secondary component that is erased at fairly low temperatures, usually by 300-400°C (Figure 6b). Reliable archeointensities can still be estimated at higher temperature (Table S2). The last third has NRM/TRM diagrams either with no linear segment over a large enough temperature range or with two slopes (Figure 6c) and no archeointensity was estimated.
For baked clays, anisotropy effects on TRM intensity have to be investigated (e.g., Aitken et al., 1981;Veitch et al., 1984), even though they are generally low for such burnt structures (e.g., Kovacheva et al., 2009). The TRM anisotropy tensor is defined through six successive heatings and coolings at 530°C or 540°C, with the same applied laboratory field aligned along the specimen directions +x, −x, +y, −y, +z, and −z, followed in most cases by a stability check. Archeointensities were corrected for anisotropy effects only when the factor of evolution (%Evol ATRM ) was within the range 0%-5% (Table S2). Only 15 specimens could be corrected for the anisotropy effect. However, the corrected intensities are within 10% of the uncorrected values without systematic trend that gives us confidence in the reliability of the other intensity estimations.
TRM intensity also depends on the duration of cooling (e.g., Fox & Aitken, 1980) and the 1-h cooling used routinely during the archeointensity experiments is likely faster than the archeological cooling rate. The procedure of Gómez-Paccard et al. (2006) was carried out to correct this effect with a slow cooling (5 h) at the same temperature as used for the anisotropy correction, followed by a stability check. The duration of the slow cooling rate was constrained by technical reasons and may still underestimate the archeological cooling. However, it was shown that such an imprecision does not significantly disturb the average archeointensity values (Hervé, Chauvin, et al., 2019). Cooling rate correction factors (CR factor, Table S2) averaged around a mean value of 5% that was used to correct the intensities when the cooling rate evolution factor was too large (%Evol CR ≤5%, Table S2). Averaged archeointensity results per kiln (Table 1) show a slight decrease of the standard deviation after corrections.

Discussion
Altogether nine mean directions and eight mean intensities were obtained in this study. The differences, especially in declination, suggest that the activity of the kilns did not stop at the same time. These new data are the first full vector data obtained on archeological material in Mexico and Western North America. We compared them to intensities ( Data were relocated to Casas de Fuego via the Virtual Geomagnetic Pole (VGP) and the Virtual Axial Dipole Moment (VADM) assumptions for directions and intensities respectively. The non-dipolar component of the geomagnetic field introduces a relocation error that could reach a few degrees or a few μT at the maximum (Casas & Incoronato, 2007;Goguitchaichvili et al., 2020). We also compared the predictions, at Casas de Fuego, of global models that have been recently constructed with archeological and volcanic data (e.g., ARCH10k.1, Constable et al., 2016;COV-ARCH, Hellio & Gillet, 2018;BIGMUDI4k.1, Arneitz et al., 2019;SHAWQ2k, Campuzano et al., 2019). The predictions of these models being reasonably consistent, we choose to present on Figures 7 and 8 only the most recent SHAWQ2k model that has the advantage of weighting data according to their quality.

Comparison With Intensity Data and Models
For the last 2 kyrs, 192 intensity averages were found within 2,000 km from Casas de Fuego: 141 from Mexico (see Alva-Valdivia et al., 2021;Garcia et al., 2021;Hervé, Perrin, Alva-Valdivia, Tchibinda Madingou, et al., 2019;Mahgoub, Juárez-Arriaga, Böhnel, Manzanilla, & Cyphers, 2019) and 51 from Western North America (e.g., Bowles et al., 2002;Champion, 1980;Jones et al., 2020;Sternberg, 1982). Archeological baked clays constitute 80% of the data set. Consistently with Jones et al. (2020), we corrected the Southwestern USA mean values by a cooling rate effect of 6.73%, experimentally determined in their study. 32 volcanic data are available for Mexico and one for Utah. However, it can be noted that some Mexican flows have been studied several times, especially the Xitle lava sequence (Alva-Valdivia et al., 2020) that has been intensively sampled in Mexico City but also the Jorullo (Alva-Valdivia et al., 2019) and Paricutin lava flows .
The whole data set is plotted in Figure 7a, together with the new intensities that are within the wide range of VADM values (70-150 ZAm 2 ) observed in the first half of the second millennium CE. As several studies (Hervé, Perrin, Alva-Valdivia, Tchibinda Madingou, et al., 2019;Jones et al., 2020;Mahgoub, Juárez-Arriaga, Böhnel, Manzanilla, & Cyphers, 2019) showed the large role of the uneven data quality in this dispersion, we applied the following selection criteria: (a) only data acquired with the original or derived Thellier-Thellier and microwave methods, (b) at least three specimens, (c) a standard deviation (sd) lower than 10%, and (d) an age uncertainty less than ±200 years. Furthermore, paleointensity values on ceramics must have been corrected for anisotropy effects using the TRM tensor when the laboratory field was not applied parallel to the NRM. Data corrected for TRM anisotropy by the Mean (XYZ) correction were discarded as a bias up to 10-15 μT is possible (Hervé, Perrin, Alva-Valdivia, Tchibinda Madingou, et al., 2019). Finally, the cooling rate correction must also be considered for archeological baked clays, either by a coefficient directly estimated on the specimen or by applying an average rate. This set of criteria corresponds to the C1 set of Hervé, Perrin, Alva-Valdivia, Tchibinda Madingou, et al. (2019) to which we add one microwave determination. It is also closely similar to the one used by Mahgoub, Juárez-Arriaga, Böhnel, Manzanilla, and Cyphers (2019) to construct their secular variation curve for Mexico that we consider most reliable than the curve proposed by Garcia et al. (2021), not as selective for anisotropy correction.  Circles and triangles correspond to archeological baked clays and lava flows respectively. All data are relocated in Casas de Fuego via VADM. In (b), the red curve is the relocated Mexican curve of Mahgoub, Juárez-Arriaga, Böhnel, Manzanilla, and Cyphers (2019) and the black curve is the prediction at the archeological site of SHAWQ2k global model (Campuzano et al., 2019). Both curves are plotted with their 95% error envelope. Figure 8. Comparison of directions from Casas de Fuego (red stars) with data within 2,000 km from Western North America (a and b) and Mexico (c and d).
All compiled data are plotted in (a-c), whereas (b-d) presents the selected data set. Circles and triangles correspond to archeological baked clays and lava flows respectively. In (a), the blue data come from Southwestern USA and the green one from the Lower Mississippi river region. In (c), the brown data were acquired on burnt material carrying a TRM and the beige data on unburnt lime-plasters. In (b-d), the red curves with their 95% error envelope are the regional curves of Hagstrum and Blinman (2010) for Western North America and of Soler-Arechalde et al. (2019) for Mexico. The prediction at the archeological site of SHAWQ2k global model (Campuzano et al., 2019) is plotted in black with its 95% error envelope. All data and curves are relocated to Casas de Fuego via VGP. Böhnel, Siebe, & Pavón-Carrasco, 2019;Morales et al., 2006;Pérez Rodriguez et al., 2019) and 39 baked clays (Alva-Valdivia et al., 2021;Hervé, Perrin, Alva-Valdivia, Tchibinda Madingou, et al., 2019;Mahgoub, Juárez-Arriaga, Böhnel, Manzanilla, & Cyphers, 2019). Other high-quality archeomagnetic data exist in Mexico, but they are either beyond the 2,000 km limit (e.g., Fanjat et al., 2013) or older than 2 kyrs (e.g., Duran et al., 2010;Hervé, Perrin, Alva-Valdivia, Rodríguez-Trejo, et al., 2019). For Western North America, 47% of the data set passed the selection criteria, with one lava flow (Champion, 1980) and 23 data on ceramics (Bowles et al., 2002;Jones et al., 2020;Sternberg, 1982). As well in Mexico as in Western North America, the most discriminant criterion was the anisotropy correction, usually not performed in early studies (e.g., Bucha et al., 1970;Lee, 1975).
While the selection discards the highest and lowest intensity values (Figure 7b), the Mexican and Western North American data sets remain consistent, even though data are now missing for some periods as around 500-600, 1200, and 1500-1700 CE. Between 1200 and 1500 CE, besides Casas de Fuego, five Western North American and six Mexican data are available with intensities ranging from 80 to 130 ZAm 2 . The closest results from our sites are the Paquimé results (pink circles in Figure 7; Alva Valdivia et al., 2021) that are in better agreement with the Mexican secular variation curve (Mahgoub, Juárez-Arriaga, Böhnel, Manzanilla, & Cyphers, 2019), while our results are more consistent with the SHAWQ2k global model (Campuzano et al., 2019). The difference between regional and global models can be explained by the absence of high intensities in the selected data set of Mahgoub, Juárez-Arriaga, Böhnel, Manzanilla, and Cyphers (2019), but it was also shown that the global models are strongly influenced by low-quality data (Hervé, Perrin, Alva-Valdivia, Tchibinda Madingou, et al., 2019) that tend to smooth secular variation. We think that none of the models properly described the likely rapid secular variation at that time. New reliable data, especially in Northern Mexico, will be necessary to confirm this hypothesis.

Comparison With Directional Data and Models
Many directions (589) were published from sites within 2,000 km of Casas de Fuego for the last 2 kyrs, a large majority coming from archeological fireplaces or ovens. Four-hundred two archeological data are from Southwestern USA (Arizona, Colorado, New Mexico, Utah), 61 archeological data are from the Lower Mississippi river (Arkansas, Mississippi, Tennessee, and Missouri), and 10 volcanic data are from North California and Utah (e.g., Hagstrum & Champion, 2002). The Western North American data set appears homogeneous with a rather uniform temporal distribution over the last two millennia (Figure 8a).
The Mexican data set contains 116 data (without Casas de Fuego): 100 on archeological material and 16 on lava flows (see Soler-Arechalde et al., 2019). Most of them are located in Central Mexico along the Trans Mexican Volcanic Belt and in the Maya area. There are no directions from North of 21°N, emphasizing even further the interest of our work for the Mexican database. The Mexican data set is less homogeneous than in Western North America, with a wider range of inclination values and some data with declinations up to 40°W (Figure 8c). To test the geomagnetic origin of this dispersion, we also performed a critical analysis of the directional data with (a) at least five specimens; (b) an α 95 lower than 5°; and (c) an age uncertainty less than ±200 years. A peculiarity of the Mexican data set is the presence of lime-plasters covering house floors, which mix water, lime, and lithic clasts (Hueda-Tanabe et al., 2004). Due to the volcanic environment of Central Mexico, the lithic clasts are rich in ferromagnetic particles, and a remanent magnetization is acquired by the alignment of the particles parallel to the local geomagnetic field when the plaster is mixed, therefore carrying a detrital remanent magnetization that does not record the ancient geomagnetic field as accurately as a TRM. Furthermore, lime-plasters data have systematically higher uncertainties (Figure 8c), so we choose to consider only directions from burnt structures. Finally, we discarded the direction of Tula415 (Figure 8c) that passes the criteria but was considered as unreliable in the original article of Wolfman (1990).
The selection process preferentially removed data with the westernmost declinations and highest or lowest inclinations, and the remaining 49 Mexican directions (39 and 10 from archeological and volcanic material, respectively) form a more homogeneous data set ( Figure 8d). Unfortunately, the selection also removes almost all data that could be compared to Casas de Fuego. The Casas de Fuego data are in good agreement with the Western North American data set, with selection (420 data, Figure 8b) or without selection (473 data, Figure 8a).
Casas de Fuego data are basically in agreement with the predictions of the SHAWQ2k global model and with the regional secular variation curves of Western North America (Hagstrum & Blinman, 2010) and Mexico (Soler-Arechalde et al., 2019), relocated at the archeological site (Figures 8b-8d). The Western North America regional and global curves show a plateau of the inclination over the last millennium, while the curve of Soler-Arechalde et al. (2019) exhibits a rapid and short increase of the inclination between 1200 and 1400 CE. Our CF2-7 kiln result could seem to support this feature but, as said before, it should be considered with care. As can be seen from Figure 8d, all data with very high inclinations have been discarded from our selection, therefore this peak seems to result from the inclusion of low-quality data (e.g., lime-plaster data and Tula415) in the calculation of the mean directions by the moving windows technique, and likely does not reflect a geomagnetic feature. The Mexican secular variation curve over the last millennium needs to be re-established with new reliable data. For the moment, we suggest to use the Western North America directional curve (Hagstrum & Blinman, 2010) for archeomagnetic dating in Northern Mexico. If applied to Casas de Fuego, the fairly steady increase in declination between 1100 and 1800 CE may indicate a relative chronology between our kilns, with CF1-1 being the earliest in activity and CF1-6 the latest.

Conclusions
Casas de Fuego results are the first full vector determinations obtained in Northern Mexico. The difference between kilns, especially in declination, suggests that the kilns were used at different times. The nine directions and the eight intensities are consistent as well with selected data from Western North America than from Mexico. Considering the directional secular variation, the SHAWQ2k global model and the regional curves for Western North America (Hagstrum & Blinman, 2010) and Mexico (Soler-Arechalde et al., 2019) show the same tendency over the last two millennia if the last part of the Mexican curve is not considered. For the last millennium, the regional Western North America curve seems more suitable for archeomagnetic dating in Northern Mexico. For the intensity secular variation, SHAWQ2k model gives predictions higher than Mahgoub, Juárez-Arriaga, Böhnel, Manzanilla, and Cyphers (2019) intensity curve but none seems to reflect the rapid secular variation that seems to occur between 1200 and 1500 CE.