Petrogenesis of voluminous silicic magmas in the Sierra Madre Occidental large igneous province, Mexican Cordillera: Insights from zircon and Hf-O isotopes

This is a study of the petrogenesis of voluminous rhyolites in the Sierra Madre Occidental (SMO) silicic large igneous province (SLIP) of northwestern Mexico (Fig. 1) through the lens of Lu-Hf and O isotopes in zircons dated by U-Pb geochronology. The SMO is the largest Cenozoic SLIP (Bryan and Ferrari, 2013) and is therefore an important natural laboratory to examine the petrogenesis of voluminous silicic magmas. Moreover, the SMO is an atypical SLIP not associated with continental breakup (e.g., Chon Aike) nor with a mantle hotspot (e.g., Yellowstone). Specifically, this is a regional-scale study of the northern and central SMO that tests competing hypotheses for the formation of the eponymous Eocene to Miocene ignimbrite flare-up events and offers insights into the formation of voluminous silicic magmas in an extending continental arc setting (e.g., Taupo, etc.). Through analysis of 49 rhyolitic to dacitic ignimbrites and lavas, we reveal ~25 m.y. of monotonous silicic magma production associated with fractionation of mantle-derived juvenile melts mixing with lower-crustal anatectic melts throughout the SMO. We demonstrate that Hf and O isotopic values are compatible with a largely primitive source for melts and heat in the mantle and only ~20% recycling of paragneissic continental crust. This analysis supports one set of models and strongly refutes others that argue for a dominantly crustal recycling petrogenesis and, therefore, emphasizes the importance of material and heat transfer from the mantle to the continental crust during arc rifting.

The middle Eocene to early Miocene SMO (Fig. 1) is the most voluminous (~3.9 × 10⁵ km³; Bryan et al., 2008) Cenozoic SLIP on Earth (Bryan and Ferrari, 2013) and is host to the world's largest epithermal Ag province (Camprubí et al., 2003). The ash dispersed during the SMO ignimbrite eruptions has been suggested as a driver of global cooling at the Eocene–Oligocene transition and the onset of Antarctic glaciation (Cather et al., 2009; Jicha et al., 2009). The SMO buries ~400,000 km² beneath a 1–1.5 km thickness of silicic ignimbrites and lavas (Fig. 1); McDowell, 2007; Bryan et al., 2008), forming a central plateau with mean elevation 2000 m above sea level (Fig. 2). The margins of the SMO are extended, tilted, and truncated by normal fault-bounded graben and half-graben of the southern Basin and Range province in the east (Henry and Aranda-Gomez, 1992, 2000) and the Gulf of California extensional province in the west (e.g., Ferrari et al., 2013). The elevated core is less extended and is characterized by subhorizontal ignimbrite sheets (>10 km³ each) exposed for tens of kilometers in canyon walls and on peaks and ridges, with isolated clusters of contemporaneous lava domes and several large calderas (e.g., Swanson et al., 2006).

The SMO forms the southern and largest part of an ignimbrite flare-up that extended from British Columbia to Mexico from the Eocene to the earliest Miocene, widely interpreted to result from rollback of the Farallon slab that had been subducting under North America since the Jurassic (e.g., Sedlock et al., 1993). The silicic SMO arc developed on the remnants of the Paleocene–Eocene continental arc that had formed during slab flattening, which in turn, was constructed on the Mesozoic continental margin (Fig. 1). However, the SMO is understudied compared to other parts of the Cordilleran arc, and much of it remains unexplored: even basic knowledge of its internal architecture, stratigraphy, and geochronology is lacking. What understanding there is comes from studies along the few roads that exist, and a few regional-scale geophysical studies (e.g., Bonner and Herrin, 1999; Purucker et al., 2007). Many published studies of the northern and central SMO are constrained by either a small number of samples (e.g., Bryan et al., 2008), a small study area (e.g., Murray et al., 2013), or both (e.g., Mahar et al., 2019), limiting understanding of regional variations and patterns. The available age data are relatively sparse compared to the total area and the total age range and are probably biased by accidental over-sampling of some widespread units (e.g., Vista Tuff; McDowell and McIntosh, 2012) and under-sampling in less accessible areas and in poorly exposed, non-welded units. The few regional-scale compilations available are focused on K/Ar and Ar/Ar geochronology and the first-order geological history of the SMO (Ferrari et al., 2007; McDowell and McIntosh, 2012; Ferrari et al., 2018).

Silicic volcanism generally migrated southwestward after slab rollback began at ca. 40 Ma (Damon et al., 1981). The northern and central SMO experienced several localized, ~2-m.y.-long, ignimbrite flare-ups throughout the period with ca. 46–42 Ma and 38 Ma (mainly north of Chihuahua City; Fig.2A), 36–28 Ma (widespread), and ca. 24–18 Ma (southwestern margins) being significant timespans although not necessarily discrete events (McDowell, 2007; McDowell and McIntosh, 2012). The largest flare-up pulse was at 36–28 Ma, with a total volume flux of ~3.0 × 10⁵ km³ (Bryan and Ferrari, 2013). The late Oligocene pulse (25–23 Ma) encroached on the western margins of the northern and central SMO where it is unconformable upon and recycles zircon from the much more voluminous middle Eocene–middle Oligocene (36–28 Ma) volcanic succession (e.g., Murray et al., 2013; Mahar et al., 2019). This pulse is contemporaneous with the oldest parts of the predominantly early Miocene (24.5–18 Ma) ignimbrite flare-up in the southernmost SMO (Ferrari et al., 2007; Aguirre-Díaz et al., 2008; Bryan and Ferrari, 2013; Ferrari et al., 2013), which we will not discuss further.

Volcanism appears to have developed as ~100–150-km-wide clusters of similar-aged calderas, lava dome complexes, and fissure vents at different times along each transect (McDowell and Keizer, 1977; Swanson and McDowell, 1984, 1985; Swanson et al., 2006; McDowell and McIntosh, 2012; Ferrari et al., 2013; Murray et al., 2013). The relative paucity of calderas identified in the SMO compared to the number of thick and extensive ignimbrites suggests that either many calderas remain to be discovered (e.g., Swanson and McDowell, 1984; Henry et al., 2012; McDowell and McIntosh, 2012) or that many eruptions were fed by fissure vents bounding or within graben (Aguirre-Díaz and Labarthe-Hernández, 2003; Torres-Hernández et al., 2006; Aguirre-Díaz et al., 2008; Tristán-González et al., 2008; Murray et al., 2013).

The majority of SLIPs are associated with either intracontinental hotspot volcanism (e.g., Yellowstone–Snake River Plain), continental rifting (e.g., Whitsunday), or both (e.g., Etendeka-Paraná), and may or may not be associated with voluminous mafic volcanism (Bryan and Ferrari, 2013). The SMO is unusual among SLIPs for being an arc; although large calderas and large-volume rhyolitic eruptions occur in arc settings (e.g., Altiplano-Puna, Bolivia; Taupo volcanic zone, New Zealand; Kyushu arc, Japan), none match the scale and duration of the SMO. Attempts to explain the tectonic influences on the petrogenesis of SMO magmas fall into two categories: (1) waning of the Cordilleran arc accompanied by rollback of the Farallon slab (e.g., Wark, 1991; McDowell and Mauger, 1994; McDowell and McIntosh, 2012); and (2) breakup of the Farallon slab and opening of a slab window beneath western Mexico (Mahar et al., 2019), accompanied by waxing lithospheric extension leading to seafloor spreading in the Gulf of California (Ferrari et al., 2018, and references therein).

