Abstract
Four monogenetic Holocene lava flows located within the Michoacán-Guanajuato Volcanic Field, Mexico, were sampled for paleomagnetic dating. These flows (namely, El Infiernillo, Malpaís Las Víboras, El Capaxtiro, and Malpaís Prieto) are within the heartland of the pre-Hispanic Tarascan civilization and were inhabited repeatedly since at least 100 BC, but no relation with the volcanic evolution has been proposed so far. The stratigraphically oldest lava flow, El Infiernillo, has a radiocarbon age of 3200 ± 30 yr BP (cal. 1525–1420 BC), and it was used to validate the method. Using full-vector paleomagnetic data from three sites as input for paleomagnetic dating applying the global paleosecular variation model SHA.DIF.14k, an age range of 1500–1370 BC was obtained. Two age ranges of 1340–1230 and 1030–940 BC were obtained for Malpaís Las Víboras. A younger age range of 200–80 BC was obtained for the Capaxtiro lava flow and, finally, the Malpaís Prieto lava flow erupted within the range of AD 830–960. The human occupation history of these flows started around 100 BC during the late Pre-Classic, probably shortly after the Capaxtiro eruption. Archeological records indicate an abandonment of the entire area around AD 900 (late to terminal Classic), which coincides with the paleomagnetic age of the Malpaís Prieto eruption. Interestingly, this area was heavily repopulated again only few hundred years later around AD 1250 and belongs to the core region in which the Tarascan civilization has its roots. The eruption recurrence interval of roughly 1000 years indicates that a new monogenetic eruption should be expected to occur again in the future and that this area deserves to be studied in greater detail with particular emphasis on the impact of past eruptions. This could help to better evaluate volcanic hazards and design preparedness strategies to minimize the impact of a future eruption.
Introduction
Paleomagnetic secular variation (PSV) refers to temporal and spatial variations of the Earth’s magnetic field (EMF) due to internal processes (flow of liquid in the Earth’s outer core) in geologically short timescales (≪1 Ma). In recent decades, numerical models describing the observed temporal evolution of the geomagnetic field during the past centuries have been developed, for example, GUFM (Jackson et al., 2000) and UFM1 (Bloxham and Jackson, 1992). For longer periods of time, the evolution of the EMF can only be derived from the record in different archives (lake sediments, heated archeological artifacts, and lavas). Data obtained from archeological and volcanic products are the preferred inputs for the regional and global modeling of the EMF variations (Korte and Constable, 2003, 2005, 2011; Korte et al., 2009; Pavón-Carrasco et al., 2009, 2014), as these are based on high fidelity thermoremanent magnetization (TRM). Obtaining high-resolution PSV records in a particular region is of interest because they can be applied in paleomagnetic dating (e.g. Arrighi et al., 2006; Böhnel et al., 2016; Hagstrum and Blinman, 2010; Roperch et al., 2015; Speranza et al., 2006, 2008; Tanguy et al., 2003), among other fields. Paleomagnetic age constraints can be helpful in solving various types of scientific problems. For example, by dating lava flows, it can contribute to assess hazards in active volcanic regions (Cassidy, 2006). But it can also help in dating archeological remains, especially in areas where written historic records are absent. Both of these themes are addressed in this case study in which the Malpaís de Zacapu Holocene lava flows (namely, El Infiernillo, Malpaís Las Víboras, El Capaxtiro, and Malpaís Prieto), located at the W margin of the Zacapu lacustrine basin in the Michoacán-Guanajuato Volcanic Field (MGVF) in central-western Mexico were targeted. This field (Figure 1) has the largest concentration of monogenetic volcanoes in the Trans-Mexican Volcanic Belt (TMVB; Siebe et al., 2014) with more than 1000 monogenetic vents. Of these, the majority are scoria cones including the historical Jorullo (1759–1774) and Paricutin (1943–1952) volcanoes. In addition, ~400 medium-sized volcanoes (shields and domes), fewer viscous flows and lava domes, and rare maars occur (Chevrel et al., 2016a, 2016b; Guilbaud et al., 2011, 2012; Hasenaka and Carmichael, 1985, 1987; Kshirsagar et al., 2015, 2016). Notably, pre-Hispanic humans have continuously inhabited the territory of the present State of Michoacán since at least 5000 yr BC (Beekman, 2010; Faugère, 2006; Watts and Bradbury, 1982). Earliest sedentary populations around 2000 BC were succeeded by small-village societies, and from AD 1350 onward, the Tarascan Empire became established around the shores of Lake Pátzcuaro (Pollard, 1993). Considering the abundance of young volcanoes in this region, the impact of volcanism on pre-Hispanic human development and population migrations should not be ruled out. Unfortunately, the absence of written sources from archeological sites complicates evaluating the impact of volcanic eruptions on the pre-Hispanic populations of Michoacán (Pereira et al., 2013), as it can be established only based on archeological findings and their interpretation.

