Petrology and geochronology of Cretaceous–Eocene plutonic rocks in northeastern Washington, USA: Crustal thickening, slab rollback, and origin of the Challis episode Open Access

Spatial and/or temporal variations in the composition, geographic distribution, and production rate of magmas are recorded in the rocks of both modern and ancient continental arcs. These variations have been linked to a number of factors, including changes in the age, geometry, or composition of the subducting plate, differences in the nature of the mantle wedge and overlying crust, and shifts in the regional tectonic environment (Ducea et al., 2015). Examples of many such variations are seen in rocks of the ca. 45 Ma to recent Cascade arc, where spatial geochemical patterns have been linked to differences in slab age (e.g., Green and Sinha, 2005) or geometry (e.g., Guffanti and Weaver, 1988), heterogeneities in mantle composition (e.g., Leeman, 2020), or variable inputs of slab-derived fluids (Leeman et al., 2004; Pitcher and Kent, 2019), and decreasing magma production over time, which has been attributed to a slowing subduction rate (Verplanck and Duncan, 1987).

However, the most profound change in arc magmatism in the Pacific Northwest occurred immediately prior to the establishment of the modern Cascade arc. Its cause(s) remains enigmatic. Between ca. 52 and 45 Ma, during an interval known as the Challis episode (Armstrong, 1974), the distribution of igneous activity transitioned from an ~800-km-wide belt extending at least as far east as Montana (Fig. 1) into the narrow N-S–trending present-day arc. Rocks associated with the Challis episode have chemical and mineralogical traits of subduction-related systems, but their wide distribution and association with orogenic collapse and core complex formation are atypical of an arc setting, and thus their origin has been a long-standing topic of debate (e.g., Ewing, 1980; Morris et al., 2000; Madsen et al., 2006). A goal of this study was to combine U-Pb geochronology, trace-element data, and Sr-Nd-O isotopic compositions with petrology to elucidate the cause(s) of the Challis episode.

A related challenge is that rocks with arc-like petrologic traits are not necessarily the product of contemporaneous subduction: Those geochemical signatures can be inherited from mantle or crustal sources affected by earlier subduction events (Hooper et al., 1995; Morris et al., 2000). The Pacific Northwest, with its long history of subduction-related magmatism, provides an excellent opportunity to assess whether it is possible to distinguish granitoid rocks generated by contemporaneous subduction from those produced by remelting of preexisting arc crust. We approached this question by comparing the geochemistry of three groups of plutonic rocks in Washington—two that are associated with contemporaneous subduction (Cretaceous–Paleocene intrusions in northeastern Washington and those of the modern Cascade arc) versus Eocene plutons emplaced during the Challis episode.

The study area encompasses ~41,000 km² in northeastern Washington bounded to the south by the Miocene Columbia River Basalts and to the west by rocks of the Cascade arc (Fig. 2). Cretaceous plutonic rocks, some now orthogneisses, make up ~11% of the field area (~4500 km²), Paleocene–Eocene intrusions constitute another ~11% (~4400 km²), and Eocene volcanics make up ~3% (~1300 km²). The intrusions were emplaced into Proterozoic- to Mesozoic-age basement rocks that extend from the edge of the ancient North American craton across a series of accreted terranes. From east to west, Lund et al. (2015) recognized five basement domains within the study area: (1) the 2670–2650 Ma Pend Oreille sialic protocontinent, (2) the 780–500 Ma Kootenay passive-margin sediments, (3) the 320–250 Ma Okanagan oceanic arc rocks, (4) the 385–175 Ma Western Quesnellia oceanic arc rocks, and (5) material associated with mid-Cretaceous accretion of the Insular terrane, an assemblage of island-arc rocks that includes the Wrangellia terrane. The Okanagan and Quesnellia rocks are part of the composite Intermontane terrane, accreted during the Jurassic. These terrane-accretion events caused crustal thickening and development of an orogenic plateau, which in northeastern Washington is estimated to have been up to 2 km higher than present-day elevations (Wolfe et al., 1998) prior to its collapse in the Eocene. Assuming average densities for lower-crustal and upper-mantle rocks, this elevation difference implies 10–15 km of total crustal thickening, which is consistent with a total crustal thickness >60 km as inferred from petrologic data (Whitney et al., 2004).

Cretaceous arc magmatism associated with the arrival of the Insular terrane began ca. 115 Ma, producing a belt of batholiths that extends from the Canadian Coast Ranges, across the North Cascades and northeastern Washington, and into Idaho. Plutonism in this belt occurred in pulses with peaks at 114–102 Ma and 85–70 Ma in the southern Coast Ranges (Cecil et al., 2021) and at 96–89 Ma and 78–71 Ma in the North Cascades (Miller et al., 2009).

The ca. 52–45 Ma Challis episode began following a ca. 60–53 Ma magmatic lull that is attributed to flat-slab subduction of the Farallon (and Kula) plate (Schellart, 2020). Igneous rocks that formed during this episode define the Challis-Kamloops belt, which extends over 2000 km from central British Columbia to central Idaho (Fig. 1). In Canada, this belt is narrow (<200 km wide) and encompasses six discrete volcanic centers, all active ca. 55–45 Ma, which consist predominantly of medium- and high-K calc-alkaline lavas (Dostal and Jutras, 2022). South of the international border, the belt widens dramatically to >800 km and is farther from the continental margin. There, it comprises four major volcanic fields, including the Colville igneous complex in Washington (Morris et al., 2000), the Challis volcanics and associated intrusive rocks in Idaho (Fisher et al., 1992), and the Absaroka and Montana alkalic provinces in Wyoming and Montana (Sundell, 1993). The rocks are mainly high-K calc-alkaline basalts to rhyolites, but in Idaho, Wyoming, and Montana, shoshonitic and high-K alkaline rocks are also present. In northeast Washington, broad exposures of Challis-age calc-alkaline plutons (the focus of this study) are coeval with volcanic rocks preserved in grabens (Fig. 2).