Geochemical studies of the SMO rhyolites, coeval mafic lavas, and Cenozoic xenolith localities in the 1980s and 1990s were dominated by whole-rock elemental and isotopic (Sr and Nd) analyses and modeling to test between petrogenesis through assimilation of continental crust, fractionation of mantle-derived mafic magmas, or a combination of the two (assimilation-fractional crystallization [AFC]). Such studies produced contradictory interpretations for rhyolite petrogenesis. For example, Ruiz et al. (1988b) inferred up to 100% melting and recycling of Paleozoic and older lower continental crust, and Bryan et al. (2008) interpreted zircon solubility and age data as evidence in favor of crustal melting. Alternatively, relatively closed-system fractionation of arc-derived mafic and intermediate magmas with ≤20% lower-crustal contamination satisfies AFC mixing models in several studies (e.g., Lanphere et al., 1980; Cameron and Cameron, 1985; Cameron et al., 1986; Wark, 1991; Smith et al., 1996). Analysis by Mahar et al. (2019) of Hf isotopes in zircons from southwestern Chihuahua suggests that 90–50 Ma plutons and early Miocene SMO rhyolites share a common source in the subcontinental lithospheric mantle and were produced by AFC processes. Analysis by Ferrari et al. (2018) of Oligocene to Pliocene mafic and silicic rocks throughout northwestern Mexico reveals a transition from calc-alkaline magmas that incorporated continental crust (low-negative εNd) before 20 Ma and increasingly tholeiitic, uncontaminated, mantle-derived magma after that. They interpret this to reflect the cessation of arc magmatism and the increasing dominance of crustal extension and eventual tapping of the asthenospheric mantle.

This paper focuses on the isotopic geochemistry of zircon to elucidate the petrogenesis of SMO magmas and to assess competing tectonomagmatic models. Combined zircon U-Pb, Hf, and O studies are excellent tools for studying the origin of igneous and metaigneous rocks and have been applied to high-grade metamorphic assemblages (e.g., Zheng et al., 2006; Hiess et al., 2009, 2011) and batholith complexes (e.g., Kemp et al., 2007; Appleby et al., 2010; Bolhar et al., 2012; Zhao et al., 2013) to evaluate crustal growth mechanisms and rates (e.g., Hawkesworth and Kemp, 2006; Kemp et al., 2006). This O-Hf-U/Pb in zircon method has been applied to only a few silicic volcanic systems including Yellowstone–Snake River Plain SLIP (e.g., Stelten et al., 2013; Colón et al., 2018) and Iceland (e.g., Banik et al., 2018; Carley et al., 2020); this is the first application to the SMO.

Our samples come from the northern and central SMO collected from three transects: northern, central, and southern (Figs. 1 and 2), all in the ancestral homelands of the displaced Rarámuri people. From north to south: (A) the ~250-km-long northern transect, which extends from Chihuahua City to the Hermosillo area (Cochemé and Demant, 1991; Swanson et al., 2006; McDowell, 2007; Murray et al., 2013, 2015; Murray and Busby, 2015); (B) the ~150-km-long central transect in the Chihuahua-Durango border region west of Parral; and (C) the ~100-km-long southern transect from Durango to just east of Mazatlán (Aguirre-Díaz and McDowell, 1991, 1993; Ferrari et al., 2007; Aranda-Gómez et al., 2015). This is the first study of the central transect. All three transects include samples from the high-elevation SMO plateau and the adjacent rifted eastern margins of the southern Basin and Range (Fig. 2), thus allowing us to look for variations across the complete extent and timespan of the SMO and across the structural boundaries of the non-extended core.

Zircon (ZrSiO₄) is uniquely able to record isotopic changes during rhyolite petrogenesis because it is abundant, refractory, and preserves useful textural information. It is often one of the only phases present in magmas erupted close to their liquidus. Zircon readily incorporates Hf over Lu (Lu/Hf ~0.002) in substitution for Zr (0.5–2.0 wt%; Kinny and Maas, 2003) and has a strong affinity for U but not for Pb, making it ideal for both U-Pb geochronology and Hf isotope analyses. Importantly, zircon allows for only very slow diffusion of O and cations; therefore, O and Lu-Hf isotope ratios in zircon are typically not reset by metamorphism or metasomatism (Valley et al., 1994; Kinny and Maas, 2003; Page et al., 2007; Bowman et al., 2011).

Oxygen and Hf isotopes inform on the sources of magmas (Hawkesworth and Kemp, 2006). Hafnium isotopes are radiogenic and are sensitive to the ages of rocks being melted (Kinny and Maas, 2003); conversely, O isotopes are stable and reflect the specific oxygen isotope compositions of the melts from which they crystalized (Valley, 2003). Mantle-derived magmas are expected to have strongly positive εHf values and δ¹⁸O_zircon values ~5.3‰ ± 0.3 (1σ; Valley, 2003). Anatectic magmas are expected to have widely variable δ¹⁸O_zircon values and εHf values correlated to the nature of the protolith: low-positive to low-negative for melting of metaigneous rocks and strongly negative for metasedimentary rocks.

Zircon separates were prepared by standard methods at the University of California Santa Barbara (UCSB) and California State University Bakersfield (CSUB). Samples were crushed by hand with a mortar and pestle and then by jaw crusher and BICO disk mill. Sand-sized crush was sieved to remove grains greater than 200 μm and then was density separated on a water table. The densest fractions underwent magnetic separation, followed by separation of apatite from zircon by flotation in methylene iodide (MEI). Zircon concentrates were mounted in Epofix epoxy and imaged by scanning electron microscopy and cathodoluminescence (CL) at CSUB, the University of Wisconsin Madison, and the University of Idaho.

Zircon recovery varied from minimal (≤5 grains per sample) to abundant (1000 s), with no systematic correlations between abundance and age or location. Zircon are uniformly prismatic and oscillatory zoned under CL (Fig. 3). Grains typically range from ~40–150 mm in length and are usually euhedral or subhedral. Many grains contain melt inclusions. High CL contrast cores are common and subhedral, but only very rarely display internal zoning discordant to the outer rims. Zircon grains were analyzed for U-Pb, Hf, and O isotopes in the same core and rim locations determined by CL imaging.

Oxygen isotopes in zircons were measured at the University of Wisconsin-Madison (Kita et al., 2009; Valley and Kita, 2009) and the University of California Los Angeles secondary ion mass spectrometry (SIMS) laboratories.