Digital elevation model (modified from Kshirsagar et al., 2015) of the MGVF showing the location of the Malpaís de Zacapu lava flows at the western margin of the Zacapu lacustrine basin. The major fault system is Cuitzeo Fault System (CFS), the most prominent fault system affecting the Zacapu basin. Black rectangle indicates area covered by the geological map shown in Figure 2. Inset map at lower right corner shows location of the MGVF within the TMVB.
Geological background: The MGVF and the Zacapu basin
The Holocene Malpaís de Zacapu lava flows (malpaís means ‘badland’ in Spanish) are located in the MGVF that forms the western-central segment of the TMVB (Figure 1). This continental arc stretches across central Mexico for 1200 km in an East-West direction from the coast at the Gulf of Mexico to the Pacific Ocean. It traverses the Mexican Altiplano, a highland characterized by active normal faulting and horst-and-graben structures that result in the formation of basins often occupied by broad (but shallow) lakes, such as the Pátzcuaro and Cuitzeo lakes in Michoacán.
The TMVB is related to the subduction of the oceanic Cocos Plate underneath the continental North American Plate (e.g. Blatter and Hammersley, 2010; Gómez-Tuena et al., 2003; Kim et al., 2012; Pardo and Suárez, 1995) and consists of a large number of late Tertiary to Quaternary maars, scoria cones, domes, calderas, and strato-volcanoes, the chemical and mineralogical composition of which is largely calc-alkaline (e.g. Carmichael, 2002; Demant, 1978; Ferrari et al., 2012). One notable feature of the TMVB is the abundance (>3000) of scoria cones and other types of small monogenetic volcanoes, which outnumber by several orders of magnitude the few dozens of the much larger strato-volcanoes (Guilbaud et al., 2012; Hasenaka and Carmichael, 1985).
Within the boundaries of the MGVF, which occupies an area of ~40,000 km2, the surface topography is dominated by andesitic Plio-Quaternary volcanic landforms. Intense recent volcanic activity conceals most of the older rocks; hence, the underlying basement can only be inferred from outcrops located beyond the limits of the MGVF, which must include plutonic rocks as evidenced by partially fused granodiorite xenoliths found in products of the Paricutin (McBirney et al., 1987; Wilcox, 1954) and Arocutin (Corona-Chávez et al., 2006) scoria cones.
The Holocene Malpaís de Zacapu lava flows are located at the W margin of the ENE-WSW trending Zacapu tectonic basin (Figure 1). The lowest part of the basin (1980 m a.s.l.) is today occupied by a cultivated flat surface of lacustrine origin that is surrounded by Plio-Quaternary volcanoes that are mostly basaltic andesite to andesite in composition (Demant, 1992; Siebe et al., 2013). A few dacites and rhyolites occur forming domes and ignimbrite sheets. The area is characterized by active normal domino-style fault systems that form the western extent of the seismically active Cuitzeo Fault Zone (Johnson and Harrison, 1990) also called the Morelia-Acambay Fault System (Garduño-Monroy et al., 2009; Suter et al., 2001). Faulting around the basin follows a general direction of N65°E to N85°E with dip-directions toward the NNW and SSE. Noticeably, faulting is less evident toward the west of the basin, because of a denser coverage by younger (<40 ka) volcanics (Pasquaré et al., 1991).
Hasenaka and Carmichael (1985) and Ban et al. (1992) reported first radiometric ages (14C and K-Ar, respectively) for the MGVF, including a few for volcanoes in the Zacapu basin. New radiometric dates of volcanoes in the eastern and northern parts of the basin (14C and 40Ar-39Ar) were published by Kshirsagar et al. (2015, 2016).
The only previous work focusing on the Malpaís de Zacapu lava flows was undertaken by Demant (1992) and includes a geologic map with stratigraphic data, as well as chemical and mineralogical analyses of the main lava flows. In his seminal work, the author established that the flows comprising the malpaís are andesitic and must be Holocene in age due to their youthful morphology. They were issued from four different eruption sites in the following chronological order: El Infiernillo, Malpaís Las Víboras, El Capaxtiro, and Malpaís Prieto. Our study confirms in essence these findings and was primarily aimed at dating the lava flows. Fieldwork aided by digital elevation model–based morphological observations and geochronological data were integrated to construct a new geological map (Figure 2), which shows the overlapping small-volume Holocene monogenetic lavas underlain by Pleistocene volcanic units. The oldest unit (Mesa El Pinal) is cut by a prominent normal fault and probably not much older than 2 Ma (early Pleistocene).

Geological map showing sampling locations of the young lava flows dated in this study.
The oldest of the four Holocene flows, El Infiernillo, was issued from the Las Vigas volcano, a small (~100 m high) scoria cone. Its olivine-bearing basaltic andesite (SiO2 = 56.3–58.7 wt.%) lavas are ~50 m thick and cover a minimum surface of 5.8 km2 with a volume estimated at 0.3 km3. Mineralogically, they are composed of forsteritic olivine (~1 mm), pargasitic hornblende (frequently opacitisized), augite, and hypersthene phenocrysts in a glassy matrix with feldspar microlites and opaque oxides. One paleosol sample (ZAC-1524A, laboratory no. B-411354, latitude 19°51′36.6″, longitude 101°51′08.9″, altitude 2192 m a.s.l.) directly underlying a 45-cm-thick deposit of coarse ash to fine lapilli fallout layers from the Las Vigas scoria cone (Figure 3) was dated by the accelerator mass spectrometry (AMS) method at Beta Analytic (Miami, Florida) at 3200 ± 30 yr BP (∂13C = −22.2). The date was calibrated to calendar years by applying the Stuiver and Reimer (1993) procedure and with the help of the CALIB Computer Program (version 7.1, n.d.; IntCal13 calibration curve), which yielded a 95% probability range (2σ) of cal. 1525–1420 BC. This is the only radiocarbon sample that could be found during the course of several field campaigns and its dating turned out to be crucial for carrying out with confidence the paleomagnetic dating of the remaining younger lava flows (see below).