Armstrong (1978) recognized that Challis magmatism was accompanied by widespread resetting of K-Ar ages and formation of extensional structures, including core complexes, grabens, and dike swarms. Core complexes in northeastern Washington include the Okanagan and Kettle Domes and the Priest River Complex (Fig. 2). The Okanagan Dome records melting at midcrustal depths between 61 and 49 Ma, followed by rapid exhumation and cooling, which was complete by 47 Ma (Kruckenberg et al., 2008). Rocks of the Priest River Complex record a similar history, with peak metamorphism 74–54 Ma, followed by rapid uplift, although the timing of uplift (54–44 Ma) is less tightly constrained (Stevens et al., 2015). Core complex exhumation was coeval with widespread dike emplacement that represents a minimum of 12 km of extension across northeastern Washington (Miller et al., 2022). The dikes are predominantly NE-striking, consistent with a regional strain field that also resulted in >300 km of dextral offset along the Straight Creek, Ross Lake, Leavenworth, and Entiat faults during this interval (Miller, 2022, personal commun.).

In summary, the Challis episode was associated with the collapse of an orogenic plateau that had formed in response to Cretaceous terrane accretion. In northeastern Washington, this collapse resulted in exhumation, extension, basin subsidence, and “arc-like” magmatism distributed over a broad area. This atypical association has led to diverse tectonic explanations, including collapse of a shallowly dipping slab (Ewing, 1980), decompression melting during extension (Morris et al., 2000), passage of a slab window (Thorkelson and Taylor, 1989; Madsen et al., 2006), intra-arc rifting (Dostal et al., 2005), slab breakoff and/or rollback (Schmandt and Humphreys, 2011; Dostal and Jutras, 2022), and lower-crust delamination (Bao et al., 2014).

Rocks included in this study are divided into four groups: (1) Cretaceous and Paleocene plutons, (2) Cretaceous and Paleocene dikes, (3) Eocene plutons, and (4) Eocene dikes. In this section, we summarize the field and lithologic traits of each group. Pluton ages are based on our U-Pb dating or on published mapping and geochronology. Based on previous mapping and dating (Miller et al., 2022), dike ages are assumed to be Eocene, aside from two we measured, and excluding a few coarser-grained bodies as described below.

This group includes samples from 24 map units older than 52 Ma (Table 1). Rock types range from hornblende diorite to two mica granite, but the most common lithologies are granodiorite, granite, and quartz monzonite. Grain size ranges from fine to coarse; K-feldspar phenocrysts are common and up to 5 cm in length. Color indices are typically low (2–15) with biotite representing the dominant mafic mineral; hornblende occurs in only eight of the plutons and generally only in minor amounts. Primary muscovite occurs in some Cretaceous and Paleocene plutons, as does, less commonly, magmatic epidote; both are indications of crystallization at pressures ≥4 kbar (Anderson, 1996). Other accessory minerals include sphene, apatite, zircon, allanite, opaque oxides, and (mainly in associated aplite dikes) garnet.

Although Eocene dikes are known to be widespread in northeastern Washington, one coarser-grained dike in the Duncan Hill domain (sample 30414-25) yielded a zircon U-Pb age of 89.3 Ma (see below), suggesting that other coarser-grained dikes may also be Cretaceous in age. Consequently, we have assumed that four coarse-grained foliated dikes near the margins of the Okanagan core complex are older than Eocene. These dikes are also wider (3–20 + m), more heterogeneous at the outcrop scale, and more felsic than most other dikes, consisting dominantly of plagioclase (up to 1 cm long), K-feldspar, and hornblende. Based on their texture and proximity to major faults, these dikes may be ductile shear zones rather than intrusions.

Samples from 13 intrusions younger than 52 Ma (Table 1) range from quartz diorite to granite, with quartz monzonites, monzodiorites, and granodiorites most prevalent. These plutons differ from the Cretaceous and Paleocene plutons in having higher color indices (7–40) and abundant hornblende, in some cases with pyroxene cores. Peraluminous minerals (muscovite, garnet) are absent, and K-feldspar phenocrysts are rare and smaller. Sphene, apatite, zircon, and opaque oxides are ubiquitous accessory phases; allanite and magmatic epidote are absent. Field evidence (miarolitic cavities, andalusite in contact aureoles, chilled border zones, intrusion into hypabyssal dikes) indicates the Eocene plutons were emplaced at shallow depth.

Eocene-aged dikes are widespread throughout northeastern Washington, but the majority are found within one of several clusters (Fig. 2), including those occurring in: (1) the Republic and Kettle Falls grabens (Holder et al., 1990), (2) Corbaley Canyon, (3) the Teanaway swarm, and (4) domains that accompany batholiths in the North Cascades, particularly the Cooper Mountain, Duncan Hill, and Golden Horn batholiths (Miller et al., 2022). Zircon U-Pb ages for dikes in clusters 2–4 range from 49.6 to 46.2 Ma (Miller et al., 2022; this study).