Samples from the northern and southern transects were analyzed on the CAMECA IMS-1280 ion microprobe at the WiscSIMS Laboratory, University of Wisconsin-Madison against the KIM-5 zircon standard (δ¹⁸O = 5.09 ± 0.06‰ (1σ) Vienna standard mean ocean water [V-SMOW]; Valley, 2003). Samples were cast within 7 mm of the center of 25.4-mm-diameter epoxy mounts, ground to approximate midsection, and polished. Sample relief was evaluated by profilometer. Data were acquired in four 12–15 h sessions following the procedures described previously (Kita et al., 2009; Valley and Kita, 2009). Three Faraday cup detectors were used to simultaneously analyze ¹⁶O⁻, ¹⁶O¹H⁻, and ¹⁸O⁻. Ratios of OH/O were background corrected against values on bracketing anhydrous KIM-5 and provide a monitor of “water,” which might result from inclusions or metamict zones within a zircon (Wang et al., 2014). Analyses were conducted using a 1.9–2.2 nA primary Cs⁺ beam with a spot diameter of 10 μm. The KIM-5 standard was analyzed four times between every batch of 10–15 sample spots to calibrate analyses and evaluate possible drift. The mean precision (1σ) of analyses is 0.15‰. Usually eight to ten analyses were made on up to 12 different grains.

Oxygen isotopes from samples from the central transect were analyzed on the CAMECA IMS-1290 ion microprobe at UCLA against the R33 (δ¹⁸O = 5.55 ± 0.04‰ (1σ) V-SMOW; Valley, 2003) and Mud Tank (δ¹⁸O = 5.03 ± 0.10‰ (1σ) V-SMOW; Valley, 2003) standards. Data were acquired in two 8–12 h sessions following the procedures described by Trail et al. (2007). Analyses were conducted using a focused 2–3 nA primary Cs⁺ beam with a spot diameter of 10–15 μm. Usually 10–12 analyses were made on up to 12 different grains.

Oxygen isotopes were measured in quartz phenocrysts from eight samples (five from the northern transect and three from the southern transect) by laser fluorination at the University of Wisconsin-Madison (Spicuzza et al., 1998; Valley, 2003). Anhedral, 100–200-μm-diameter quartz crystals were concentrated by HF digestion of a coarse, sand-sized aggregate, and then were handpicked. Samples were run against the UWG-2 standard (≈5.8 ± 0.15‰ (1σ) V-SMOW; Valley et al., 1995). Six quartz samples (four from the northern transect and two from the southern transect) correspond with analyses of coexisting zircon (samples CH03, CH97-7, N-99, SJ-54, K-LP-T, and SL-41).

Zircons were subsequently analyzed by laser ablation–inductively coupled mass spectrometry (LA-ICPMS) and SIMS to measure U-Pb and Lu-Hf analyzes.

Archived samples from the northern and southern transects were analyzed for U-Pb and Hf simultaneously in ~8–11 grains by laser ablation split-stream (LASS) ICPMS at Washington State University (40 μm); ~10–12 additional U-Pb ages were determined by LA-ICPMS (30 μm) for each sample. The LASS setup at Washington State University measures U/Pb and Lu-Hf isotopes simultaneously from a single ablation; the method is described in detail in Fisher et al. (2014). Samples are ablated using a NewWave 213 nm Nd:YAG laser using laser spot diameters of either 30 or 40 μm. The laser was operated at 10 Hz with a fluence of ~9–10 J/cm³. The ablated material is transported out of the laser cell using helium carrier gas and is mixed with a small amount (5 mL/min) of N₂ gas, which is added using a plastic “Y” split ~30 cm after exiting the laser ablation system. The ablated sample+He+N₂ mixture is then split into two separate paths, again using a “Y” split, an additional 30 cm downstream from the point of N₂ gas addition. One path carries the ablated aerosol to a ThermoScientific Element2 high-resolution ICPMS for U-Pb isotopic measurement, and a second path carries the ablated aerosol to a ThermoScientific Neptune MC-ICPMS for Lu-Hf isotopic measurement. Ar sample gas is added to each path ~20 cm from the torch.

Samples from the central transect were analyzed at UCLA and at the University of Arizona. U-Pb zircon geochronology was carried out in several analytical sessions on the CAMECA ims 1270 ion microprobe at UCLA, using established protocols for the analysis of Cenozoic zircon (Schmitt et al., 2003). Primary beam intensities for mass-filtered ¹⁶O⁻ were ~20 and ~40 nA with a spot size of ~20–30 μm with depths after analysis being approximately one-tenth of the crater diameter. Oxygen flooding was used at the saturation pressure to enhance sensitivity for Pb. A ²⁰⁴Pb-based common Pb correction was used for reference zircon, whereas ²⁰⁷Pb was used as a proxy for common Pb for unknowns. Reproducibility of ²⁰⁶Pb*/²³⁸U ages for AS3 (1099 Ma; Paces and Miller, 1993) used as a primary reference for the UO⁺/U⁺ versus U-Pb relative sensitivity calibration was ~2%–3% (relative), and accuracy of the weighted-average ²⁰⁶Pb*/²³⁸U ages based on comparison with replicate analyses (n = 10 each) on 91500 reference zircon was −1.5% (relative deviation from reported age of 1065 Ma; Wiedenbeck et al., 1995). Th/U in unknown zircon was determined from measured ²³²Th⁺/²³⁸U⁺ using a relative sensitivity factor calibrated on measured ²⁰⁸Pb*/²⁰⁶Pb* of AS3 reference zircon. Uranium abundances were estimated from U⁺/⁹⁴Zr₂O⁺ in relation to zircon 91500 with 81 ppm U (Wiedenbeck et al., 2004).

U-Pb analyses were corrected in IsoplotR (Vermeesch, 2018) for common Pb based on measured ²⁰⁷Pb, assuming ²⁰⁷Pb/²⁰⁶Pb = 0.8283 (Sañudo-Wilhelmy and Flegal, 1994) and disequilibrium based on ΔTh/U = 0.17 (Blundy and Wood, 2003). IsoplotR was used to calculate mean ages and to plot Tera-Wasserberg concordia.

Hf isotopic analyses were conducted on grains already dated at UCLA using a Nu Instruments HR multi-collector ICPMS connected to a Photon Machines Analyte G2 excimer laser at the Arizona LaserChron Center at the University of Arizona, following established analytical protocols described in Gehrels and Pecha (2014). Instrument settings were optimized for laser ablation analysis, and in-run analyses of seven different natural zircon standards (Mud Tank, Temora, FC52, R33, Plesovice, and Sri Lanka) were analyzed and monitored every 15–20 unknowns.