Photographs showing the outcrop where paleosol sample ZAC-1524 dated at 3200 ± 30 yr BP was obtained: (a) Roadcut 600 m NW of Las Vigas scoria cone (shown in background). (b) Close-up view of the Las Vigas coarse ash fallout with sampled ochre-brown paleosol directly underneath. Spatula for scale is 25 cm long.
To the west of the Las Vigas cone and stratigraphically above is the Malpaís Las Víboras lava flow. This flow is also quite thick (~70 m) with steep margins attesting to a high viscosity during emplacement. It covers an area of 5.9 km2 with a volume of ~0.4 km3 and is andesitic in composition (SiO2 = 61.1–62.5 wt.%). Observation under the polarizing microscope revealed plagioclase, pargasitic hornblende (mostly opacite), augite, and hypersthene phenocrysts in a glassy matrix with abundant feldspar microlites and opaque oxides.
The most voluminous (3.2 km3 covering an area of 20.9 km2) of the eruptions corresponds to El Capaxtiro whose multiple ~150-m-thick overlapping lava flows were issued from a small cone (El Capaxtiro proper). They reached as far as 6 km to the east toward lake Zacapu, where they directly cover lacustrine deposits, as observable in the immediate vicinity of the city of Zacapu. Capaxtiro lavas are the most silicic (SiO2 = 61.1–64.2 wt.%) and also cover the southern portion of the Infiernillo lavas. Petrographically, they include plagioclase, pargasitic hornblende (mostly opacite ‘ghosts’), augite, and hypersthene phenocrysts in a glassy matrix with abundant feldspar microlites and opaque oxides.
Finally, Malpaís Prieto is not only the youngest but also the smallest of the studied lava flows. It is quite thick (~90 m) and covers an area of 5.7 km2 with a volume of ~0.5 km3. As its name implies (prieto means dark colored in Spanish), its rough blocky surface is pitch black and resembles a moonscape, being almost devoid of vegetation. Just from its appearance, it can be judged that it must be one of the youngest lava flows in the entirety of Michoacán. Its composition (SiO2 = 61.5–62.8 wt.%) is quite similar to Malpaís Las Víboras and El Capaxtiro and also contains plagioclase, pargasitic hornblende (mostly opacite), augite, and hypersthene phenocrysts in a glassy matrix with abundant feldspar microlites and opaque oxides.
Finally, it is worth mentioning that rounded quartz grains (1 mm in size) with reaction coronas of augite, as well as 1- to 2-mm partially resorbed plagioclase crystals displaying polysynthetic twinning are ubiquitous in all lavas. They are in disequilibrium and xenocrystic in origin. We suspect that they either stem from shallow Tertiary plutons or from younger shallow partly crystallized silicic magma chambers from which they were incorporated by ascending mafic magmas as attested by the forsteritic olivines found in the Infiernillo lavas. All of these observations point toward a complex petrologic origin of these monogenetic magmas, as indicated in various recent studies elsewhere in the MGVF (e.g. Rasoazanamparany et al., 2016).
Environmental and archeological background
The lake basins of Michoacán have been of interest for their archeological records. It seems that favorable aquatic and riparian conditions attracted early nomadic humans and promoted the development of agriculture in this region, which eventually became a major dwelling hub for the Pre-Hispanic Tarascan populations. Also known as Purépechas, the Tarascans erected around Lake Pátzcuaro one of the largest Mesoamerican empires during the late Post-Classic period, prior to the conquest by the Spaniards. The Tarascans and their predecessors resided mostly near riverine valleys and lakes (O’Hara et al., 1993; Pollard, 1993, 2012), including the area of the Zacapu basin. In this context, several archeological excavations (e.g. Arnauld et al., 1994; Arnauld and Faugére-Kalfon, 1998; Michelet and Carot, 1998; Michelet et al., 2005; Pereira, 2005) have been undertaken within the basin. At the same time, and in order to define the environmental factors that fostered the rise of human civilization in this region, paleoclimate studies focusing on the analysis of the lake sediments (and particularly on their microfossil contents) have been carried out (e.g. Metcalfe, 1992, 1995; Newton et al., 2005; Pétrequin, 1994; Telford et al., 2004; Xelhuantzi-López, 1994). Although particular attention has been paid to the Holocene, some of these records go back as far as 52,000 yr BP (e.g. Correa-Metrio et al., 2012; Metcalfe and Harrison, 1984; Ortega-Guerrero et al., 2002; Tricart, 1992). These studies indicate that the extent and depth of the lake Zacapu have changed over time and that conditions even turned marshy during the late Pleistocene (Ortega-Guerrero et al., 2002; Tricart, 1992). The lake possessed a natural shallow discharge toward the North (Siebe et al., 2012) until it was artificially drained in the late 19th century to gain fertile land for agricultural purposes (Noriega and Noriega, 1923).
Archeological investigations carried out in the Zacapu area since the early 1980s within the frame of the Michoacán Project (Michelet, 1992; Michelet et al., 1989) have documented a human occupation sequence that starts around 100 BC and continues until the conquest by the Spaniards. Figure 4 shows the ages of ceramics found in 18 excavation sites, which indicates continuous occupation between AD 550 and about AD 1550, with a clear interruption between AD 900 and 1250. These studies as well as more recent investigations that initiated in 2010 under the umbrella of the Uacusecha arqueological project show that the Malpaís de Zacapu area was densely populated during the middle/late Pre-Classic (AD 1250–1450). Surface prospection as well as the acquisition of a high-resolution light detection and ranging (LIDAR) image indicates that the Infiernillo, Capaxtiro, and Malpaís Prieto lava flows were occupied by four large urban settlements. The remains of thousands of domestic units as well as 50 ceremonial complexes with pyramidal mounds are still observable (Forest, 2014; Michelet, 1998; Migeon, 1998; Pereira et al., in press) These recent arqueological studies also show important demographic fluctuations, which probably correspond to migration events (Arnauld et al., 1993; Arnauld and Faugére-Kalfon, 1998; Pereira et al., 2013). As discussed further on it is possible that some of these ruptures in the history of the pre-Hispanic settlements could be related to episodes of volcanic activity in this area.