Most dike samples included in this study were collected from road cuts in Corbaley Canyon and in the southern and eastern domains adjacent to the Duncan Hill batholith (Fig. 2). One Corbaley Canyon dike had a zircon U-Pb age of 46.2 ± 0.9 Ma (this study), in good agreement with earlier K-Ar ages of 48.4 ± 2.2 and 47.8 ± 1.9 Ma from this swarm (Miller et al., 2022). These dikes are aphanitic or sparsely porphyritic (phenocrysts of feldspar, amphibole, and/or quartz), with 7–10 m average widths, WNW-ESE strikes, and steep dips. Compositions range from basalt to rhyolite, and some dikes are composite. Dikes associated with the Duncan Hill swarm are also generally aphanitic or porphyritic with phenocrysts of plagioclase or hornblende. They typically strike NE-SW, have steep dips, and are <4 m wide. These dikes range from basalt to rhyolite, but mafic and intermediate compositions are more typical. One of the hornblende andesite dikes (sample PL15-005) yielded a zircon U-Pb age of 55.7 ± 1.0 Ma (this study).

Over five field seasons (2005–2015), we collected a total of 224 samples, primarily from road outcrops. For whole-rock analysis, 144 of these samples, chosen on the basis of freshness and location/map unit, were trimmed or crushed and handpicked to remove any weathered material and then pulverized in an alumina or tungsten carbide shatterbox. Major elements were measured by inductively coupled plasma–optical emission spectrometry (ICP-OES) following lithium metaborate fusion, either at the University of Puget Sound (UPS) or at ALS Global Geochemical Laboratory in North Vancouver, British Columbia, Canada. Samples analyzed at UPS were run against a suite of natural rock powder standards, including BHVO-1, BIR-1, BCR-1, AGV-1, RGM-1, and G-2. Loss-on-ignition (LOI) was determined gravimetrically after firing at 1000 °C for 1 h. Trace elements, except for B, Be, and Li, were determined by inductively coupled plasma–mass spectrometry (ICP-MS) at the GeoAnalytical Laboratory, Washington State University (prefix CC05 and EI07 samples) or at ALS, again following lithium metaborate fusion. B, Be, and Li were analyzed by ICP-OES at UPS using the HF-mannitol method (Nakamura et al., 1992) and dedicated glassware (spray chamber, tubing, etc.) used only for light element analyses. For analyses at UPS, based on replicate in-run analyses of standards, analytical error was <0.6 wt% for SiO₂, <0.2 wt% for all other major oxides except K₂O (<0.3 wt%), and <0.02 wt% for MnO₂ and P₂O₅, while ALS reported analytical error <1.5 wt% for SiO₂ and <0.2 wt% for all other major oxides except Al₂O₃ (<0.4 wt%). From both laboratories, analyses with unnormalized totals ≥96.0% and ≤103.0% were considered acceptable. Analytical uncertainty was <10% relative for B, Be, and Li, based on replicate in-run analyses of natural rock standards AC-E, JA-2, JG-1a, JR-1, and SCO-1.

Sr and Nd isotopic analyses were performed at the University of New Mexico (UNM) on a Thermo Scientific Neptune multicollector (MC) ICP-MS and at the University of Washington (UW) using a Nu instruments MC-ICP-MS. At UNM, column separations employed Eichrom Sr-Spec resin for Sr and Eichrom TRU-Spec resin followed by LN-Spec resin for Nd, following the method of Mitchell and Asmerom (2011). Separation procedures at UW followed Nelson (1995); analytical procedures have been described by Brach-Papa et al. (2009) for Sr and by Gaffney et al. (2007) for Nd. Sr and Nd isotope ratios were corrected for mass fractionation to ⁸⁶Sr/⁸⁸Sr = 0.1194 and ¹⁴⁶Nd/¹⁴⁴Nd = 0.7219, respectively. Sr isotopic compositions were normalized to NIST987 (⁸⁷Sr/⁸⁶Sr = 0.710240) and Nd results to La Jolla Nd (¹⁴³Nd/¹⁴⁴Nd = 0.511843).

Oxygen isotopes were determined on quartz separates from 12 samples by laser fluorination at the University of Oregon. Picked clean quartz crystals (2–8 individual grains, 0.9–1.7 mg) were dissolved in purified BrF₅ reagent under a 35 W CO₂ laser. Gases were purified using liquid N₂ cryogenic traps and a Hg diffusion pump before conversion to CO₂ with a platinum-graphite filament. Analyses were on a MAT 254 10 kV gas source isotope ratio mass spectrometer, and measured ratios were corrected to repeat analyses of Gore Mountain garnet (UOG, accepted δ¹⁸O = 6.52‰). UOG analyses were within 0.08‰ 2σ (Loewen and Bindeman, 2015). All sample locations and elemental and isotopic data are provided in Supplemental File S1¹ and are available from the EarthChem Library (Tepper et al., 2023).

U-Pb zircon ages for 33 pluton and dike samples were measured by laser ablation (LA) MC-ICP-MS at the Arizona LaserChron Center, University of Arizona. Samples were chosen to provide broad geographic coverage of the study area and thereby facilitate recognition of spatial trends over time. Zircons were obtained by heavy liquid density and magnetic separation procedures, mounted in epoxy, polished to reveal grain interiors, and then imaged with cathodoluminescence (CL) to identify spots for analysis. U-Th-Pb isotopic data were measured in two sessions, the first (2008) using a GVI Isoprobe connected to a New Wave UP193 HE laser, and the second (2015) using a Nu ICP-MS connected to a Photon Machines Analyte G2 excimer laser (Gehrels and Pecha, 2014). Ages reported in this paper are the weighted means of 4–39 individual points per sample calculated using AgeCalcML (Sundell et al., 2021) and incorporating only points that were concordant within 2σ. Age uncertainties, reported with their accompanying mean squared weighted deviations (MSWD), incorporate both the uncertainty arising from scatter and precision in individual sample analyses and the systematic uncertainty associated with analyses of the standards. Additional method details and all geochronologic data are provided in Supplemental Files S2 (ranked order ²⁰⁶Pb/²³⁸U age plots and concordia diagrams) and S3 (method details and all individual point analyses).