Ablation with a laser beam diameter of 40 μm targeted zircon domains directly over the earlier U-Pb analysis craters. All measurements were made in static mode using Faraday collectors equipped with 3 × 10¹¹ Ω resistors. Each acquisition consisted of one 40-second integration for backgrounds (on peak with laser idle) followed by 60 1-second integrations with the laser firing. A 30-second delay between laser analyses provided adequate time to ensure the previous sample was completely purged from the collector block. The laser was run in constant energy mode with an output of 7.5 mJ (fluence ~8 J/cm²) and a pulse rate of 7 Hz. Isotope fractionation (β) was accounted for by incorporating the method of Woodhead et al. (2004) with βHf determined from the measured ¹⁷⁹Hf/¹⁷⁷Hf and βYb from the measured ¹⁷³Yb/¹⁷¹Yb (except for very low Yb signals); βLu is assumed to be the same as βYb; and an exponential formula is used for fractionation correction. Yb and Lu interferences are corrected by measurement of ¹⁷⁶Yb/¹⁷¹Yb and ¹⁷⁶Lu/¹⁷⁵Lu (respectively), as advocated by Woodhead et al. (2004). Isotope ratios of ¹⁷⁹Hf/¹⁷⁷Hf = 0.73250 (Patchett and Tatsumoto, 1981); ¹⁷³Yb/¹⁷¹Yb = 1.132338 (Vervoort et al., 2004); ¹⁷⁶Yb/¹⁷¹Yb = 0.901691 (Vervoort et al., 2004; Amelin and Davis, 2005); ¹⁷⁶Lu/¹⁷⁵Lu = 0.02653 (Patchett, 1983) were used and all corrections are done line-by-line. Data reduction protocols account for all standards and unknowns analyzed during an entire session, and the βHf and βYb cut-offs were determined by monitoring the average offset of the standards from their known values resulting in a final standard offset of 0.0000001. ¹⁷⁶Hf/¹⁷⁷Hf uncertainty on the standard analyses for the entire session was determined to be 0.000037 (2σ).

A total of 49 samples (Table 1) were described (Supplemental Material, Section A¹) and analyzed including 21 from the northern transect, 22 from the central transect, and six from the southern transect (Fig. 2). The samples include silicic lavas and ignimbrites and encompass the full spatial and temporal range of volcanism in the northern and central SMO. The only unit knowingly duplicated was the Quemada Tuff (N-12 and N-88; Swanson et al., 2006).

Archived samples and mineral separates were analyzed from collections at the University of Texas Austin and the University of Texas San Antonio. Twenty-nine samples have a complete suite of δ¹⁸O, εHf_T, and U/Pb age measurements on individual zircon grains.

New samples along the central transect were collected during field work along Federal Highways 23 and 24 in 2013 and 2015. Figure 4 is a graphic log of a previously undocumented ~1.2-km-tall SMO section in a side canyon of the 1.5-km-deep Barranca de Rio Guerachi, west of Guachochi (Figs. 1 and 2); this section is exemplary of geology in the previously under-explored Parral-Guachochi-Guadalupe y Calvo region and the northern and central SMO generally (e.g., Swanson et al., 2006; Murray et al., 2013).

The base of the SMO is unconformable on Eocene (?) andesites and late Cretaceous (?) granite. In this section, the lower third of the SMO is composed of >500 m of silicic lavas and volcanic breccias, with minor debris flow deposits and sandstones. The middle ~650 m is composed of horizontal, ~40–120-m-thick ignimbrites separated by rare, thin, and altered fall deposits and paleosols. The lower units are weakly to moderately welded, lithic-rich lapilli tuffs. The two uppermost ignimbrites are very strongly welded with small fiamme and few lithic clasts; the informally named “San Miguel Arcangel” ignimbrite (Fig. 4) is weakly rheomorphic at its base and has a prominent basal vitrophyre. The uppermost ~350 m contains several welded dacitic ignimbrites, interbedded felsic volcaniclastic layers with paleosols, and ~150 m of felsic and andesitic breccias. The top ~50 m is a sequence of mafic lavas. To the best of our knowledge, this is the first detailed description and sampling of such a thick stratigraphic section in the SMO.

All central transect samples and a small number of archived samples from the northern and southern transects were available for whole-rock and petrographic analysis. Additionally, several archived samples have published data and descriptions available. All but five units are rhyolites, one of which is peralkaline (Acantilado tuff; CH97-7); the others are rhyodacites (3-71-2, SMO15_11, and _52), dacite (SMO15_10), and latite (SMO15_41). The rhyolites are typically 10%–50% phyric with ≤30% lithic fragments, and a typical crystal assemblage of plagioclase-quartz-biotite, ± alkali feldspar, ± hornblende. The groundmass is usually thoroughly devitrified to a cryptocrystalline mass of interlocking quartz and albite. Plagioclase, biotite, and hornblende are usually subhedral, broken crystals; many plagioclases are sieve-textured. Quartz content is highly variable within each sample, and all quartz crystals are anhedral and often embayed. Hydrothermal alteration is pervasive in many samples with sericitized feldspars, chlorite replacement of amphibole, and vesicles filled with chalcedony.

Zircon crystals are characteristically euhedral and oscillatory zoned and typically vary from 50–150 μm in length (e.g., Fig. 3). Anhedral, metamict cores are very rare. Melt inclusions are often abundant (>5 per crystal) and are typically subparallel to the crystal margins. No anhedral or rounded, strongly embayed, or metamict crystals were recovered.

U-Pb age data are presented in Figure 5. ²⁰⁶Pb/²³⁸U zircon crystallization ages (n = 578) were determined to constrain εHf_T values and mean square of weighted deviates (MSWD) for WMAs are high, ranging from 2–9 and 0.1–4 for LA-ICPMS and SIMS, respectively. Mean weighted ages (WMA) of zircon calculated in IsoplotR (Vermeesch, 2018) are presented in Table 1.

²⁰⁶Pb/²³⁸U zircon crystallization ages are unimodal within each sample and are overwhelmingly concordant (Fig. 5). Age-distinct antecrysts ~3–8 m.y. older than the mean age are common, usually two or three per sample. Xenocrysts are very uncommon, and only 20 grains (~3%) older than 56 Ma were analyzed in ten samples. Eighteen Cretaceous to Paleocene grains were identified in samples 3-20-4 (n = 3), CH01 (1), CH03 (3), J406 (2), MAJ (2), N-12 (2), RED (1), SJ-3 (1), SMO15_09 (1), _10 (1), and _14 (1). One Devonian xenocryst (ca. 392 Ma; SMO15_09) and one Triassic xenocryst (ca. 237 Ma; CH03) were measured; no Lower Paleozoic or Precambrian xenocrysts were measured.

Mean crystallization ages are within error of, or are up to 2 m.y. older than, the corresponding published K/Ar or ⁴⁰Ar/³⁹Ar cooling ages (Table 1) except for samples SMO15_17 and _44, which are 3 and 4 m.y. older, respectively, and 3-20-4, MAJ, SJ-3, SMO15_46 and _52, which are all significantly (2–4 m.y.) too young. These anomalous ages are nevertheless used to calculate the corresponding εHf_T value because of the insensitivity of the long half-life Lu-Hf system to 1–10 m.y. age variations.