Time-graph showing pre-Hispanic human occupation as revealed by systematic excavations at different sites around the Malpaís de Zacapu area (modified after Michelet, 1992, with new data obtained during the ongoing Uacusecha archeological project). Site 31 is located on the Malpaís Prieto. Note the sudden abandonment of several sites around AD 900, followed by a hiatus that lasted 350 years until AD 1250.
Sampling procedures
The four lava flows were sampled at 21 different sites in order to check the within-flow paleomagnetic consistency (e.g. Böhnel et al., 2016; Hagstrum and Champion, 1994; Speranza et al., 2006). Sites are listed in Table 1, and their distribution is shown in Figure 2. No road cuts exhibiting the interior of the flows were found, and field observations indicate that their thickness ranges from 10 to almost 200 m with blocky surfaces and prominent marginal levées. Finding suitable drilling sites where the lava was in situ since it cooled down was, therefore, not an easy task. In order to avoid erroneous site mean directions due to moved blocks, multiple sites distributed as far as possible from each other were sampled from each flow. Based on the above, sampling was carried out with a gasoline-powered drill with a 25-mm diamond drill bit. Cores were between ~8 and 15 cm long and oriented using both magnetic and sun compasses to check for local magnetic anomalies (e.g. blocks struck by lightning). The difference between the magnetic and sun readings did not exceed 9° with an average of ~4°. From the four flows, a total of 141 cores (21 sites) were sampled representing an average of 7 cores per site with a minimum and maximum of 5 and 14, respectively.
Site mean paleomagnetic directions: latitude and longitude of the sampling coordinates.
n: number of samples used in the calculation of the site mean direction; N: total number of samples measured; k: precision parameter; α95: 95% confidence level; Dec: declination; Inc: inclination.
Shaded gray rows represent rejected sites or sites where no meaningful site mean could be calculated (for interpretation, see text). Results on site level are shaded in light blue color.
Laboratory experiments
In the laboratory, it is customary to cut each drill core into several specimens 22 mm in length, providing at least three specimens. In order to obtain reliable mean paleomagnetic direction measurements of the lava flows, one specimen from every drill core was utilized. Low-field volume susceptibility (k) was measured for each specimen using a susceptibility meter (MS2B, Bartington Instruments Ltd), and remanence magnetization measurements were made with an AGICO JR-5 spinner magnetometer (noise level ~5 × 10−6 A m−1). The direction of the characteristic remanence directions (ChRM) was determined by means of stepwise progressive alternating field (AF) demagnetization. For AF demagnetization, an AGICO LDA-3 equipment was used, and the samples were demagnetized in 14 steps: 3, 5, 7, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, and 90 mT. Demagnetization results were analyzed with the aid of the program PMGSC 4.2 (Enkin, 2005). For statistical analysis and graphical representation of the data, the program PMag Tool 4.2b by Mark W. Hounslow was used. Characteristic remanence directions were calculated applying principal component analysis (PCA) as described by Kirschvink (1980). The site mean directions were determined using statistical procedures delineated by Fisher (1953). Thermomagnetic curves (Tmax = 650°C) were done on two cores per flow in order to define their magnetic mineralogy configuration and thermal stability. For paleointensity (PI) determinations, we carried out double-heating Thellier experiments (Thellier and Thellier, 1959), using the IZZI protocol (Tauxe and Staudigel, 2004), which combines the in-field/zero-field (Aitkin et al., 1988) and zero-field/in-field (Coe, 1967) protocols. Both partial TRM (pTRM) checks (Coe, 1967) and pTRM tail checks (Riisager and Riisager, 2001) were executed in order to analyze the reliability of the forthcoming results. PI experiments were conducted with an ASC Scientific TD48 furnace. Temperature steps used were 100°C, 200°C, 250°C, 300°C, 340°C, 370°C, 400°C, 430°C, 460°C, 490°C, 510°C, 530°C, 560°C, and 580°C. pTRM checks were done at 100°C, 250°C, 340°C, 400°C, 460°C, and 510°C, while pTRM tail checks were done at 250°C, 340°C, 400°C, 460°C, 510°C, and 560°C. Results of the IZZI-Thellier data were analyzed using the ThellierTool4.11 software (Leonhardt et al., 2004). In order to judge the credibility of our PI estimates, acceptance criteria sets A and B as given in the Thellier tool (Leonhardt et al., 2004) with the modifications of Paterson et al. (2014) were used in this study.
Results
Rock-magnetic analyses
Figure 5 shows the variation of the natural remanent magnetization (NRM) intensity and magnetic susceptibility for the 21 studied sites, together with lines of constant Königsberger’s factor Q. From the three sites of El Infiernillo, we observe that both Inf1 and Inf2 have similar NRM intensities and susceptibilities (Figure 5a). Together they have a geometrical average NRM intensity of 8.80 A m−1 (excluding one specimen with NRM < 1 A m−1) and an average susceptibility of 24.4 × 10−3 (SI). For Inf3, we observe that the susceptibility is nearly half of the aforementioned sites but has approximately the same NRM intensities. Most of the Infiernillo specimens have Q values between 6 and 28 (Figure 5a). In contrast, Malpaís Las Víboras sites display a large scatter in NRM intensity values, while their susceptibility values are coherent (Figure 5b). Sites Mlv4 and Mlv5 have the highest and lowest NRM intensities, respectively. For all sites, susceptibility values are on geometrical average 10.9 × 10−3 (SI). In this context, it is worth noting that sites Mlv1, Mlv2, Mlv3, and Mlv6 were sampled roughly along the same level (middle part of the lava exposures), Mlv4 was sampled from the uppermost part, while Mlv5 samples were taken from the lowermost part. Therefore, their NRM intensities may have resulted from cooling-rate disparities. The cooling rate was lower in Mlv5 samples and lead to the growth of larger pseudo-single-domain (PSD) to multi-domain (MD) magnetic grains. However, the presence of small PSD to single-domain (SD) grain sizes is expected from Mlv4 samples because of their relatively faster cooling rate. However, near-surface Mlv4 samples might have also been affected by lightning strikes, which produce strong isothermal remanence magnetizations (IRMs), which may account for rather high Q values. In comparison with Malpaís Las Víboras, much smaller scatter in the NRM intensities is noticeable in the Capaxtiro sites, in particular considering that samples were obtained at different levels within the lava exposures (Figure 5c). Average values for the NRM intensity and susceptibility are 9.0 A m−1 and 10.2 × 10−3 (SI). Both properties together define Q values between 10 and 100 (Figure 5c), with the possibility of lightning strikes affecting sample Cpx7. Finally, convergent values were observed in Malpaís Prieto samples (Figure 5d) with average NRM intensity and susceptibility values of 11.2 A m−1 and 9.3 × 10−3 (SI), respectively. Most Malpaís Prieto specimens have Q values ranging between 10 and 100 (Figure 5d).