Previous K-Ar dating of igneous rocks in northeastern Washington established that Eocene magmatism was widespread (e.g., Pearson and Obradovich, 1977; Armstrong, 1978), but those dates showed evidence of resetting and lacked the precision needed to discern spatiotemporal patterns. To address this shortcoming, we measured U-Pb zircon ages, most with 2σ uncertainties <1 m.y., on 33 samples collected throughout the study area (Table 2; Supplemental Files S2 and S3). Eight of the Cretaceous and Paleocene samples (40%) and two of the Eocene samples (15%) contained grains with cores up to 1730 Ma that were clearly inherited. More common (present in 70% of Cretaceous and Paleocene and 50% of Eocene samples) were grains marginally older (<10 m.y.) than the rest (Fig. 3) that we interpreted as antecrysts, cannibalized from slightly older but likely related material. All antecryst points, points that were discordant (at 2σ level), and points with anomalously young ages that we attributed to Pb loss were excluded from the age calculations. In general, this had minimal impact (<0.7 m.y. difference) on the final age calculations because the total number of grains dated (average 27 per sample) was large. One exception was dike PL15-005, in which 8/17 grains were clearly inherited, and the remaining nine points clustered in two groups, seven that yielded an age of 55.75 ± 1.04 Ma (MSWD = 0.98) and two that yielded 49.57 ± 1.59 Ma (MSWD = 0.10). The younger age is closer to that of another dike associated with the Duncan Hill batholith (47.5 Ma; Miller et al., 2022), but because it was based on two points that were U-rich and thus susceptible to Pb loss (Mattinson, 2005), we favor the older date as the emplacement age. Two samples (EI07-25, EI07-47) gave weighted mean ages with MSWD >4 (Table 2), indicating greater scatter than can be attributed to analytical uncertainty (MSWD ~1). Both of these samples contained a high proportion of inherited grains that were excluded from the mean age determination, but their high MSWD values suggest some inheritance remained.

The new ages ranged from 113 to 46 Ma (Fig. 4), but for purposes of discussion, we divided them into Cretaceous and Paleocene (>52 Ma) and Eocene (<52 Ma) groups. Samples in the Cretaceous and Paleocene group were widespread but showed no spatiotemporal trends within the study area (Fig. 5A). Among Cretaceous and Paleocene plutons, three age clusters were recognizable in a kernel density estimation (KDE) plot (Fig. 4; Wickham, 2016): (1) 113–110 Ma, (2) 92–87 Ma, and (3) 80–77 Ma. The 92–87 Ma group correlates in time with the 96–89 Ma magmatic flare-up in the North Cascades (Miller et al., 2009) as well early phases of activity (98–85 Ma) in the Idaho batholith (Gaschnig et al., 2010). Similarly, the 80–77 Ma episode overlaps in time with a second flare-up in the North Cascades (78–71 Ma; Miller et al., 2009) and with emplacement of the voluminous Atlanta lobe of the Idaho batholith (ca. 80–67 Ma; Gaschnig et al., 2010).

Eocene samples were also widely distributed across the study area, and their ages become progressively younger to the southwest, ranging from 52 Ma in northern Idaho to ca. 46 Ma in central Washington (Fig. 5B). Rocks of similar age occur in southern British Columbia (Ickert et al., 2009; Dostal and Jutras, 2022) and in Idaho (Gaschnig et al., 2010), but in neither case is any systematic age progression discernible.

Major-oxide contents of 114 Cretaceous–Paleocene and Eocene samples followed the trends expected for calc-alkaline suites: As SiO₂ increases, MgO, Fe₂O₃^t, CaO, and TiO₂ decrease, Na₂O and K₂O rise, and Al₂O₃ defines a hump pattern (Fig. 6). Average major-element contents of Eocene dikes and plutons were very similar, but the dikes (47–79 wt% SiO₂, 0.01–10.5 wt% MgO) showed a wider range of compositions (Table 3; Fig. 6). In contrast, with the exception of the gabbroic Trapper Peak stock, most Eocene pluton analyses fell between 60 and 72 wt% SiO₂ and 0.5–4.0 wt% MgO. The Cretaceous and Paleocene rocks, both plutons and dikes, were on average more felsic than the Eocene rocks (70 vs. 65 wt% SiO₂; 1.0 vs. 2.9 wt% MgO for plutons; Table 3), and at equivalent silica contents, they had slightly higher CaO and Fe₂O₃^t, but otherwise the Eocene and Cretaceous–Paleocene rocks displayed overlapping trends on major-oxide variation diagrams (Fig. 6). All samples were calc-alkaline with the exception of eight (11%) of the Eocene dikes (Fig. 7). Roughly 30% of the samples were classified as peraluminous with alumina saturation indices (ASI) ≥1.04; the remainder of the samples were metaluminous. ASI increases with SiO₂, and thus peraluminous compositions occurred in a higher proportion of Cretaceous and Paleocene (48%) than Eocene (22%) samples, but there was no evidence that ASI was influenced by location.