Sierra Madre Occidental rhyolites have published mean whole-rock Hf contents of ~9 ± 4 ppm (n = 157), ranging from 2–20 ppm (http://portal.earthchem.org; Supplemental Material, Section C [^{footnote 1}]). Lu-Hf isotopes are presented as εHf_T values based on the corresponding ²⁰⁶Pb/²³⁸U crystallization age of the same analytical spot (e.g., Fig. 3), grouped by sample (Fig. 6); Supplemental Material, Section D). The highest single crystal analysis measured is +7.7 epsilon units from a 393 Ma xenocryst in SMO15_09, and the lowest is −16 epsilon units from a 39 Ma grain in SMO15_55. Most samples show less than 2.5 epsilon units of variation between phenocrysts, and 1σ uncertainties are often <1.0 epsilon units (Table 1). Including εHf_T and age outliers, 1σ uncertainties are usually less than 1.5 epsilon units, and the largest is 2.6 epsilon units in SMO15_52. The highest mean εHf_T value is +4.6 epsilon units ± 1.0 (N-99; 29.93 Ma), and the lowest is −4.1 epsilon units ± 1.4 (CH88-17; 37. 23 Ma), both in the Chihuahua transect (Table 1). The average of all εHf_T values is +1.5 epsilon units.

δ¹⁸O_zircon values range from 4.4‰ to 10.4‰ (Fig. 7); Supplemental Material, Section E [^{footnote 1}]), and the range within individual samples varies by <4.5‰ and commonly <2.5‰. There is no systematic variation in δ¹⁸O_zircon between cores and rims in the same grain (e.g., Fig. 3). The lowest δ¹⁸O_zircon single-grain values measured are 4.4‰–4.7 ‰ in SMO15_01, _41 (two grains), _46, and _52; these and an additional seven grains are all less than the 5.3‰ ± 0.3‰ (1σ) mantle zircon of Valley et al. (2005). The highest mean δ¹⁸O_zircon value (Table 1) is 8.0‰ ± 2.2‰ (SMO15_52), and the lowest is 5.5‰ ± 0.2‰ (CH03). Two anomalously high single-grain δ¹⁸O_zircon values measured 23.9‰ and 15.8‰, respectively, both from SMO15_44.

Eight δ¹⁸O_quartz measurements are available, six for samples with coexisting zircon (CH03, CH97-7, N-99, SJ-54, K-LP-T, and SL-41; Table 1), and two without (CH06 and 2-71-1). δ¹⁸O_quartz values range from 8.89‰–10.23‰ and are ~3.0‰–4.0‰ higher than for coexisting zircon. Oxygen isotope data from coexisting quartz-zircon pairs suggest subsolidus oxygen isotope fractionation in quartz between ~530–685 °C (Fig. 8).

Zircon saturation temperatures (T_Zircsat; Watson and Harrison, 1983; Hanchar and Watson, 2003), following the approach of Boehnke et al. (2013), were calculated from whole-rock analyses of 22 samples collected along the central transect and from three published analyses from the northern transect (Supplemental Material, Section F [^{footnote 1}]). Zirconium abundances for these and an additional 67 central transect samples (Moreno, 2016; Pack, 2017) are between ~75 and ~630 ppm and M-values ([Na + K + 2Ca]/Si × Al) are between 0.9 and 1.83 (Fig. 9A). Zircon saturation temperatures (T_Zircsat) range between 680 and 840 °C (Fig. 9B), with a median value of ~750 °C (Fig. 9C). Utilizing the model of Watson and Harrison (1983) results in temperatures ≥100 °C higher than that of Boehnke et al. (2013); taking this into account, our data are comparable (Fig. 9D) with those calculated for six SMO rhyolites by Bryan et al. (2008) and 12 from Siégel et al. (2018).

The estimated liquidus temperatures for SMO rhyolites were calculated using the rhyolite-MELTS and MELTS_Excel packages (Gualda et al., 2012; Gualda and Ghiorso, 2014). Liquidus temperatures were calculated for 25 rhyolites at 125 MPa (~5 km depth) and 0–6 wt% H₂O. A typical SMO rhyolite has a liquidus temperature of 944 °C with 1 wt% H₂O and 861 °C with 3 wt% H₂O (Fig. 9C).

Sierra Madre Occidental tuffs and lavas exhibit limited ranges of εHf_T and δ¹⁸O_zircon within individual samples and throughout the suite (Fig. 10). Based on these data, the SMO rhyolites have δ¹⁸O_zircon values ranging from crystallization under mantle conditions (~5.3 ± 0.3‰ (1σ); Valley et al., 2005) to 7‰–8‰, and εHf_T values from +6 to −6 epsilon units with a mean of ~+2 epsilon units. SMO rhyolites are isotopically homogeneous and form a distinct cluster in εHf_T-δ¹⁸O space. There is no systematic temporal (Fig. 10) or spatial variation in mean δ¹⁸O_zircon and εHf_T values, or T_Zircsat (Fig. 9D). This implies a uniform petrogenesis throughout the duration (~20 m.y.) and extent (4 × 10⁵ km²) of SMO silicic volcanism. Combining their homogeneous zircon isotopic character with published whole-rock elemental and isotopic data leads us to conclude that the SMO silicic magmas were remarkably monotonous.

The SMO likely buries a thick (25–30 km), Phanerozoic succession dominated by marine sedimentary and volcaniclastic lithologies (Fig. 11); Bryan et al., 2008). Magmatic assimilation of significant volumes of Phanerozoic sediment derived from Precambrian basement should have generated strongly negative εHf_T values. Assimilation of upper-crustal rocks, including pelites and carbonates (e.g., Bindeman, 2008), can raise the δ¹⁸Ozircon value significantly. In contrast, assimilation of high-temperature hydrothermal systems (e.g., intracaldera successions; Bindeman et al., 2007) and epithermal alteration zones (e.g., Boroughs et al., 2005) can significantly lower the δ¹⁸Ozircon value. Therefore, the absence of low- (<4‰; Valley et al., 2003) and high-δ¹⁸O_zircon (≥8‰), and strongly negative εHf_T zircons is inferred to indicate that SMO magmas did not assimilate large volumes of upper- and middle-crustal rocks, including carbonaceous marine sedimentary rocks exposed around Parral (Fig. 1).

The narrow ranges of δ¹⁸O_zircon and εHf_T values (Figs. 6, 7, and 10) are typical of magmas that have evolved, at least in part, by fractional crystallization of mafic melts (cf. Balsley and Gregory, 1998) in “hot” arc regimes with significant input from the mantle (e.g., Attia et al., 2020). Assimilation of some crustal material is required to explain the few xenocrysts identified and the rhyolites with weakly negative εHf_T values (Fig. 10). The paucity of xenocrystic zircon corroborates the results of Bryan et al. (2008) and suggests that magmatic temperatures (~700–850 °C) and xenocryst residence times were great enough to dissolve most older zircon. It is unlikely that the elevated magmatic temperatures needed to allow fractional crystallization processes and xenocryst destruction could be sustained for very long (10⁶ years) in the upper crust (e.g., Annen et al., 2006; de Silva and Gosnold, 2007) without cooling and crystallization (e.g., Miller et al., 2003).