Variation of NRM intensity and magnetic susceptibility for the four studied lava flows. Dashed lines define constant values of Königsberger’s factor Q.
Thermomagnetic analyses indicate the presence of two types of behaviors with characteristic Curie temperatures (Tc) and distinctive degree of reversibility between the heating and cooling curves (Figure 6). Samples of type A were characterized by single, high Curie temperature (Tc ~ 530–560°C) and with a decrease <10% in magnetization upon cooling to room temperature (Figure 6a). This behavior is typical for Ti-poor titanomagnetite, which partially oxidized during heating. We note that type A was observed in seven out of eight cores. Thermomagnetic curve of type B (I core) points to the presence of both Ti-rich (Tc ~ 300–330°C) and Ti-poor (Tc ~ 500–540°C) titanomagnetite minerals (Figure 6b). During cooling, the low-Tc component was largely suppressed, probably by exsolution of this mineral, which produced low-Ti titanomagnetite, with a slight increase in the room temperature magnetization.

Two types of thermomagnetic curves: (a) type A and (b) type B. Black and gray curves indicate heating and cooling curves, respectively.
Paleomagnetic directions
In order to obtain reliable information on the flow mean directions, one specimen of every drill core (total of 141) was utilized. The direction of the ChRM was defined for each sample using PCA (Kirschvink, 1980; program PMGSC 4.2). Predominantly, the ChRM directions have been calculated by 7–10 vector end points and are characterized by maximum angular deviation (MAD) values between 3.2° (few specimens from El Infiernillo and Capaxtiro) and 0.5° and often <1°. Flow mean directions on core level and site level were calculated using Fisher statistics (Fisher, 1953, using PMag Tools Version 4.2) and are listed in Table 1.
For El Infiernillo, we point out that single components of magnetization were recorded in most of the investigated cores (Figure 7a, Inf2-5z). Samples displaying this behavior have relatively high median destructive fields (MDFs) of 40–55 mT. Few specimens of Inf3 showed minor secondary magnetization components that are probably of viscous origin and were removed by AF amplitudes of 5–10 mT (see Inf3-3b in Figure 7b). An overall flow mean direction was determined from all three sites of the El Infiernillo flow: Dec = 356.7°, Inc = 52.0°, and α95 = 3.2° (Table 1 and Figure 8a).

Orthogonal vector plots of AF demagnetized samples from the four studied lava flows. Labels along curves denote the maximum AF amplitude applied during the demagnetization steps. (a, b) El Infiernillo, (c, d, e) Malpaís Las Víboras, (f, g) Capaxtiro, and (h, i) Malpaís Prieto.

Characteristic remanent magnetization directions for the four lava flows. Inset in (b) shows the dispersed ChRM directions of Mlv4 (green) and Mlv6 (gray; for interpretation, see text). (a) El Infiernillo, (b) Malpaís Las Víboras, (c) Capaxtiro, and (d) Malpaís Prieto.
In regard to the Malpaís Las Víboras flow, no results were obtainable from site Mlv4 where we note that all seven specimens show large secondary components during the demagnetization. It is important to mention that some specimens of this site were also characterized by high Q values indicating lightning strike effects (see strong magnetic overprint in Figure 7c). Although individual demagnetization diagrams were stable and of high quality (small MAD values for the best fits), we could not obtain a site mean direction as the ChRM were strongly scattered (see inset in Figure 8b). At site Mlv6, all 11 specimens yielded dispersed ChRM directions (which probably was related to moved blocks) with a site mean direction of high α95 = 16.5° (Table 1 and inset in Figure 8b). Thus, it is clear that directions provided by both sites are not reliable and were not included in the calculation of the flow mean direction. Univectorial magnetization trends toward the origin were observed in most samples from Mlv2, Mlv3, and Mlv5 (see Mlv5-10b in Figure 7c). Their MDFs range from 45 to 60 mT. In most cases, samples from Mlv1 have a small viscous component that was easily removed by AF amplitudes of 3–7 mT (Mlv1-9c in Figure 7d). Varying MDF values displayed by samples from this site point toward the simultaneous presence of soft and hard magnetic minerals. These four sites have similar ChRM directions (Table 1) and define an overall flow mean direction with a low dispersion: Dec = 1.3°, Inc = 18.3°, and α95 = 2.3° (Figure 8b).
In the case of Capaxtiro, five sites (Cpx2, Cpx3, Cpx4, Cpx5, and Cpx6) were used to calculate the flow mean direction (Table 1). At these sites, a small percentage of samples show low-stability secondary components probably of recent origin that were easily removed by AF amplitudes of 5–15 mT (Cpx2-5a in Figure 7e). These samples have relatively high MDF of 50–70 mT. No site mean direction could be calculated for Cpx7 because it was probably affected by lightning strike and/or moved blocks (see Cpx7-4b in Figure 7g). Cpx1 gave a site mean direction with a relatively large uncertainty of α95 = 7.9° (Table 1), which is also significantly different from the site mean directions calculated from the aforementioned sites. After excluding Cpx1 (shaded row in Table 1, gray dots in Figure 8c), we were able to obtain a well-defined flow mean direction for Capaxtiro (Dec = 353.0°, Inc = 30.7°, and α95 = 3.5°).
Most sites of Malpaís Prieto are characterized by a single-component NRM (Figure 7g and h). Samples from these sites show similar NRM intensities, susceptibilities (Figure 4d), and MDFs of 30–45 mT during AF demagnetization. ChRM directions calculated for sites Mps1, Mps3, and Mps4 are very similar and define a flow mean direction with a lower dispersion (Dec = 356.8, Inc = 27, and α95 = 2.1°; see Table 1 and Figure 8d). However, site Mps2 displays a more easterly and dispersed mean direction (Dec = 38, Inc = 18.6, and α95 = 8.3°; gray dots in Figure 8d). Hence, we deduce that Mps2 samples were obtained from moved blocks and that, therefore, their ChRM directions are not reliable and were thus rejected. Because of the dense vegetation, it was difficult to determine during fieldwork whether this site contained displaced blocks. For Mps5, no site mean direction was obtained, as all ChRM directions were highly dispersed.
Paleointensities
IZZI-Thellier PI experiments were conducted on 40 samples obtained from the four lava flows. Pre-selection of specimens was done based on two criteria: (1) specimens from the same drill core must possess a single-component NRM during demagnetization, and (2) specimens must show a relatively high MDF of >20 mT. Laboratory fields of 50 or 60 µT were applied during heating and cooling of infield steps. According to Biggin (2010) and Paterson et al. (2015), the greater the angle between the laboratory field and the primary NRM, the greater the influence of MD grain size. Therefore, inside the oven, selected specimens for PI measurements were oriented so that their NRM directions would be parallel to the applied field with a precision better than 5°. Representative examples of successful Arai plots are shown in Figure 9a and b. PI results for each flow are listed in Table 2 together with different quality parameters. The Arai plots of accepted specimens yielded well-constrained straight lines with 6–14 data points (N), the fractions factor (f) ranges from 39% to 96%, and the quality parameter (q) varies between 6.50 and 76.9. In regard to the directional variations during the PI experiments, the ChRM produced MADs (MADanc; anc = anchored to the origin) between 4.1° (Inf1-8b) and 1.5° (Mlv2-4z; Table 2). This range is only slightly higher than the MAD range obtained during the AF demagnetization. In total, 20 specimens passed the Thellier tool A or B acceptance criteria resulting in an overall success rate of 50%. But when compared for the different flows, samples from Malpaís Las Víboras have the highest success rate (70%) while those from Capaxtiro have the lowest (40%) success rate. In most cases, failure of the test was due to alteration that resulted from the repeated heating steps as indicated by the pTRM cumulative check difference (dpal) criterion (Figure 9c). A few samples from Capaxtiro failed the test because the fraction of unblocked NRM was too small as indicated by the fraction (f) factor (Figure 9d). More importantly, our investigated samples did not show sagging in the Arai plots, which in turn mean that the MD effect was absent or has been suppressed, as suggested by Biggin (2010) and Paterson et al. (2015).