Bulk chemical distinctions between rocks of the two age groups were more apparent in trace-element data, where the Eocene samples tended to have higher Sr (Fig. 8A) and lower Yb at equivalent SiO₂ (Fig. 8B). Over 50% the Eocene samples were classified as adakites (Sr > 400 ppm, Sr/Y > 20, Yb < 1.8 ppm, La/Yb_N > 10; Figs. 8C and 8D; Defant and Drummond, 1990), and nearly 80% had at least three of these traits, which are also those identified for magmas associated with slab breakoff (Hildebrand and Whalen, 2014). With increasing SiO₂, the Cretaceous and Paleocene plutons showed trends of decreasing Eu/Eu* (indicative of feldspar fractionation) and decreasing Dy/Yb_N (indicative of amphibole fractionation) that were not seen among the Eocene plutons (Figs. 8E and 8F). All samples were moderately light rare earth element (LREE) enriched (Fig. 9), and on spidergrams, all samples displayed large ion lithophile element (LILE) enrichments and high field strength element (HFSE) depletions characteristic of subduction-related magmas (Fig. 10); this was also reflected in their high Ba/Nb values (Table 3).

Initial Sr-Nd isotopic ratios of 35 samples varied greatly (Table 4; Figs.11A and 11B), trending toward values indicating significant incorporation of older continental crustal material. Isotopic compositions of Cretaceous and Paleocene plutons (⁸⁷Sr/⁸⁶Sr_i = 0.7032–0.7108; ε_Ndi = −13.9 to +4.6) overlapped with those of the Eocene plutons (⁸⁷Sr/⁸⁶Sr_i = 0.7044–0.7132; ε_Ndi = −18.5 to +2.7), and, in general, the plutonic rocks became more isotopically diverse (extending to more crust-like values) with increasing SiO₂ or decreasing Mg number (Fig. 11C). Conversely, the Eocene dikes, despite having the broadest range of SiO₂ contents (Fig. 11C), showed the least isotopic variation (⁸⁷Sr/⁸⁶Sr_i = 0.7036–0.7053; ε_Nd = −1.2 to +4.4) and no correlation with degree of differentiation. For plutons, there was a clear correlation between isotopic composition and location, with the compositions becoming more crust-like toward the east (Fig. 11A). On an Sr-Nd isotope plot (Fig. 11B), the samples define a trend that closely mirrors that seen in data for the Saddle Mountains member of the Columbia River Basalts (Hooper and Hawkesworth, 1993), which has been attributed to crustal assimilation, although no specific crustal end member(s) has(have) been identified. Sample EI07-37, from the Silver Point Quartz Monzonite, lies off the array defined by the other samples, with lower ε_Ndi (−17.9; Fig. 11B). Whitehouse et al. (1992) reported similar Nd isotopic data from this pluton.

Oxygen isotope ratios of quartz gave δ¹⁸O = 8.4‰–10.0‰ for all but one sample in both Cretaceous–Paleocene and Eocene plutons (Table 4; Figs. 11D and 11E). These values are significantly higher than expected for rocks fractionated from mantle-derived melts: Mantle δ¹⁸O values are 5.5‰–6.0‰ (e.g., Eiler, 2001; Bindeman, 2008), and fractionation from basaltic to rhyolitic compositions should increase δ¹⁸O by <1‰ (Bindeman et al., 2004). Quartz-melt isotopic fractionation should be +1‰ or less for evolved rhyolitic compositions, and refractory quartz is generally considered to be resistant to secondary exchange (Bindeman et al., 2004). Thus, δ¹⁸O_quartz values >7‰ indicate crystallization of “high-δ¹⁸O” magmas.

Although Challis-age plutonic rocks have many geochemical and mineralogical traits of subduction-related granitoids, they also display subtle differences that point to generation under drier and hotter conditions. We summarize two such differences in this section, comparing data for Challis rocks with data from older (>52 Ma) plutonic rocks from northeastern Washington and from Cascade arc plutonic rocks.

B/Be ratios of subduction-related magmas (4–186) are higher than those of magmas from other tectonic settings. This is attributable to incorporation of a slab-derived hydrous fluid enriched in B, a fluid-mobile element (Morris et al., 1990; Leeman, 1996). B contents and B/Be also vary among arcs, with lower values in younger, hotter arcs, such as the Cascade arc, where slab devolatilization occurs at more shallow levels (Leeman et al., 2004). Eocene plutons and dikes in northeastern Washington have extremely low B contents (2.2–8.0 ppm) and B/Be (0.8–2.2), lower than values observed in lavas (3.4–14.0 ppm B, B/Be = 2.4–7.8) or plutonic rocks (6.0–32.9 ppm B, B/Be = 2.9–33.2) of the Cascade arc (Fig. 12). These low values, comparable to data for ocean-island basalts (Ryan et al., 1996), are indicative of melting under drier conditions than are typical of subduction zones.

Crystallization of amphibole in magmas of intermediate composition requires ~4 wt% or greater H₂O (Naney, 1983), and amphibole fractionation is thus typically restricted to hydrous magmatic systems. Geochemical consequences of amphibole fractionation include higher SiO₂ and depletion in middle rare earth elements (MREEs), which can be expressed as (Dy/Yb)_N (Davidson et al., 2013). Among data for northeastern Washington Cretaceous and Paleocene plutons and Cascade arc batholiths (Fig. 8F), a decrease in (Dy/Yb)_N is evident for samples with ≥70 wt% SiO₂. Conversely, no similar decrease is observed in data for Eocene samples, suggesting that those rocks formed from drier magmas that achieved amphibole saturation only in the late stages of crystallization. Late saturation is also consistent with textural data, which show that amphibole in the Eocene plutons is generally a late-formed phase, commonly with pyroxene cores. The lack of amphibole fractionation, which would drive residual melts to higher silica contents, may also explain why samples with ≥72 wt% SiO₂ are scarcer among Eocene than Cretaceous and Paleocene rocks (Fig. 6).