Instead, eruptible SMO magmas evolved from melts in excess of their zircon saturation temperatures (>850–900 °C). The moderately crystal-poor and xenocryst-free rhyolites of the SMO were likely “hot” magmas (e.g., Miller et al., 2003; Halder et al., 2021) that required significant advective thermal input into the lower crust to sustain the necessary high magmatic temperatures. This anomalously high thermal flux was sustained throughout the duration and extent of SMO volcanism to account for the monotony of the rhyolites. Bryan et al. (2008) attributed the sustained high thermal flux to the repetitive emplacement of mafic bodies into the lower crust, thereby transferring both mantle melt and advected heat. The lower crust evolved into a region of melting-assimilation-storage-homogenization (Fig. 11); “MASH zone”; Hildreth and Moorbath, 1988), where melts of the country rock protolith and older mafic intrusions are mixed and assimilated into juvenile magmas undergoing fractional crystallization (e.g., Lee et al., 2007). Thus, AFC processes are inextricably linked. Thermal models demonstrate the efficiency of repetitive mafic intrusions to elevate the temperature of the surrounding lower crust (Huppert and Sparks, 1988; Annen and Sparks, 2002; Dufek and Bergantz, 2005; Annen et al., 2006). The same models indicate that temperature can reach >300 °C above the solidus (950–1000 °C) in 1–1. 5 m.y. (Annen et al., 2006) provided that the intrusion rate is 50 m (sill thickness) per ≥0.01–0.1 m.y. Super-liquidus silicic magmas formed in the lower or middle crust will have some concentration of dissolved water in them and will likely have viscosities in the 10⁴–10⁵ Pa.s range (e.g., Scarfe et al., 1987; Giordano et al., 2008), many orders of magnitude less than at the surface. Such low viscosities allow for efficient transport through the remaining crust from where the magmas are either erupted directly or assembled over short timescales (e.g., Annen, 2009) in shallow sills associated with calderas.

Xenolith localities and a very few outcrops around the northern SMO allow sampling of basement lithologies likely representative of the lower crust where SMO magmas are inferred to have been generated. Lower-crustal xenoliths from La Olivina (Fig. 1); Ruiz et al., 1988b; Rudnick and Cameron, 1991; Cameron et al., 1992; Nimz et al., 1995) include mafic granulites, intermediate to felsic orthogneisses, and paragneisses. Xenolith U/Pb zircon age populations are diverse and complex (Rudnick and Cameron, 1991), and all gneissic samples show evidence for granulite-facies metamorphism in the Neogene, Late Triassic, and at ca. 1.1 Ga (Grenville orogeny). La Olivina is placed south of the Laurentian craton margin in the Late Paleozoic Ouachita suture zone (Fig. 1) by Ortega-Gutiérrez et al. (2018). The presence of xenoliths with ≥1.0 Ga metamorphic ages requires that either (1) the basement beneath La Olivina is the very northernmost tip of Oaxaquia (ca. 1.2–0.9 Ga; Ortega-Gutiérrez et al., 2018), or (2) the Laurentian margin has to be extended at least 20 km farther southeast (Ruiz et al., 1988a; Rudnick and Cameron, 1991; Cameron et al., 1992).

The only mafic granulite xenolith dated has a mean age of 1.3 ± 0.3 Ma, within error of the K/Ar age of the surrounding alkali basalt lavas (Rudnick and Cameron, 1991). The mafic granulites are isotopically primitive compared to the gneisses (⁸⁷Sr/⁸⁶Sr₃₀ ~0.705; εNd₃₀ ~0) and indistinguishable from the SMO and Neogene mafic volcanic rocks (Cameron et al., 1992). They are typical of mafic lower-crustal rocks throughout northern Mexico (Ruiz et al., 1988b). A single granulite-facies, sillimanite-paragneiss has a maximum depositional age of 1.0 Ga and was metamorphosed at 40–20 Ma, broadly contemporaneous with SMO magmatism. The paragneisses have very strongly negative εNd₃₀ (−8 to −16 epsilon) and high ⁸⁷Sr/⁸⁶Sr₃₀ (>0.720). Differently aged but compositionally similar intermediate orthogneisses are present below La Olivina, where xenoliths with inferred magmatic ages of ca. 1.4 Ga, ca. 300 Ma, and 37–25 Ma have been recovered (Rudnick and Cameron, 1991). Rudnick and Cameron (1991) interpreted 37 −25 Ma sample MN-40 (64.5 wt% SiO₂) to be deep-seated relation of SMO rhyolites, and they have identical strontium, neodymium, and lead isotopes values (Cameron et al., 1992). The older orthogneisses are less primitive with strongly negative εNd₃₀ values and form a seamless trend with the paragneisses.

Xenoliths from Potrillo in northern Chihuahua and basement outcrops at Los Filtros north of Chihuahua City (Fig. 1) provide additional constraints on plausible lower crust beneath the northern SMO. Non-dated paragneiss xenoliths from Potrillo are less isotopically evolved (εNd₃₀ ~−7 epsilon; ⁸⁷Sr/⁸⁶Sr₃₀ ~0.706) than those at La Olivina (Ruiz et al., 1988b). High-grade orthogneisses and granitoids spanning 1.3–1.0 Ga at Los Filtros have initial ⁸⁷Sr/⁸⁶Sr values of ~0.708 and are inferred to be part of a belt of Grenville-age Laurentian basement extending to and beyond the Llano Uplift in central Texas (Ruiz et al., 1988b; Ortega-Gutiérrez et al., 2018).

Cameron et al. (1992) interpreted Cenozoic xenoliths at La Olivina as evidence for wholesale or significant modification of the lower crust by mafic underplating and associated granulite-facies metamorphism (Fig. 11). However, because there are no basement outcrops or xenoliths from within or adjacent to the central SMO between La Olivina and the Mesa Central (Fig. 1), this interpretation cannot be rigorously tested. The same area is underlain by the Mesozoic accreted Alisitos-Guerrero composite terrane (Fig. 1); Centeno-García et al., 1993); the composition and age of the lower crust of the Guerrero terrane and its subcontinental lithospheric mantle are unknown (Fig. 11). The lack of sub-SMO information contributes to a “chicken and egg” paradox where geochemical data yield insights into the buried basement age and composition, but the same data are modeled to understand potential sources of the SMO magmas and their relative contributions (e.g., Cameron and Cameron, 1985; Ruiz, et al., 1988b; Cameron et al., 1989; Wark, 1991).

The observed difference in εHf_T (~+14 epsilon) and δ¹⁸O (~+1.2‰) values in SMO zircons from mantle values, where 35 Ma depleted mantle is +15.5 epsilon units and 5.3‰, is best explained by AFC processes (e.g., DePaolo, 1981). Lower-crustal xenoliths from La Olivina and Potrillo provide potential end-members with which to model assimilation of different crustal contributions with both enriched and primitive mantle-derived melts. We applied Equation 15a of DePaolo (1981) to model AFC paths in εHf_T-δ¹⁸O space (Fig. 12) with the same fixed rate of fractional crystallization twice that of assimilation (r = 0.5) used by Cameron and Cameron (1985) and Cameron et al. (1989). We modeled a range of scenarios using three possible crustal end-members and four different mantle end-members (Supplemental Material, Section G [^{footnote 1}]). The crustal end-members are a mafic orthogneiss and a paragneiss from La Olivina (Cameron et al., 1992) and a paragneiss from Potrillo, northern Chihuahua (Ruiz et al., 1988b). Epsilon hafnium values for 35 Ma zircon (εHf₃₅) were calculated from published whole-rock εNd₃₀ values using the terrestrial array equation of Vervoort et al. (2011). Whole-rock δ¹⁸O values were calculated from δ¹⁸O_zircon using an equilibrium temperature of 800 °C (Bindeman et al., 2007).