Examples of (a, b) successful and (c, d) unsuccessful paleointensity (Arai) plots for the four lava flows obtained by the IZZI version of the Thellier method. NRM and pTRM are normalized. NRM versus pTRM data are shown as circles, with the black best-fit line. pTRM checks are shown as triangles. The analyses were carried out using ThellierTool.
The Thellier-IZZI paleointensity results and associated statistics.
N: number of points included in the linear best fit; T: the temperature ranges used for the best fit; f: fraction of the natural remanent magnetization (NRM); g: the gab factor; q: quality factor; MADanc: anchored maximum angular deviation; α: angular difference between anchored and nonanchored best solution; δCK: relative check error; δpal: cumulative check difference; δTR: tail check; δt*: normalized tail of pTRM; PI: paleointensity; σ (µT): standard deviation.
Paleomagnetic dating
To date, a well-constrained regional paleosecular variation curve for central Mexico for the last few millennia, suitable for paleomagnetic dating, does not exist. For this reason, we used the global model SHA.DIF.14k (Pavón-Carrasco et al., 2014), which was developed by using archeomagnetic and lava flow data obtained from sites around the globe, including Central America, Mexico, and the United States, covering the last 14 millennia. Paleomagnetic dating was carried out with the help of the MATLAB tool archaeo_dating (Pavón-Carrasco et al., 2011). In a recent paleomagnetic dating study (Böhnel et al., 2016), this model was successfully applied on two lava flows from Ceboruco volcano in western Mexico. This motivated us to use this dating technique on the four lava flows under consideration in this study. Knowledge of the age of El Infiernillo (1525–1420 cal. BC, see above) together with the paleo-direction and intensity data listed in Tables 1 and 2 allows constraining the time interval of archeomagnetic dating to the period 2000 BC to AD 1900. For the paleomagnetic dating, flow mean directions on the core level were used, as individual sites were sampled with variable numbers of cores and also showed large variations of within site dispersion: in the case of the El Infiernillo flow the number of ChRM directions varied between 5 and 10, and precision parameters k between 64.25 and 438.94. For the other flows, similar variations were observed. The use of unweighted site mean directions would overestimate those sites with small number of ChRM directions and/or large dispersion. Finally, the number of sites per flow is small for all flows (three or four sites), and only for Capaxtiro, five site means are available. As all sites could be assigned beyond any doubt to the corresponding lava flow, the dispersion of the individual ChRM directions should be caused by random processes (e.g. Böhnel and Schnepp, 1999), and the overall mean then may also average out minor block movements that may have occurred along the lava flow. Thus, we argue that the use of individual ChRM directions for calculating a flow mean direction is justified in this specific situation. Nevertheless, in Table 1, we list both flow mean directions, and it is evident that in most cases, the uncertainty α95 based on site mean directions is about twice compared with that based on individual ChRM directions. The exception is Malpaís Prieto, where both overall mean directions have a very similar α95. Moreover, in each flow, the site mean directions were statistically evaluated using the F-distribution test (McFadden and Lowes, 1981) in order to analyze their directional independence at the 95% confidence level. Applying the F-test on El Infiernillo sites proved that the site mean directions are indistinguishable at the 95% confidence, which is also valid for the Malpaís Prieto sites. In contrast, the F-test applied to Malpaís Las Víboras mean directions shows that they are different at the 95% confidence. It is likely that the difference in declination values of Mlv5 is the main cause of this negative outcome, but we also note that some specimens of Mlv1 have similar declinations as Mlv5. Applying the F-test on the Capaxtiro also reveals that the site mean directions are different at the 95% probability level. Previous studies (e.g. Hagstrum and Champion, 1994) have shown that such differences between site mean directions of the same lava flow may occur, due to a number of possibly contributing factors. In case of this study, these are most likely due to undetected relative movements between different sites. Therefore, we attribute these within-flow paleomagnetic variations to the natural dispersion (see Böhnel and Schnepp, 1999). As no information is available regarding the most reliable site mean direction, we assume that the overall mean directions calculated from all individual directions may best correspond to the accurate paleomagnetic direction of that particular flow. The flow mean direction calculated this way has a larger uncertainty, which thus will only increase the range of the resulting paleomagnetic dating age. In conclusion, we consider the flow mean paleomagnetic directions calculated for Malpaís Las Víboras and Capaxtiro to be reliable within the obtained confidence limits and to be acceptable input data for the paleomagnetic dating procedure which will provide an age range on the safe side.
Paleomagnetic dating provided for El Infiernillo a well-constrained age range of 1500–1370 BC (Figure 10a; to take into account the limitations of this dating method, here we will round up the ages to the nearest 10 years). This age range matches well the calibrated radiocarbon age range and lends further credibility to the paleomagnetic dating method. Using the mean direction based on site means (Table 1) would increase the age range to 1590–1360 BC, and additional possible age ranges of 1800–1700 BC, AD 10–60, AD 300–650, and AD 1890–1900 would result. Even then, there would be an agreement between the C-14 age and one of the possible paleomagnetic ages. Based on the stratigraphy, the three other lava flows are younger than El Infiernillo, and for their paleomagnetic dating, we constrained the time interval of their emplacement to the period 1500 BC to AD 1900. Two narrow age ranges of 1340–1230 BC and 1030–940 BC (Figure 10b) were obtained for the Malpaís Las Víboras. This means that it erupted shortly after or up to 600 years after El Infiernillo. For Capaxtiro and Malpaís Prieto, two ages were obtained in both cases. Based on the stratigraphic relationships between these flows, it can be proposed that Capaxtiro erupted within the period 200–80 BC (Figure 11a), followed by Malpaís Prieto during the interval AD 830–960 (Figure 11b). The age range of this youngest flow would remain almost unchanged when using the mean direction based on site means for the dating.