The wide range of Sr-Nd isotopic compositions that become more radiogenic to the east (Fig. 11A) indicates significant involvement of continental crustal material of varying age. Crustal contributions are also evident from the geographic distribution of depleted mantle model ages (T_DM) calculated from Nd data. These ages, taken to represent the time since the crustal source material was extracted from a depleted mantle, were determined by projecting Nd isotopic ratios of individual samples back in time to the point of intersection with an evolving depleted mantle reservoir having ¹⁴³Nd/¹⁴⁴Nd₀ = 0.51315 and ¹⁴⁷Sm/¹⁴⁴Nd = 0.2137. The results delineate three geographic regions having distinct T_DM that become younger toward the west (Fig. 13) and closely align with the crustal domains of Lund et al. (2015). The western region, which corresponds to the Insular and Quesnellia domains and has T_DM <950 Ma, is characterized by 16 samples: 10 dikes (923–455 Ma), five Cretaceous and Paleocene plutons (899–501 Ma), and one Eocene pluton (621 Ma). The central region, which corresponds to the Okanagan and Kootenay domains and has T_DM = 1300–1000 Ma, includes eight samples: five Cretaceous and Paleocene plutons (1279–1008 Ma), and three Eocene plutons (1282–1107 Ma). The eastern region, which corresponds to the Pend Oreille domain and has T_DM >1300 Ma, is defined by data from five Cretaceous and Paleocene plutons (2251–1432 Ma) and four Eocene plutons (2000–1380 Ma), as well as 23 T_DM ages for Cretaceous and Eocene plutons analyzed by Whitehouse et al. (1992). This region has the greatest spread in T_DM and likely encompasses more than one basement unit. The cluster of older T_DM ages (2200–1800 Ma) in the vicinity of the Silver Point Quartz Monzonite was attributed by Whitehouse et al. (1992) to an older than 2600 Ma deep crustal source that experienced early high-grade metamorphism (Rb and Pb depletion). The extent of this source is unknown, but there is no isotopic evidence for its involvement in any other Washington plutons included in this study. A few samples in the central and eastern regions have anomalously young T_DM ages (e.g., sample PL15-025 with T_DM = 745 Ma) that may reflect instances of mixing where there was involvement of a mantle end member. In none of the regions is there any systematic difference between T_DM of Eocene versus Cretaceous and Paleocene rocks or of dikes versus plutons, suggesting that all magmas within each region, regardless of age, shared the same deep crustal source(s).

Uniformly high quartz δ¹⁸O values (>8‰) over the spectrum from low to high SiO₂ bulk compositions support derivation from a high-δ¹⁸O source(s), and not fractionation or mixing with mantle compositions (Fig. 11E). Conversely, most Columbia River Flood Basalt isotopic values define a normal δ¹⁸O mantle array (Carlson, 1984; Brandon et al., 1993), with only a few higher δ¹⁸O (and ⁸⁷Sr/⁸⁶Sr) values that track crustal contamination. Brandon et al. (1993) noted that rocks from the Blue Mountain province and accreted terranes in southern British Columbia have similar elevated oxygen isotopic values, along with high ⁸⁷Sr/⁸⁶Sr, suggesting that crustal contamination is responsible for the high δ¹⁸O of rocks in our study as well.

Our new Sr isotopic data also allow us to refine the location of the “⁸⁷Sr/⁸⁶Sr = 0.706 line” (Armstrong et al., 1977) in northeastern Washington to ~150 km farther east than originally drawn. This revised location (Fig. 12) closely traces the boundary between the ensialic Pend Oreille domain to the east and accreted terranes to the west (Lund et al., 2015), but it does not correspond to an abrupt change in whole-rock chemistry.

Multiple lines of evidence indicate that, within each geographic region, Cretaceous and Paleocene and Eocene magmas were derived from the same lower-crustal source(s): They overlap in Sr-Nd-O isotopic values, show the same regional patterns in T_DM ages, and are broadly similar in major-element compositions. However, trace-element data suggest the mineralogy of the source regions changed over time, with the Cretaceous and Paleocene magmas originating from amphibole-bearing gabbroic sources, whereas the Eocene magmas originated from eclogitic sources. This change, which is reflected in MREE depletions [i.e., (Dy/Yb)_N < 1] in the Cretaceous and Paleocene rocks and in the lower average SiO₂ contents and lack of negative Eu anomalies in the Eocene rocks (Figs. 6, 8E, and 9), indicates generation of the latter at greater depth, below the stability limit of feldspar. There is no indication the Eocene rocks formed by a higher degree of partial melting: Incompatible trace-element contents, aside from those controlled by feldspar or garnet, are similar for Cretaceous and Paleocene and Eocene samples (e.g., La, Rb, U, Th; Table 3). The generally less felsic character of the Eocene plutons is thus dominantly a consequence of derivation from an eclogitic versus gabbroic source.

Values of (La/Yb)_N, which can be a proxy for crustal thickness (Farner and Lee, 2017; Shea et al., 2018), show a gradual increase over time among Cretaceous and Paleocene plutons, rising from <10 at 113 Ma to 48 at 51.5 Ma (Fig. 14A). We suggest this increase reflects development of the orogenic plateau, which is thought to have formed in response to Cretaceous terrane accretion (Whitney et al., 2004). In the central Andes, a similar history of crustal thickening and plateau uplift, attributed to ductile thickening of lower crust that had been heated by mantle-derived magmas, is also reflected most clearly in changing rare earth element (REE) patterns of felsic magmas (Kay and Mpodozis, 2001). In northeastern Washington, subsequent collapse of the plateau is recorded by an abrupt drop in (La/Yb)_N, which dropped to an average of 18 for plutons ca. 52 Ma or younger in a span of <1 m.y. (Fig. 14A). Sr/Y data, another proxy for crustal thickness (Chapman et al., 2015), tell a similar story of plateau development and collapse (Fig. 14B). However, this ratio is sensitive to feldspar fractionation, which probably explains why scatter is greater on this plot than on the (La/Yb)_N plot.