Figure 12A depicts AFC mixing lines between enriched mantle sources and the three crustal end-members. Twenty percent assimilation of either paragneiss end-member with moderately enriched mantle (EMII₃₅) reproduces the δ¹⁸O_WR and εHf₃₅ values most characteristic of the SMO rhyolites. Furthermore, assimilation of ~10% paragneiss reproduces the isotopic characteristics of the coeval mafic granulite xenoliths inferred by Cameron et al. (1992) to be cumulates formed during fractional crystallization of the SMO silicic magmas. Negative epsilon hafnium and δ¹⁸O_WR ≈8 ‰ SMO rhyolites (Fig. 12A) are best modeled by significant (i.e., ≥50%) assimilation of intermediate composition orthogneiss. Although there are relatively few of these rhyolites (~4; Table 1), they all occur at the eastern ends of transects, on inferred, Proterozoic continental basement (i.e., not on the accreted Guerrero terrane). Furthermore, samples CH01 and CH88-17 of the northern transect are both less than 50 km west of exposed Proterozoic basement (1.3–1.1 Ga) at Los Filtros, Chihuahua (Fig. 1); Ruiz et al., 1988a, 1988b).

High epsilon hafnium (>+4) samples are modeled poorly by EMII₃₅ (Fig. 12A); so Figure 12B presents some alternative mantle end-members: depleted mantle (DM₃₅) and a putative δ¹⁸O-enriched (~9 ‰) form of EMII₃₅ based on Liu et al. (2014) representing hydrated arc mantle. The high epsilon hafnium samples are modeled successfully when using a mix of 10%–20% paragneiss or orthogneiss and a combination of EMII₃₅ and the δ¹⁸O-enriched arc mantle, or alternatively, a form of EMII₃₅ not as highly δ¹⁸O enriched as that described by Liu et al. (2014) from Tibet. Another explanation is satisfied if the EMII₃₅ mantle beneath the SMO was enriched in hafnium from subducted sediments; however, such sediments cannot have been derived from Precambrian crust without lowering the εHf₃₅ values, and must, therefore, have been relatively young (i.e., Phanerozoic protoliths). Finally, depleted mantle is not a viable end-member for SMO silicic magmas (Fig. 12B) due largely to its very low hafnium concentration; however, it is potentially a source of melt for the coeval southern Cordilleran basaltic andesite (SCORBA) mafic lavas with ~2%–3% paragneiss assimilated.

Our data suggest a juvenile-recycled crust ratio of ~4:1 for the SMO. Assuming an erupted volume of ~3.9 × 10⁵ km³ (Bryan et al., 2008) implies ~3.1 × 10⁵ km³ of silicic magma fractionated from a mafic melt in the lower crust. Assuming this was 100% of the available residual silicic melt fraction after 90% crystallization (e.g., Dunbar et al., 1993; Bower and Woods, 1998) implies fractionation from a mafic source with a volume of ~2.8 × 10⁶ km³. These values are underestimates by ~3–10 times (e.g., White et al., 2006) if inferred upper-crustal SMO reservoirs (Fig. 10) are included.

Estimating the magmatic and eruptive fluxes over the ~20 m.y. duration of apparently uninterrupted SMO volcanism (McDowell, 2007; McDowell and McIntosh, 2012) yields values of ~0.14 km³ yr⁻¹ and ~0.02 km³ yr⁻¹, respectively. White et al. (2006) calculated a volcanic flux of 0.03 km³ yr⁻¹ for the SMO over a 15 m.y. time interval. Time-averaging underestimates the peak volcanic fluxes of discrete pulses. Seventy percent of SMO eruptive ages reported in McDowell and McIntosh (2012) are in a 7 m.y. range between 36–29 Ma. If 70% of the volume of the SMO erupted at this time, the peak magmatic and volcanic fluxes are 0.28 km³ yr⁻¹ and 0.03 km³ yr⁻¹, respectively.

Mean global magma production since the Jurassic is ~26–34 km³ yr⁻¹ (Crisp, 1984); magmatism associated with the SMO was ~1%–1.5% of this long-term average. Mean global volcanic output is ~0.01–0.1 km³ yr⁻¹ (Crisp, 1984; White et al., 2006). Silicic and intermediate volcanism both average 0.004 km³ yr⁻¹, basaltic volcanism (excluding flood basalts) averages 0.026 km³ yr⁻¹, and flood basalts average 1 km³ yr⁻¹ (White et al., 2006). Therefore, the volcanic part of the SMO had a similar output to the long-term global mafic output and was ten times that of intermediate or silicic output. When Quaternary silicic and intermediate volcanic outputs are combined, similar ~10⁻² km³ yr⁻¹ rates are estimated at the Ethiopian and Kenyan portions of the East African Rift, the Kamchatkan, Kyushu, and Taupo intra-arc rifts (White et al., 2006).

Petrogenesis. We estimate the juvenile component of SMO magmas to be ~80% (Fig. 12). This is in excellent agreement with several AFC models based on trace-element and Sr and Nd isotopic analyses (e.g., Lanphere et al., 1980; Cameron and Cameron, 1985; Cameron et al., 1986; Wark, 1991; Smith et al., 1996). Fractional crystallization of a primitive melt to produce a silicic magma, enhanced further by partial melting of the surrounding lower crust, also provides a mechanism to advect-in the necessary thermal input to sustain and drive the AFC process (e.g., Annen and Sparks, 2002). This would be enhanced further by latent heat of crystallization from the fractionating mafic cumulate.

Our data are incompatible with the SMO magmas forming by assimilation of >50% continental crust (e.g., Ruiz et al., 1988b; Bryan et al., 2008) except along the eastern margins of the SMO where there is thick Proterozoic crust and an absence of accreted Mesozoic terranes (Fig. 12A). This implies that voluminous mafic melt and heat were transferred from the mantle into the base of the crust (cf. Ducea and Barton, 2007) causing melting, assimilation, and fractional crystallization under the thinner, younger crust below the western and central SMO. Thermal models indicate that conduction of thermal energy from the mantle is insufficient to drive and sustain melting and assimilation of continental crust. However, it is difficult for hafnium and oxygen compositions to discern fractionation of juvenile mafic magma and re-melting and subsequent fractionation of a slightly older mafic body (e.g., Bryan et al., 2008), and while neither are mutually exclusive, we cannot rule the latter out entirely. In contrast, in the thicker, older crust below the eastern SMO, heat transfer was dominant and drove anatectic melting of fertile Proterozoic orthogneisses and paragneisses with less mixing-in of mafic melts.