Paleomagnetic dating of (a) El Infiernillo and (b) Malpaís Las Víboras. The combined probability density derived from the declination, inclination, and PI data is shown as shaded peaks and the minimum (95%) confidence level by horizontal green lines.

Paleomagnetic dating of (a) Capaxtiro and (b) Malpaís Prieto. For further details, see caption of Figure 9.
Discussion
Rock-magnetic studies were carried out on the four flows under study in order to demonstrate the NRM stability and also to define the magnetic carriers and thermal stability, which is crucial for the PI experiments. These studies revealed significant differences between the studied lavas in terms of their bulk magnetic properties (NRM intensities and susceptibilities). It was also observed that each studied flow showed variations of these properties to different degrees. Malpaís Las Víboras shows the largest disparities in its values, whereas Capaxtiro and Malpaís Prieto display a greater coherency. This is probably due to differences in composition and crystal size and to processes such as lightning strikes affecting the TRM. As a whole, their Königsberger’s factors were located within the 100 > Q > 10 domain, which indicates the high NRM stability over geological time periods and its suitability for paleomagnetic dating. Thermomagnetic analyses showed that stable Ti-poor titanomagnetite is the predominant magnetic mineral in these lavas. This conclusion is admittedly based on a limited number of samples, but a more detailed rock-magnetic study was beyond the main goal of this study.
As mentioned above, for most sites, ChRM directions could be determined. Nevertheless, some sites had to be eliminated because either these directions were dispersed, or they differed strongly from other sites of the same flow. Mean directions of accepted sites correlate well within their confidence limits. The mean directions for the four lava flows used for paleomagnetic dating are always based on the ChRM directions of all accepted samples from that flow, but in Table 1, we also indicate for comparison mean directions based on the site mean directions, which are very similar but have a larger α95. We always used individual directions, as there was no doubt that the corresponding sites belong to one single lava flow.
In the case of El Infiernillo lava flow, site Inf3 has a more westerly mean direction than Inf1 and Inf3, but this site also has a larger dispersion with a α95 of 9.6°, which makes this difference insignificant. Site Mlv4 from Malpaís Las Víboras seems to have been affected by lightning strike, as indicated by high Q values, and ChRM directions were dispersed, and no site mean was determined. Partly, this dispersion may be also the effect of moved blocks, which is also suspected for site Mlv6. Both sites were rejected for further analysis. Hence, the overall mean direction for Malpaís Las Víboras is based on four sites, which are distributed over a distance of more than 2 km. For the Capaxtiro lava flow, six out of seven sites yielded well-defined ChRM directions (Table 1). Five sites (Cpx2 to Cpx6) distributed over a distance of ~120 m gave consistent site mean directions, whereas nearby site Cpx1 (less than 10 m away) has a significantly different direction, most probably representing a moved block. Cpx7 has large Q values and seems to be affected by lightning strike, and hence, its site mean direction was not determined. The overall flow mean direction is, therefore, based on five sites distributed over a distance of about 120 m. From the above, it seems that Cpx1 could represent a moved block and, therefore, was excluded from the overall flow mean calculations. Similarly, for Malpaís Prieto, two of five sites had to be rejected: Mps5 had dispersed ChRM directions and was sampled from several apparently moved blocks. Mps2 directions are less dispersed but define a significantly divergent mean direction compared to the rest of the flows. ChRM directions calculated from Mps1, Mps3, and Mps4 are concordant and best represent the Malpaís Prieto lava flow mean direction. These experiences lead us again to emphasize the utmost importance of sampling lava flow units at multiple sites and over the largest possible distance, as already suggested by Böhnel et al. (2016).
El Infiernillo was studied in this work in order to validate the archeomagnetic dating method by using the global model SHA.DIF.14k. Among the four studied lava flows, this flow is stratigraphically the oldest with a 14C age of cal. 1525–1420 BC. The concordant paleomagnetic age of 1500–1370 BC thus suggests that this volcano erupted during the early Pre-Classic, and its birth was certainly observed by early settlers in this region. Malpaís Las Víboras erupted shortly or up to 560 years after El Infiernillo, incrementing the volcanic impact on the Middle Pre-Classic local societies. The voluminous Capaxtiro lava flows were emplaced about 740–1260 years after the Malpaís Las Víboras, during the late Pre-Classic. This would imply a large disruption for the more evolved human communities inhabiting these lands. Approximately 910–1160 years later, between AD 830 and 960 (during the late Classic period), the Malpaís Prieto flow erupted. Then, extended settlements were distributed all over central Mexico, and while this lava flow covered a much smaller area than the previous volcanoes, it must have had an important impact on these more developed and extended societies. Regarding recurrence rates of volcanic eruptions, these varied around 1000 years, with the exception of Malpaís Las Víboras, which erupted much sooner after El Infiernillo. If these recurrence rates are valid, another eruption could happen soon.
The different eruptions that formed the Malpaís de Zacapu probably had an impact on the settlement process, provoking displacements within the Zacapu basin area. In regard to the Infiernillo and Las Víboras eruptions, it is difficult to assess their repercussions on early/middle Pre-Classic populations because archeological evidence for that time is still lacking. But with respect to the Capaxtiro eruption (dated here at 200–80 BC), it is interesting to note that it coincides with the earliest archeological sequence in the Zacapu basin, the Loma Alta phase, dated at 100 BC to AD 400. During this phase, a series of islets known as Las Lomas just east of the Capaxtiro lavas were densely inhabited (Arnauld et al., 1993). So far, we can only guess that a direct relation may exist between the Capaxtiro eruption and the rise of the important site of Loma Alta.
Finally, in the case of the Malpaís Prieto eruption, more consistent data are available in regard to its possible impact on local populations. Recent excavations and surveys around this lava flow show that the area was densely populated between AD 500–800/900 (six important sites have been reported so far). All of these sites were abandoned abruptly at the end of this period, which was then followed by a hiatus that lasted from AD 900–1250. Hence, a close relationship between this eruption and the abandonment of the area is quite likely. Around AD 1250, the Malpaís de Zacapu becomes again an important population hub. We note here that many population centers in Mesoamerica collapsed around AD 900–1000, a phenomenon that has been explained by famine due to the effects of long-term regional droughts (e.g. Bhattacharya et al., 2015; Hodell et al., 2005). In the case of Malpaís de Zacapu such a climate-related impact on ancient societies remains to be proven, while the possible volcanic impact is hard to deny in light of the above-mentioned evidence.
Conclusion
Lava flows from four different Holocene monogenetic eruptions collectively known as Malpaís de Zacapu (El Infiernillo, Malpaís Las Víboras, Capaxtiro, and Malpaís Prieto, in stratigraphic succession) located at the western margin of the Zacapu lacustrine basin in Michoacán were subjected to paleomagnetic dating. Because this region (including the surface of the lava flows) was inhabited since at least the Pre-Classic and eventually became part of the heartland of the Tarascan Empire during the Post-Classic, it was desirable to determine the timing of these eruptions and whether they could have impacted human development, that is, by triggering population migrations. Only the oldest eruption, El Infiernillo, could be dated by the radiocarbon method, which yielded an age of cal. 1525–1420 BC. A total of 21 sites distributed as far as possible were sampled from four lava flows in order to provide reliable paleomagnetic site mean directions which are not affected by blocks that moved after they acquired their remanence or by lightning strike–induced remagnetizations. Flow mean directions are of small dispersion with 2.1 ⩽ α95 ⩽ 3.5. Robust estimates of flow mean PIs were obtained by using the IZZI-Thellier method. Full-vector paleomagnetic results were used for paleomagnetic dating applying the MATLAB tool archaeo_dating and the global paleosecular variation model SHA.DIF.14k. For El Infiernillo, the dating resulted in an age ranging 1500–1370 BC (95% probability level) which coincides well with the radiocarbon age data, and both correspond to the early Pre-Classic period. For the next younger lava flow of Malpaís Las Víboras, a possible age between 1340 and 940 BC was obtained, corresponding to the Middle Pre-Classic, and reflecting the sustained volcanic activity in this region, which may have affected the ancient population. After 740–1260 years of quiescence, the voluminous lava flows of Capaxtiro were emplaced between 200 and 80 BC during the late Pre-Classic, implying a disaster for communities living nearby. This period coincides with the rise of the important site of Loma Alta in the lowlands east of Capaxtiro. Finally, between AD 830 and 960 and during the late Classic period, Malpaís Prieto erupted, which apparently led to the abandonment of the heavily populated surrounding lava flows as observed in archeological profiles which indicate that this happened around AD 800–900. Long drought periods as observed in eastern Mexico and the Yucatán peninsula might have amplified the volcanic impact. Finally, recurrence periods center around 1000 years in this area and thus indicate that another monogenetic eruption might occur in the near future. This work shows that future volcanic hazard mitigation efforts could benefit from results obtained from interdisciplinary studies that include archeological research.
Footnotes
Funding
Ing. J. Escalante supported studies with the Curie balance. Field and laboratory costs of ANM and HB were covered by Consejo Nacional de Ciencia y Tecnología (CONACyT-180032) and the Dirección General de Asuntos del Personal Académico (UNAM-DGAPA IN-111915) granted to HB and for NR-G and CS were defrayed from projects funded by the Consejo Nacional de Ciencia y Tecnología (CONACyT-167231) and the Dirección General de Asuntos del Personal Académico (UNAM-DGAPA IN-101915) granted to CS. We thank Américo González-Esparza, Adriana Briseño, and Ernesto Andrade at UNAM Campus Morelia for providing lodging facilities at the Mexican Array Radio Telescope (MEXART) near Coeneo during field campaigns. The archaeological investigations carried out in the Zacapu region were funded by the French Ministère des Affaires Etrangères et du Développement International (Mission Archéologique Uacusecha) as well as by the Agence Nationale de la Recherche (Programme Mesomobile).