The Eocene dikes, which occur only in the western domain, likely originated from the same lower-crustal source(s) as the collocated Eocene plutons: Their average major- and trace-element contents are very similar, and they have similar T_DM ages. The broader range of chemical compositions defined by the dikes (Fig. 6) suggests they are the products of unblended smaller melting events that more fully reflect the compositional diversity of the source regions and melting processes. Conversely, plutons probably represent either larger-scale melting events, such that source heterogeneity was blended prior to magma ascent, or aggregations of smaller melt batches that accumulated and homogenized at shallower crustal levels. In either scenario, the absence of more crustal Sr-Nd isotopic compositions in the Eocene dikes (e.g., lowest ε_Nd = −1 for dikes vs. −15 for plutons; Fig. 11B) reflects their geographic restriction to the western part of the study area, where plutons have similar Sr-Nd isotopic values. Notably, despite their more mantle-like Sr-Nd isotopic values, these dikes have high δ¹⁸O values similar to those of high-⁸⁷Sr/⁸⁶Sr rocks to the east. These high values further support the contention that the dikes, although less radiogenic, also represent melts of preexisting lower crust with little or no evidence of mixing with mantle melts.

The presence of primary muscovite and magmatic epidote indicates some, if not most, Cretaceous and Paleocene plutons were emplaced at midcrustal levels, deeper than the Eocene intrusions. This difference is attributable to the change in stress regime: Ascent of Cretaceous and Paleocene magmas would have been inhibited by the compressional stress regime, whereas ascent of Eocene magmas would have been facilitated by the shift to an extensional setting. With collapse and spreading of overthickened crust, more rapid ascent, and thus less time for differentiation, may also explain the absence of negative Eu anomalies in most Eocene rocks. Negative Eu anomalies (Eu/Eu* < 0.8) are more common in Cretaceous and Paleocene rocks, but these could be a result of either residual feldspar left in the source and/or greater fractionation during passage through the crust.

Magmatism in northeastern Washington between 113 and 46 Ma occurred over an interval during which the dip of the subducting plate was changing. Widespread Cretaceous plutonism, with peaks at 113–110, 92–87, and 80–77 Ma, was part of the arc that extended from British Columbia to Idaho and was presumably associated with a moderate slab dip. This was followed (64–53 Ma) by a time of rapid, flat-slab subduction during which magmatism was scattered and sparse. Small-volume felsic magmas in modern-day flat-slab settings have been attributed to crustal melting (Ramos and Folguera, 2009), and a similar process likely explains Paleocene–early Eocene plutonism in our study area, perhaps facilitated by leaky transform faults (Cosca et al., 2021) or “edges” in the subducting plate (Yogodzinski et al., 2001) that allowed asthenospheric melts to impinge upon the crust.

Following this quiet interval, and coincident with orogenic plateau collapse, igneous activity increased between 52 and 46 Ma and migrated ~300 km southwestward across the study area (Fig. 5B). This corresponds to an average migration rate of 5 cm/yr, although the rate may have been slower during the final 2–3 m.y. based on the closer spacing of samples with ages between 49 and 46 Ma. Two possible causes for a migrating belt of magmatism are passage of a slab window and slab rollback. In this case, the oblique orientation of the Kula-Farallon or Farallon-Resurrection ridge relative to the margin, in combination with northward migration of the oceanic plates, should have produced a north- or northwest-migrating slab window (Madsen et al., 2006). Conversely, rollback should have proceeded roughly along the azimuth of convergence but toward the trench, resulting in southwest-directed magmatic progression. The observed pattern of ages among Eocene plutons in northeastern Washington clearly points to the model of slab rollback magmatism. We suggest that rollback, and consequent upwelling of hot mantle, probably facilitated plateau collapse by heating and weakening the lower crust.

Tomographic studies of the mantle beneath northwestern North America provide additional support for the slab breakoff and rollback model (Schmandt and Humphreys, 2011; Fuston and Wu, 2021). These studies delineate the geometry of the subducting Cascadia (Juan de Fuca) slab as well as the existence of detached fragments of older slabs and of areas where there is no slab. The Cascadia slab, which connects upward to the modern subduction zone, consists of two segments that meet at ~45°N latitude (Fig. 15). The longer, southern segment extends as far east as Georgia and represents >75 m.y. of Farallon plate subduction, whereas the shorter, northern segment extends only as far east as central Montana and represents ~52 m.y. of subduction (Fuston and Wu, 2021). Schmandt and Humphreys (2011) imaged a high-velocity “curtain” beneath northeastern Washington and Idaho that they interpreted as a detached segment of the Farallon slab; if this “Idaho slab” is added to the northern segment of the Cascadia slab, a total of ~60 m.y. of subduction is represented (S. Fuston, 2022, written commun.). The northern boundary of the Cascadia slab, located at the northern end of Vancouver Island, is abrupt and marks the beginning of a slab-free region that underlies much of western Canada (Fig. 15). This region originated as the Farallon-Resurrection slab window and has been enlarged by thermal erosion (Fuston and Wu, 2021).