Crustal Thickness. Any AFC model needs to accommodate the addition of enormous volumes of mafic melt in the lower crust, most likely as a mafic underplate. Ruiz et al. (1990)'s most compelling argument against contemporary AFC models was the strong negative gravity anomaly co-spatial with the SMO. They inferred the presence of a low density, granitoid batholith in the upper and middle crust beneath the SMO. They discounted the alternative of a 4-km-thick mafic underplate (Cameron et al., 1980) in the lower crust as this would be expected to produce a positive gravity anomaly. However, there is no supporting evidence for an upper-crustal batholith; for example, there are no Cenozoic granitoid xenoliths reported. Furthermore, our hafnium and oxygen isotope data are inconsistent with assimilation of large volumes of middle- and upper-crustal rocks (e.g., δ¹⁸O > 10‰).

This inconsistency can be reconciled if some or all of the mafic underplate has been removed and recycled into the mantle (e.g., Lee et al., 2011). Assuming that the underplate had the same area as the modern extent of the SMO, the underplate would be >7 km thick; in contrast, reconstruction of the ~40-km-thick Sierra Nevada batholith crustal section requires a now missing ~80 km of mafic and ultramafic arc root (Saleeby et al., 2003). The thickness of the crust beneath the SMO is not well constrained; estimates range from ~55 km in the north (Bonner and Herrin, 1999) and ~40 km in the central SMO (Gomberg et al., 1989), to a uniform 30–35 km (Laske et al., 2013; Zhu et al., 2016). The minimum thickness of the necessary mafic underplate, therefore, is ~7%–23% of the current crust.

The absence of pyroxenite, ultramafic cumulate, or mantle xenoliths at La Olivina is consistent with removal of a mafic underplate and/or arc root before ca. 1.5 Ma. Application of a simple Airy isostasy model to the ~2000 m asl SMO plateau is compensated by a missing 10-km-thick mafic underplate (2.86 g/m³) or a 5-km-thick underplate of equal parts mafic and pyroxenite (3.13 g/m³). Removal of either of these masses would produce a negative gravity anomaly and anomalously high topography. Thicker SMO crust compensates for high topography but not for the negative gravity anomaly.

Extension and thinning of the upper crust was broadly contemporaneous and co-spatial with Cordilleran volcanism during retreat and rollback of the Farallon slab (e.g., Damon et al., 1981; Ferrari et al., 2007). Although diffuse across the Basin and Range province, extension became localized to either margin of the SMO in Mexico, where the Paleocene–Eocene arc had not advanced so far to the east, and therefore, had less far to retreat. There was probably insufficient slab rollback to expose a large enough region of already hydrated subcontinental lithospheric mantle to allow enough up-welling of isotopically primitive asthenosphere (i.e., depleted mantle). The strongest argument in support of this is that rollback was much greater under the Basin and Range, and yet the volumes and rates of volcanism are an order of magnitude smaller than in the SMO (e.g., White et al., 2006), regardless of whether or not AFC processes were dominant.

Bryan and Ferrari (2013) and Ferrari et al. (2018) recognized the bimodal nature of calc-alkaline volcanism (basaltic andesite [SCORBA] and rhyolite) in the SMO region through the Oligocene and early Miocene. They attribute this to the gradual waning of the arc as a slab tear that widened to form a slab window was subducted beneath northwestern Mexico (Ferrari et al., 2018; Mahar et al., 2019). As a result, the local mantle potential temperature increased as up-welling asthenosphere was exposed through the slab window. Reorientation of far-field plate motions (Ferrari et al., 2018) and the local addition of heat to weaken the lithosphere gradually focused crustal extension along the western, younger margin of the SMO where the ca. 24 Ma flare-up was localized. Ultimately crustal extension over the slab window led to the opening of the active spreading center in the Gulf of California (6–0 Ma).

The onset of these gradual changes to crustal extension and increased heat from the deeper mantle began likely before the ca. 24 Ma flare-up within the SMO and are thought to have been on-going throughout much, if not all, of the Oligocene. We propose that the onset of extension in the lower crust allowed mafic melts to more readily intrude than they had in the Eocene and before. However, compared to modern continental arcs, the Oligocene SMO had a volcanic flux about double; therefore, it is likely that extra factors were at play. We suggest that the otherwise unexceptional increase in mantle melt supply was amplified by (1) the remnants of the lower crust and mantle wedge being hydrated and isotopically enriched through the late Cretaceous to Eocene, and (2) the arrival of the still hot and thin oceanic lithosphere preceding the opening slab tear and/or window. Asthenospheric melts from the slab window are a potential source of the coeval SCORBA mafic lavas (Fig. 12B). The presence of readily melted, extensively hydrated lithosphere is central to interpretations of enhanced magmatism through slab rollback (e.g., DeCelles et al., 2009), and cannot be ruled out in this scenario. Slab break-off and slab windows are widely recognized as major perturbations to arc systems that enhance magmatism. These processes are not mutually exclusive, and their coincidence can explain the anomalously high volume of the SMO. If correct, the origins of the SMO are distinctly different from the Mesozoic batholiths of the Cordillera (e.g., Ducea and Barton, 2007), even though they are of similar scales.

Most SMO magmas were generated by juvenile melts that assimilated ~20% lower crust. The mantle-derived melts fractionated in the lower crust and advected heat into surrounding orthogneisses and paragneisses. The SMO magmas are transitional in Hf-O isotope space between continental arc granites and those formed in extensional arcs and continental rifts, in agreement with interpretations where the SMO represents the first stage of magmatism over subducting slab window during the Eocene to Miocene evolution of the Cordilleran margin from convergent to extensional. The immense volume of the SMO requires an anomalously strong heat source such as a slab window coupled with extension of the upper plate to allow for transfer of magma and heat.

We recognize that this research took place throughout the historic lands of the Rarámuri, and we thank them for their hospitality. Andrews and Busby were supported by National Science Foundation (NSF) award EAR-1019559. Andrews received additional support from NSF awards EAR-1137774 and EAR-1547784. WiscSIMS is supported by NSF (EAR-2004618) and the University of Wisconsin-Madison. The ion microprobe facility at University of California, Los Angeles (UCLA) is partly supported by a grant from the Instrumentation and Facilities Program, Division of Earth Sciences, NSF. Eric Swanson kindly donated several samples. We thank John Valley and Mike Spicuzza at UW-Madison for assistance in oxygen isotope analysis by SIMS and laser fluorination. We thank Luca Ferrari, Jose Duque, Nick Moreno, Brenda Pack, Axel Schmitt, Ming-Chang Liu, Mark Pecha, Matthew Wielicki, Ryan Ickert, Mattias Smit, Charles Knapp, and Scott Boroughs, and we are indebted to William Leeman for support and encouragement. We thank two anonymous reviewers, Associate Editor Jonathan Miller and Editor Shanaka de Silva, for their constructive comments.