The locations of these mantle features relative to the overlying North America plate would have been different at the time of the Challis episode. Undoing 52 m.y. of subduction, the leading edge of the northern Cascadia slab would have been under western Washington, and its northern boundary (the Farallon-Resurrection slab window) would have been near the U.S.-Canada border. This geometry agrees well with the distribution of 52–45 Ma igneous activity: The northern limit of Farallon slab rollback appears to have been near the border, and the belt of igneous activity attributed to slab breakoff (Kant et al., 2018) is in western Washington (Fig. 15).

Major differences between the Challis belt and its continuation to the north are: (1) its greater width (Fig. 1), (2) its higher proportion of plutonic versus volcanic rocks, particularly in northeastern Washington, and (3) its systematic younging-to-the-southwest age progression. Volcanism north of the border occurred over the same ca. 55–45 Ma time interval, but with no discernible age pattern (Dostal and Jutras, 2022). Breitsprecher et al. (2003) noted a southward increase in alkalinity, with Kamloops rocks being more arc-like and Challis rocks having within-plate affinities. In recognition of the large volume of magmas produced over a short time interval, Dostal and Jutras (2022) attributed Kamloops-Challis activity to slab breakoff and inferred that regional patterns in isotopic composition (e.g., the oldest model ages occur in areas of Precambrian basement) indicate differences in the age of mantle metasomatism by slab-derived fluids and not crustal contamination. We suggest that while slab rupture occurred both north and south of the U.S.-Canada border, rollback occurred only in the south, causing the wider breadth and age progression of magmatism as well as greater uplift. Based on the correlation between Sr-Nd isotopic compositions and SiO₂ contents of Challis-age rocks in Washington and their uniformly high δ¹⁸O values, it appears that melting of isotopically heterogeneous lower crust, rather than mantle metasomatism, was the dominant control on Sr-Nd and O isotopic compositions of the rocks in this study.

Cretaceous to Eocene plutonic rocks in northeastern Washington document an ~60 m.y. history of crustal thickening and subsequent collapse that reflects two contrasting responses to terrane accretion. Arrival of the Insular terrane at ca. 100 Ma resulted in thrust faulting and ductile thickening, producing an orogenic plateau with crustal thicknesses exceeding 60 km. This was followed ~50 m.y. later by arrival of the Siletzia terrane, a young, buoyant large igneous province that was partially thrust under the edge of the continent but that was ultimately not subductable (Trehu et al., 1994; Wells et al., 2014). That collision caused the Farallon slab to rupture, perhaps under present-day Montana (Duke et al., 2014), and then roll back toward the margin. Rollback triggered upwelling of asthenospheric mantle, producing the southwest-migrating front of crustal melting and plutonism that swept across northeastern Washington between 52 and 46 Ma. This transient thermal event would also have weakened the crust, facilitating collapse of the orogenic plateau, and was likely responsible for a topographic bulge that diverted drainage patterns as it migrated across western North America between 53 and 47 Ma (Smith et al., 2014). Rollback beneath Washington may also represent the onset of the broader rollback that propagated across the Great Basin during the Eocene and Miocene (Dickinson, 2006).

Differences in the emplacement depths of Cretaceous plutons, which were generally midcrustal, versus those of Eocene plutons, which were epizonal, also provide evidence of uplift and exhumation in this time interval and suggest that rapid and voluminous near-trench sedimentation, which has been linked to erosion in the Idaho batholith area, may have been sourced from rocks in Washington as well. The prominent peak in “Challis-age” detrital zircons from several basins in Oregon and California (Dumitru et al., 2013) further supports this idea.

Rollback culminated in a second rupture of the Farallon slab that occurred ca. 51 Ma under what is today central Washington. This event is recorded in a narrow NW-trending belt of 51–48 Ma bimodal volcanic units and S-type granites that approximately parallels the farthest eastward extent of Siletzia rocks in the subsurface (Kant et al., 2018). The high-velocity “curtain” imaged by seismic tomography beneath Idaho and interpreted as a detached remnant of the Farallon slab (Schmandt and Humphreys, 2011) is probably the result of this breakoff event.

The northern limit of the region affected by slab rollback coincides roughly with the U.S.-Canada border, where rocks north of the boundary show no indication of a spatial age progression. At 50 Ma, this was the latitude of the Farallon-Resurrection slab window (Fig. 15), which would also have been the northern edge of the Farallon slab as it rolled back. The presence of alkaline lavas in the Princeton Group in southern British Columbia, which have been attributed to melting above a slab tear (Ickert et al., 2009), further supports this interpretation. More broadly, the fact that rollback in this case involved only a limited portion of the much larger Farallon slab illustrates how slab windows, fracture zones, and other zones of weakness can “divide” a subducted plate into domains that behave differently and that are expressed at the surface by different patterns of magmatism and tectonism.

Much of this work is the product of B.S. thesis research conducted at the University of Puget Sound (UPS) and of term projects by the UPS Advanced Igneous Petrology classes of 2005 and 2014. Support was provided by National Science Foundation Research in Undergraduate Institutions (RUI) grant EAR-1119252 (to Tepper) and by summer research awards to individual students from the Agricola and McCormick Funds and the University Enrichment Committee. Ilya Bindeman graciously provided access to the University of Oregon Stable Isotope laboratory. Conversations with R.B. Miller, M.P. Eddy, and K.P. Clark helped to develop the story, and C. Serdar Tepper provided field support. We thank Ryan Frazer and Rich Gaschnig for detailed and thoughtful reviews that improved the manuscript, and Robinson Cecil and Mihai Ducea for editorial handling. This article has been peer reviewed and approved for publication consistent with U.S. Geological Survey Fundamental Science Practices (https://pubs.usgs.gov/circ/1367/). Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.