The southern Coast Mountain batholith was episodically active from Jurassic to Eocene time and experienced four distinct high magmatic flux events during that period. Similar episodicity has been recognized in arcs worldwide, yet the mechanism(s) driving such punctuated magmatic behavior are debated. This study uses zircon Hf and O isotopes, with whole-rock and mineral geochemistry, to track spatiotemporal changes in southern Coast Mountains batholith melt sources and to evaluate models of flare-up behavior and crust formation in Cordilleran arc systems. Zircon Hf isotope analysis yielded consistently primitive values, with all zircon grains recording initial εHf between +6 and +16. The majority (97%) of zircons analyzed yielded δ18O values between 4.2‰ and 6.5‰, and only five grains recorded values of up to 8.3‰. These isotopic results are interpreted to reflect magmatism dominated by mantle melting during all time periods and across all areas of the southern batholith, which argues against the periodic input of more melt-fertile crustal materials as the driver of episodic arc magmatism. They also indicate that limited crustal recycling is needed to produce the large volumes of continental crust generated in the batholith. Although the isotopic character of intrusions is relatively invariant through time, magmas emplaced during flare-ups record higher Sr/Y and La/Yb(N) and lower zircon Ti and Yb concentrations, which is consistent with melting in thickened crust with garnet present as a fractionating phase. Flare-ups are also temporally associated with periods when the southern Coast Mountains batholith both widens and advances inboard. We suggest that the landward shift of the arc into more fertile lithospheric mantle domains triggers voluminous magmatism and is accompanied by magmatic and/or tectonic thickening. Overall, these results demonstrate that the magmatic growth of Cordilleran arcs can be spatially and temporally complex without requiring variability in the contributions of crust and/or mantle to the batholith.

Continental arcs are important sites of continental crust generation and modification, and the large composite batholiths produced in arcs record the combined processes of partial mantle melting and crustal recycling. Most continental or Cordilleran-type arcs are constructed over long periods of time (>100 m.y.), and although subduction may be ongoing throughout those intervals, magmatism appears to be highly episodic (e.g., Armstrong, 1988; Ducea, 2001; Paterson and Ducea, 2015). This episodic behavior is characterized by high-flux magmatic events, or “flare-ups,” that quasi-periodically punctuate the baseline rates of magmatism. The processes that control the oscillation of arcs between flare-up and baseline modes remain debated but are generally thought to be either internal to the arc (i.e., tectono-magmatic processes operating in the upper plate such as periodic crustal thickening or delamination) or external to the arc (e.g., periodic change in the obliquity or rate of plate convergence). Flare-ups are capable of producing as much as 90% of the magmatic additions to arcs (Ducea and Barton, 2007), and therefore identifying the petro-tectonic mechanisms that control them is critical to our understanding of the growth of continental crust both in individual arc systems as well as globally.

Isotopic investigation of plutons in many Cordilleran arcs, particularly those emplaced along the margins of North and South America, has revealed significant incorporation (up to 50%; Ducea et al., 2015) of preexisting upper plate lithosphere in new arc melts. This process is evidenced by higher 87Sr/86Sr and lower εNd and εHf in batholithic rocks, which points to the incorporation of older, upper plate materials in magma sources; however, these metrics cannot be used to distinguish between contributions from crustal and mantle lithospheric components. Recycling of crust is reflected in the δ18O of igneous materials elevated above mantle values, which can only occur if surficial rocks that have isotopically exchanged with the hydrosphere are partially melted (Valley et al., 2005). Examples of these isotopic systems being interpreted to indicate upper plate recycling have been documented in the Sierra Nevada (Lackey et al., 2008, 2012; Cecil et al., 2012), the Peninsular Ranges batholith (Kistler et al., 2014), the Idaho batholith (King and Valley, 2001; Gashnig et al., 2011), and the Coast Mountains batholith (Samson et al., 1991; Wetmore and Ducea, 2011; Cecil et al., 2011; Girardi et al., 2012). Furthermore, upper plate additions to batholithic melts have been temporally linked to high-flux events (e.g., Ducea and Barton, 2007), which has led to models wherein arc flare-ups are fueled by the periodic incorporation of more melt-fertile, crustal materials (Ducea and Barton, 2007; DeCelles et al., 2009). In contrast, several recent studies have suggested that flare-ups are dominated by juvenile melts and driven by episodic events in the mantle (Decker et al., 2017; Martínez Ardila et al., 2019; Attia et al., 2020). In these cases, moderately depleted εNd and εHf signatures are attributed to contributions from enriched mantle sources.

The primary goals of this study are to test the possibility that magmatic flare-ups in the Coast Mountains batholith are driven by the episodic addition of supracrustal materials to melt sources and to develop a petro-tectonic model to explain the complex spatio-temporal magmatic patterns observed in the batholith. The southern part of the batholith is an interesting and ideal site for this test because it provides continuous, arc-perpendicular exposures of intrusive rocks that have variable compositions and that span the lifetime of the active arc. Geochronologic analysis of these intrusions revealed a complex episodic magmatic history and one that is distinct from the batholith to the north, yet the processes controlling arc tempos remain uncertain (Cecil et al., 2018). It appears that magmatism is at times—but not always—synchronous with periods of crustal deformation (Rusmore et al., 2013, 2019) and metamorphism (Bollen et al., 2017; Dafov et al., 2020). By investigating the geochemical and isotopic compositions of Jurassic–Eocene intrusions of the southern Coast Mountains batholith, we can evaluate: (1) the extent to which crust and mantle sources are involved in batholith generation through time and space; (2) the depth and thermal conditions of melt generation; and (3) links between magmatism, deformation, and metamorphism.

Our results indicate that magmatism in the southern Coast Mountains batholith is dominated by mantle melting, which is recorded in primitive zircon εHf and uniformly mantle-like zircon δ18O values. Input of older and/or near-surface country rocks is limited, which indicates that melt-fertile supracrustal materials cannot be the primary drivers of magmatic flare-ups. Crustal recycling is therefore not likely to be a volumetrically significant process within the southern Coast Mountains batholith, which instead reflects the generation of new continental crust (Hollister and Andronicos, 2006). The processes controlling mantle-driven flare-ups remain uncertain, but spatial and temporal trends in the geochronologic and geochemical data suggest that periodic inboard advance of the Coast Mountains arc may lead to episodic tapping of more fertile sub-continental mantle lithosphere, which would produce the mantle-like flare-ups observed.

The Coast Mountains batholith is one of the largest exposed Cordilleran arcs; it extends >1500 km along the western margin of southeast Alaska and British Columbia. Subduction of paleo-Pacific plates beneath western North America drove the nearly continuous growth of the batholith for ca. 150 m.y. in Jurassic to early Eocene time, whereupon the margin began transitioning to a dextral transform system (e.g., Stock and Molnar, 1988; Madsen et al., 2006). The southern part of the Coast Mountains batholith investigated as part of this study is exposed between ∼50.5 °N and 52 °N (Fig. 1). This is one of the widest parts of the batholith, with very little host rock exposed, making the affinity of the basement into which the southern Coast Mountains batholith is emplaced uncertain. Vancouver Island, which comprises rocks of the Wrangellia terrane, is located to the west of the study area, whereas Triassic–Cretaceous rocks of the Intermontane superterrane outcrop to the east; the location of the boundary between the two is unknown. Paleomagnetic results, with patterns observed from geochronologic data sets, suggest that the Intermontane superterrane may extend westward to Vancouver Island and underlie much of the batholith at the study latitudes (Rusmore et al., 2013). This is also consistent with recent detrital zircon results from metasedimentary pendants preserved in Bute and Knight Inlets (Dafov et al., 2020).

The southern Coast Mountains batholith was emplaced between ca. 170–45 Ma, and although magmatism was continuous during that time, the rates of magmatic addition to the batholith varied considerably. The southern batholith experienced four magmatic flare-up events at 161–148 Ma, 114–102 Ma, 85–70 Ma, and 61–48 Ma with intervening apparent magmatic lulls at 148–114 Ma, 102–85 Ma, and 70–61 Ma (Cecil et al., 2018). The Jurassic event coincides with a Cordilleran arc flare-up documented along most of the North American margin (Gehrels et al., 2009; Kirsch et al., 2016; Beranek et al., 2017). The timing of Cretaceous and Paleocene–Eocene high flare-ups differs from those documented in the central and northern parts of the batholith (Gehrels et al., 2009), which suggests a relatively small-scale (i.e., not plate-scale) control on high-volume magmatic events (de Silva et al., 2015; Kirsch et al., 2016; Cecil et al., 2018). Jurassic to Early Cretaceous intrusions are restricted to the western part of the southern Coast Mountains batholith and magmatism generally becomes younger to the east, such that post-80 Ma plutons are only found at distances greater than 100 km inboard of the coastline. Eastward migration of the post-mid-Cretaceous Coast Mountains batholith has been recognized along the length of the arc, and although arc migration rates have been estimated at between 2.1 km/Ma and 2.8 km/Ma for all arc segments, the spatial evolution of arc magmatism may be more complex than was previously recognized (Gehrels et al., 2009; Cecil et al., 2018).

Whole-Rock Geochemistry

Many of the intrusions analyzed in this study are the same as those dated and discussed in Cecil et al. (2018). Sample information and location data are provided in Supplemental Data Table S11. During the zircon separation process, ∼10 g aliquots of crushed rock were set aside for whole-rock geochemical analysis. Additional hand samples of large or complex plutons were also collected and crushed to characterize the geochemistry of potentially heterogeneous intrusions. Unweathered rock chips selected from crushed whole-rock samples were powdered in an alumina ceramic ring and puck mill and analyzed for major and trace elements either via X-ray fluorescence (XRF) at Pomona College or by inductively coupled plasma-mass spectrometry (ICP-MS) at Activation Laboratories in Ontario, Canada. Samples analyzed by XRF at Pomona College were also digested in acid and then analyzed for rare earth elements by ICP-MS at California State University, Northridge, California, USA (CSUN). All whole-rock geochemical data and modes of analysis are given in Table S2 (see footnote 1).

Zircon Hf Isotope Analysis

Hafnium isotope analysis was performed on dated zircon grains via laser-ablation (LA)-ICP-MS at the Arizona LaserChron Center using a Nu Instruments multi-collector mass spectrometer coupled with a Photon Machines Analyte G2 excimer laser system. The analytical procedures, and mass bias and interference corrections used, follow those described in Cecil et al. (2011) and Gehrels and Pecha (2014). Hafnium isotope data were acquired for a total of 775 zircons from 54 plutonic samples. These analyses are representative subsets (∼10–15 grains per sample) of the same zircons dated by Cecil et al. (2018). Most zircon analyses reported in Cecil et al. (2018) exhibited simple U-Pb systematics and did not display evidence of inheritance or overgrowths, as reflected in discrepant ages or unusual textures in cathodoluminescence (CL) images (e.g., thin, bright rims or cores with morphologies different that those of the crystal's external shape) (Fig. 2). Nevertheless, only zircons that were concordant, had CL textures interpreted to be magmatic (oscillatory or sector zonation), and contained only a single age domain, were analyzed. Wherever possible, Hf analysis was performed on top of existing U-Pb ablation pits to ensure the greatest likelihood of analyzing the same age domain sampled for U-Pb geochronology. Hafnium isotopic data for individual grain analyses are given in Table S3 (see footnote 1).

Zircon O Isotope Analysis

Following zircon U-Pb and Hf isotope analysis, a subset of 22 samples was chosen for O isotope and trace element analysis of zircon. Samples were selected from all age groups and from a wide variety of locations and rock compositions to broadly characterize the Coast Mountains batholith in our study area. New mounts were prepared from existing zircon separates and imaged by SEM-CL. Therefore, age and Hf isotope information cannot be directly connected to a given O isotope zircon analysis spot. However, given the intra-sample homogeneity of U-Pb ages and Hf isotope ratios (see Fig. 2), we are reasonably confident that the measured O isotope values correspond to previously measured average U-Pb-Hf characteristics of a given intrusion. Oxygen isotope analyses were performed using a CAMECA IMS-1280 ion microprobe at the University of Wisconsin–Madison, Wisconsin, USA, following procedures described in detail elsewhere (Kita et al., 2009; Valley and Kita, 2009). In summary: a 133Cs+ primary beam (∼2 nA) was focused to ∼10 micron spots on the polished zircon surface; three Faraday Cup detectors simultaneously analyzed 16O, 16OH, and 18O; and eight analyses of KIM-5 zircon standard (Valley, 2003) bracket every 10–20 sample analyses. Values of δ18O are reported in standard permil notation relative to VSMOW. Ratios of OH/O, corrected for background, were monitored as an indication of inclusions or radiation damage (Wang et al., 2014), which resulted in the rejection of two out of 175 analyses. The average spot-to-spot reproducibility (external precision) for δ18O on KIM-5 for brackets of eight analyses is 0.19‰ (2 standard deviations). Oxygen isotope data for individual grain analyses are given in Table S4 (see footnote 1).

Zircon Trace Element Analysis

Zircon trace element analysis was performed on the same sample subset as that for which O isotopic data were obtained. Approximately 15 grains for each of the 22 samples were analyzed for a total of 359 analyses. Trace element data were collected via LA-ICP-MS at the University of Arizona LaserChron Center in Tucson, Arizona, USA, using a Photon Machines G2 excimer laser-ablation (LA) system coupled to a Thermo Fisher Scientific Element 2 high-resolution single collector ICP-MS. A total of 24 different masses were analyzed, including the rare earth elements (REE), as well as P, Ti, Nb, Hf, Ta, Pb, U, and Th. The natural zircon standard 91500 (Wiedenbeck et al., 2004) was used as a primary standard for calculating concentrations for all elements except for Ti, which was standardized using NIST 612 glass. Analytical methods and data reduction protocols used were the same as those described in Chapman et al. (2016). Trace element data for individual zircon analyses are given in Table S5 (see footnote 1).

Major and Trace Element Geochemistry

Whole-rock geochemical data were generated for 96 samples collected from intrusive rocks of the southern Coast Mountains batholith. Pluton compositions range from gabbro to granite, though the majority of the batholith is composed of tonalite and granodiorite with SiO2 contents >65% (Fig. 3A). Plutons emplaced during the most recent magmatic event in the arc, at 64–45 Ma, have higher average SiO2 (68.3% SiO2 as compared with an overall average of 64.9% SiO2) and are restricted to the easternmost part of the southern batholith. Southern Coast Mountains batholith intrusions are metaluminous to mildly peraluminous (alumina saturation index (ASI) = A/CNK = 0.8–1.15) and show a smooth trend of increasing ASI with increasing SiO2 (Fig. 3B). Overall, major element data show that this part of the batholith is magnesian, calcic to calc-alkalic, and generally similar to rocks of the central Coast Mountains batholith (Girardi et al., 2012) and other Cordilleran arc granitoids (Frost et al., 2001) (Fig. 3).

Although major element data reveal little spatio-temporal variability in pluton composition, there are changes in both La/Yb and Sr/Y through time. In general, plutons emplaced during high-flux magmatic events have elevated La/Yb and Sr/Y, whereas those emplaced during periods of relative arc inactivity do not (Fig. 4). For example, 46% of intrusions with ages corresponding to a documented flare-up period have Sr/Y>50 as compared with only 12% of intrusions emplaced during magmatic lulls.

Zircon Hf and O Isotope Compositions

Initial εHf values in Coast Mountains batholith intrusions range from ∼+6–+16, cluster around a mean of +11.4 ± 1.9, and are distinct from the distribution of εHf(t) values documented in the central and northern parts of the batholith (Cecil et al., 2011) (Fig. 5A). Samples from all age groups have zircons that overlap with depleted mantle values. There are no obvious arc perpendicular trends observed in the Hf isotope data, though plutons with the lowest average values (+7.7 to +8.9) are found in axial arc positions and roughly coincide with a NW-trending package of gneisses. Temporal trends in the Hf data are equally unclear. A running mean calculated for post-145 Ma εHf(t) values reveals limited variability, and mean values range from ∼+11 to +13 (Fig. 6). Contouring of individual Hf data points, however, does show values scattering to lower εHf(t) values at ca. 155–145 Ma, a time during which a high magmatic flux event is documented for the Coast Mountains batholith and most other segments of the North American Cordillera (Gehrels et al., 2009; Kirsch et al., 2016; Beranek et al., 2017; Cecil et al., 2018).

Zircon δ18O values from intrusive rocks of the southern Coast Mountains batholith range from +4.2‰ to +8.3‰, though the majority of grains analyzed cluster between +5‰ and +6‰, with an average of +5.4 ± 0.6‰ (1s) (Fig. 5B). Average δ18O values for all plutons analyzed overlap with the accepted range of mantle values (5.3 ± 0.6‰ (2s); Valley et al., 2005). Because of the relative invariability of the zircon δ18O data, they do not reveal any obvious spatial or temporal trends.

Zircon Trace Element Geochemistry

Chondrite-normalized rare earth element (REE) trends in southern Coast Mountains batholith zircons show typical patterns of increasing concentration from light to heavy REEs, with marked positive Ce and negative Eu anomalies (Fig. 7A). There are no significant changes in these trends (e.g., steepness of heavy REE pattern, size of Ce or Eu anomaly, etc.) through time or space. The majority (95%) of zircons from southern Coast Mountains batholith plutons yield U/Th ratios <10; most of the 5% of grains with U/Th ratios >10 have corresponding U-Pb ages of ca. 140–150 Ma (not shown), which is also a pattern observed in metamorphic rims of detrital zircons in metasedimentary country rocks (Dafov et al., 2020). Hafnium concentrations range from 7000–14,000 ppm, are positively correlated with U/Yb, and are generally very similar in their Hf and U/Yb variation to globally compiled data for continental arcs (Fig. 7B) (Grimes et al., 2015). A bivariate plot of Ti vs. Yb concentration again shows that our Coast Mountains batholith zircon data are generally very similar to those of arcs worldwide, but it also suggests that intrusions crystallizing during interpreted flare-up events trend toward lower Ti and Yb (Fig. 7C).

Isotopic Evidence for Mantle-Dominated Batholith Sourcing

Intrusive rocks comprising the southern Coast Mountains batholith have zircon εHf(t) and δ18O values that are relatively invariant and that overlap with documented mantle values for both isotopic systems (Fig. 5). Zircon δ18O sample averages range from 4.5‰ to 6.4‰ and, in all cases, overlap within analytical error the range of accepted values for the high-temperature mantle (5.3 ± 0.6‰; Valley et al., 2005). Zircon εHf(t) values are positive and high (between ∼+6 and +16), which is compatible with melting of a Lu-rich mantle source. In all geographic areas and in all stages of batholith development—apart from the Eocene, for which we have relatively little data—individual εHf(t) values overlap with the evolution of depleted mantle (Vervoort and Blichert-Toft, 1999; Dhuime et al., 2011) (Fig. 6). Although these data do not uniquely signify a mantle source, they do require melting of materials that have not exchanged with δ18O-enriched surface fluids and that are young and therefore characterized by juvenile εHf values. Therefore, the isotopic data presented here are consistent with partial melting of variably enriched mantle, secondary melting of mantle-derived mafic underplates (Collins et al., 2020), remelting of juvenile arc crust, or a mixture thereof. The bulk isotope modeling presented in the subsequent section, with the isotopic homogeneity and the lack of zircon inheritance observed, indicates that remelting of lower crustal arc terranes likely played an insignificant role.

The remarkable inter- and intra-sample uniformity observed in both the Hf and O data sets is one of the most distinctive features of the southern Coast Mountains batholith and suggests batholith generation from a single, relatively homogeneous source. Zircon δ18O values cluster tightly around a mean value of 5.4‰ (Fig. 5B), and although they are slightly more heterogeneous, zircon εHf(t) sample averages comprise a range of only seven epsilon units (Fig. 6). Even intrusions with SiO2 contents >70% yield homogeneous, mantle-like values (for example, 14IY11 has SiO2 = 71.2% and εHf(t) = 13.3), which suggests that fractionation of a mantle partial melt, in the absence of crustal assimilation, is sufficient to explain most of the isotope data. These results are intriguing, given the large geographic area over which the southern batholith was sampled (>15,000 km2; Fig. 1) and the large variability in pluton composition (Fig. 3).

Despite the overall isotopically primitive signature recorded in the zircon data, some deviation from purely mantle-derived melts is observed. There is a slight pull-down in εHf(t) values at ca. 150 Ma and perhaps another less defined pull-down at ca. 110 Ma (Fig. 6). Only in these times do average pluton εHf(t) values move out of the “juvenile” and into the “moderately juvenile” fields defined by Bahlburg et al. (2011) as reflecting partial melts of sources that were extracted from the mantle 0–300 m.y. and 300–650 m.y., respectively. Therefore, although there is little evidence for significant contributions of Precambrian crust to southern Coast Mountains batholith melts, contamination by small amounts of relatively young (Phanerozoic) crust is permissible. This contamination may be evident in Jurassic samples that have lower average εHf(t) and zircons with slightly elevated δ18O (+5.5‰ to + 6.8‰); one sample (14IY25, from Bute Inlet) contains a single grain with measured δ18O of +8.3‰. The amount and nature of the crustal-end-member contaminant permitted by our isotope data is explored in the next sections using a compilation of available Hf isotope data from Coast Mountains batholith basement terranes and Hf-O bulk mixture models.

Incorporation of Probable Host Rocks into Batholithic Melts?

Modeling of the potential crustal component is hindered by uncertainty as to the rocks underlying much of the southern Coast Mountains batholith. Sparse pendants of weakly to strongly metamorphosed volcanic and sedimentary rocks traditionally have been considered part of Wrangellia, a terrane that is well-exposed on Vancouver Island and the westernmost mainland (Nelson, 1979; Monger et al., 1982; Wheeler and McFeely, 1991). In contrast, magmatic patterns and paleomagnetic results suggest that the Coast Mountains batholith at the study latitudes is emplaced into Intermontane (likely Stikine) terrane (Rusmore et al., 2013). Detrital zircon data from recent work on metasedimentary pendants in the southern batholith generally support that interpretation but are also compatible with connection to Wrangellia terrane and possibly Alexander terrane as exposed in the Banks Island region to the northwest (Dafov et al., 2020). We approach estimating the contributions of partial melts of these terranes into Coast Mountains batholith magmas through bulk mixture modeling of zircon Hf-O with the understanding that this modeling is hampered by the fact that the isotopic compositions of the crustal end members are poorly known.

Stikine terrane is composed primarily of Late Paleozoic through Jurassic island arc-related volcanics with juvenile isotopic signatures (Dostal et al., 1999, 2009; Barresi et al., 2015). These arc rocks overlie and/or interfinger with carbonates and siliciclastic rocks (Nelson et al., 2006) thought to have been deposited on pericratonic basement (Jackson et al., 1991; Gehrels and Kapp, 1998). Wrangellia terrane, as exposed on Vancouver Island, is dominated by Paleozoic and Mesozoic igneous rocks, including Devonian–Mississippian bimodal volcanics of the Sicker Group and similarly aged intrusives of the Saltspring Suite (Ruks, 2015), extensive Triassic basalts of the Karmutsen Formation (Greene et al., 2009), and Jurassic bimodal volcanics and intrusives of the Bonanza arc (DeBari et al., 1999). Also present are chert, argillite, limestone, and minor clastics of the lower Paleozoic Sicker Group and the upper Paleozoic Buttle Lake Group (Brandon et al., 1986; Massey and Friday, 1989; Yorath et al., 1999). The southern part of the Alexander terrane (Banks Island region) consists of Ordovician to Permian quartzite, marble, and metapelite assemblages (Tochilin et al., 2014).

Bulk mixture modeling is used to assess the extent to which assimilation of these probable crustal contaminants (Wrangellia, Alexander, Stikine) is permitted based on the zircon Hf-O data. Results of that modeling are shown in Figure 8. Based on the dominantly igneous nature of the Wrangellia and Stikine terranes, both of which contain young, hydrothermally altered oceanic crust and volcaniclastic sedimentary rocks, we group them together and assign them a minimum zircon δ18O of ∼+10 ‰, similar to zircon δ18O values assumed for accreted ocean arc terranes in the Sierra Nevada foothills (Lackey et al., 2012). On the basis of available zircon Hf data for Wrangellia (Alberts et al., 2021), whole-rock Nd from local metasedimentary pendants (Bollen et al., 2017) and whole-rock Nd for rocks of Stikine (Samson et al., 1989; Jackson et al., 1991), εHf(100 Ma) of these terranes is estimated to have been +5 to +13 when the southern Coast Mountains batholith was generated. We therefore assign this crustal endmember an average εHf of +8 and an Hf concentration averaged from values reported in Lassiter et al. (1995) (2.9 ppm). The Hf isotopic character of Alexander terrane is heterogeneous; most igneous and metaigneous basement is relatively primitive (εHf(t) = +5 to +15), whereas metasedimentary rocks of the Banks Island assemblage, which is thought to be correlative with the northern Alexander terrane, yield zircons with εHf(t) ranging from +12 to −30 (Tochilin et al., 2014). Given that the quartz-rich metasedimentary units that yield evolved values are volumetrically less significant, we assign an εHf(t) of 0 for the Alexander terrane. Hf concentration of Alexander is estimated to be intermediate between mafic rocks of Alexander reported by Israel et al. (2014) (1.8 ppm) and the average of siliciclastic sedimentary rocks reported by McLennan (2001) (5.5 ppm). The zircon δ18O of Alexander is taken to be +11‰, slightly higher than that of Wrangellia-Stikine, because although Alexander basement is not significantly different, there is greater potential involvement of supracrustal rocks. The primitive end-member is assigned an εHf(t) of +13.6 and an Hf concentration of 1.3, which are the values measured from the most mafic rock analyzed in our batholith data set (sample 15KN38B; SiO2 = 51%). We assume a mantle-like δ18O (Zrn) for this end-member (+5.3 ‰; Valley et al., 2005).

Results of Hf-O models show that the majority of southern Coast Mountains batholith samples form a tight cluster in Hf-O space near depleted mantle and do not reveal mixing trends between mantle and any of the likely crustal sources investigated (Fig. 8). There are four samples, however, which have grains recording non-mantle-like δ18O values, which indicates that these samples involve melting of material that had a near-surface history. Two of these samples are Jurassic plutons discussed previously, which also yielded some of the most evolved zircon εHf(t) values recorded in this part of the batholith (average εHf(t) of +8.4 and +8.6). A third, however, yielded zircons with narrowly distributed δ18O slightly above mantle values (6.2–6.4‰), and this pluton has a crystallization age of 138 Ma, a time during which there is a pronounced magmatic lull recorded along much of the North American Cordilleran margin (Gehrels et al., 2009; Paterson and Ducea, 2015; Kirsch et al., 2016; Beranek et al., 2017; Cecil et al., 2018). The fourth anomalous sample—a 108 Ma granodiorite collected from our eastern-most field location—yielded uniformly low δ18O values (4.2–4.7 ‰) just below the accepted mantle range. We interpret these results as indicating variable contamination of mantle-derived melts by heterogeneous supracrustal sources. That these contaminated samples have few commonalities (they range in age from 96 Ma to 153 Ma, are from plutons emplaced during flare-ups and lulls, and are spatially and compositionally diverse) suggests that there is no systematic process controlling the limited isotopic variability in the southern Coast Mountains batholith.

The lack of mixing trajectories in our data makes it difficult to identify any of the crustal end-members as possible melt components even in the slightly contaminated samples (Fig. 8). Overall, the relatively low δ18O and juvenile εHf(t) nature of southern Coast Mountains batholith samples precludes the involvement of >∼10% upper crustal materials into any components of the southern batholith. The isotope data, however, are compatible with sourcing from mafic underplates and/or deep crustal materials, which would not have exchanged with δ18O-enriched surface fluids, but in the case of young arc-related terranes would have juvenile εHf characteristics. Indeed, it is hard to explain the generation of a large and compositionally diverse batholith from direct melting of mantle sources. Partial melts of depleted mantle ponding in the lower crust and incorporating deep, previously generated arc rocks and/or the lower crust of Stikine and/or Wrangellia terrane could produce the isotopic and geochemical signatures observed in the southern Coast Mountains batholith. As recently described in Collins et al. (2020), cooling and fractionation of basaltic magmas in the lithospheric mantle can also cause water to exsolve, leading to flux-melting of preexisting mafic underplates. This process allows for the generation of large volumes of silicic, isotopically primitive melt in arc systems, much like what is preserved in the southern Coast Mountains batholith.

Drivers of Flare-Ups in the Southern Coast Mountains Batholith

Arc flare-ups likely occur as the result of increased melt fertility driven by (1) the introduction of volatiles or different rock compositions into the sub-arc or (2) changes in temperature or pressure conditions in the mantle wedge or upper plate lithosphere. One model for the former that has received a lot of attention calls on shortening in the rear-arc leading to periodic underthrusting of older, more evolved upper plate materials (Ducea, 2001; Ducea and Barton, 2007; DeCelles et al., 2009; DeCelles and Graham, 2015). In this orogenic cyclicity model, the introduction of hydrous retro-arc crust into the sub-arc melt generation zone ignites flare-ups. Because crust that is inboard of the arc tends to be older, flare-ups in other arcs, such as the Sierra Nevada, are characterized by isotopic “pull-downs,” so called because igneous rocks generated during these times record lower εNd values (Ducea and Barton, 2007; DeCelles et al., 2009). Zircon εHf values should record similar pull-downs during flare-ups because of the strong positive correlation between Nd and Hf isotopes in the terrestrial array (Vervoort and Blichert-Toft, 1999), though some recent studies in the Sierra Nevada have not supported this relationship (Attia et al., 2020; Klein et al., 2021). Periodic flare-ups and associated isotopic pull-downs are generally seen as an internal, autocyclic process. In contrast, an episodic change in mantle melt productivity could also lead to episodicity in arc magmatism, but potentially without the concomitant isotopic pull-down. Indeed, high-flux magmatic episodes without excursions to more crust-like isotopic signatures have recently been recognized in the Sierra Nevada and other continental arcs (Decker et al., 2017; Martínez Ardila et al., 2019; Attia et al., 2020; Klein et al., 2021).

In the southern Coast Mountains batholith, isotopic pull-downs in zircon Hf data are ambiguous at best and, with the O isotope data, preclude the significant involvement of supracrustal materials or old (Precambrian), Hf isotope-depleted crust. These data argue against flare-ups occurring in response to an internal underthrusting model wherein the periodic introduction of older, melt-fertile crust from the rear-arc drives high-flux events. Interestingly, however, whole-rock and zircon trace element data indicate that flare-ups may be temporally associated with periods of crustal thickening. Evidence for crustal thickening is seen in elevated whole-rock Sr/Y and La/Yb ratios (Fig. 4) and in lower Yb and Ti in zircon (Fig. 7C) during flare-ups. Although high Sr/Y and La/Yb ratios are not unique, surveys of those ratios in modern arcs show that Sr/Y and La/Yb are positively correlated with crustal thickness (Chiaradia, 2015; Profeta et al., 2015). Increases in these chemical indices have been linked to crustal thickening in several arc systems and are being increasingly used as proxies for paleo-thickness (e.g., Kay and Mpodozis, 2001; Schwartz et al., 2011; Chapman et al., 2015). Lower Yb and Ti in zircon is attributed to garnet and/or amphibole fractionation during zircon crystallization, which may indicate increased crustal thickness given the sensitivity of garnet stability to pressure-temperature conditions (Grimes et al., 2015).

Taken together, these geochemical trends suggest that the arc crust thickens and the depth of differentiation of batholithic melts increases during flare-up events. Periodic crustal thickening in the southern Coast Mountains batholith could result from relamination of subducted buoyant material (Hacker et al., 2011), shortening and imbrication of the forearc (Pearson et al., 2017; Sauer et al., 2017), or magmatism and intra-arc shortening, as documented in several Cordilleran arc systems (e.g. Umhoefer and Miller, 1996; Cao and Paterson, 2016; Shea et al., 2018). As previously discussed, it likely does not result from underthrusting of old, supracrustal materials from the retroarc. Imbrication of supracrustal forearc materials is equally unlikely given that the observed mantle-like δ18O signatures preclude the melting of near-surface rocks. Periodic relamination of more felsic crustal components removed by subduction erosion is a possible mechanism that could explain the arc thickening and differentiation trends observed during flare-ups (Fig. 4). Intra-arc shortening is also a viable option and is supported in part by documented garnet growth and increasing pressure and temperature during the most volumetric flare-up at ca. 78 Ma (Bollen et al., 2017; Rusmore and Woodsworth, 1994). Similar P-T increases in metamorphic host rocks are not recognized for other flare-ups, however, which suggests that additional mechanisms, such as magmatic underplating, are required.

In addition to temporal changes in chemical signatures suggestive of crustal thickening, the southern Coast Mountains batholith also records complex shifts in the loci of magmatism through time. It has been observed that plutons are generally younger in the eastern Coast Mountains batholith, indicating an overall progressive inboard migration of magmatism through time (Friedman et al., 1995; Gehrels et al., 2009; Cecil et al., 2018). Greater detail provided by our dense geochronologic coverage in the southern Coast Mountains batholith, however, reveals a more complex migration pattern in synch with geochemical variations in magmas. Three notable features of these data are: (1) the inboard advance of the leading edge of the arc (shown in the red dashed line on Fig. 9) appears punctuated; (2) landward advance and trenchward retreat of the arc roughly coincide with the timing of flare-ups and lulls, respectively; and (3) the active magmatic footprint of the arc appears to widen during flare-ups and narrow during lulls, as illustrated by the colored horizontal bars in Figure 9. For example, plutons emplaced during the Jurassic flare-up (161–148 Ma) occupy nearly 120 km of arc width, whereas those emplaced during an ensuing lull at ca. 148–130 Ma are found in a <50 km arc-perpendicular swath.

We note that in the central batholith, spatial changes in the arc in the Jurassic to Early Cretaceous are difficult to interpret. This is due to the fact that intrusions emplaced during that time likely reflect two distinct arcs: a western arc that was emplaced into the Alexander-Wrangellia terrane with no magmatic activity at ca. 140–120 Ma, which is juxtaposed against an eastern arc that was emplaced into the Stikine and other related inboard terranes with continuous Early Cretaceous magmatism (Gehrels et al., 2009). In contrast to the central Coast Mountains batholith, we suggest that Jurassic–Early Cretaceous magmatism to the south represents the southern continuation of the eastern arc based on the observations that: (1) Early Cretaceous magmatism is continuous; (2) magmatism appears to sweep inboard, then outboard, as it does in the eastern arc to the north; and (3) paleomagnetic data suggest that Jurassic–Early Cretaceous magmas likely intrude rocks of the Intermontane superterrane (Rusmore et al., 2013). Latest Cretaceous–Early Tertiary dextral displacements within the batholith, such as those proposed in Rusmore et al. (2013), would be sub-parallel to the batholith axis and so would not significantly disrupt arc-perpendicular geochemical trends. Significant (hundreds of kilometers) displacement would be needed to change the width of the arc via tectonic duplication. No such displacements can be documented on shear zones within the batholith, but they cannot be ruled out on the newly recognized Scar Creek shear zone located at ∼100 km distance from the coast, which was active between ca. 85 Ma and 60 Ma (Rusmore et al., 2019). Given the reasoning above, we interpret the spatial record of pluton emplacement as a single arc system that is advancing and/or widening with respect to the trench. Landward advances of the arc similar to the Jurassic event are recognized during the mid-Cretaceous (ca. 110 Ma) and Late Cretaceous (ca. 78 Ma) magmatic flare-ups (Fig. 9). Post-75 Ma spatial patterns in the southern Coast Mountains batholith are less well-defined, but it is possible that a similar advance accompanies the Paleogene (ca. 58 Ma) flare-up. We lack the data density necessary during any of these episodes to evaluate more detailed focusing or other spatial trends, as documented by Ardill et al. (2018) for the central Sierra Nevada.

Although a single mechanism is not required to explain each magmatic event, similarities in migratory patterns and arc geochemical signatures (depletions in heavy REEs) indicate that a similar process is periodically recurring in the southern Coast Mountains batholith and driving flare-ups. The punctuated inboard advance of the arc may trigger the observed mantle-dominated magmatic flare-ups by tapping into less depleted, more melt-fertile parts of the continental lithosphere (Fig. 10A). This enhanced melt production is attributable to more landward parts of the lithosphere (1) not having been previously melted, and/or (2) having been hydrated and (re)fertilized via the accumulation of metasomatic products from prior flux melting in the main arc, and/or (3) having been (re)fertilized via the underthrusting of older and more cratonic (igneous and quartzofeldspathic) inboard lithosphere (Chapman and Ducea, 2019). As previously discussed, the latter is likely of limited importance in the southern Coast Mountains batholith based on the primitive isotopic signatures observed.

Unanswered, of course, is the question of what causes the arc to advance in the first place. One possibility is that inboard advance is brought on by temporal variations in orthogonal convergence rate. An increase in convergence rate, if associated with increased plate coupling, would not only gradually shift magmatism inboard but could also lead to tectonic thickening in the arc and the increased bulk rock Sr/Y and La/Yb(N) ratios observed. Thickening could be enhanced by the addition of voluminous mafic underplates during the flare-up event. Increased convergence and plate coupling, if accompanied by a decrease in slab dip, could also explain a widening of the arc footprint. Orthogonal convergence rates for the Coast Mountain batholith region, as presented in Kirsch et al. (2016), vary considerably throughout the ca. 150 m.y. lifespan of the arc, from slightly negative (divergent) values early in the early arc (pre-160 Ma) to almost 200 mm/yr in Cretaceous time. Although the temporal variability is highly complex, increases in convergence rate broadly coincide with all but the most recent (ca. 58 Ma) flare-up in the southern Coast Mountains batholith. One problem with this mechanism is that magmatic tempo in the Coast Mountains batholith appears to change along strike, which would require spatially variable convergence rates at unlikely scales (∼150 km of arc length) (Cecil et al., 2018).

Another possibility is that arc migration is driven by non-steady-state forearc subduction erosion, as proposed in Jicha and Kay (2018). In this scenario, the periodic erosion of forearc crust drives the arc system landward, though the trench-arc distance remains the same. Relamination of forearc crustal components—particularly if they are of the mafic, isotopically juvenile Wrangellia terrane—could promote crustal thickening and supply source rock to the subarc melt zone that is compatible with the isotopic signatures observed in southern Coast Mountains batholith magmas. This scenario has the advantage of allowing for more along-strike variability in magmatic patterns, but it is difficult to evaluate given that much of the forearc to the Coast Mountains batholith has been removed by transcurrent faulting. In either case, coeval arc advance and magmatic flare-ups are observed and are linked to crustal thickening, trends which have been documented in other Cordilleran arcs (Haschke et al., 2002; Karlstrom et al., 2014).

Eventually, flare-ups in the southern Coast Mountains batholith wane, perhaps as fertile mantle depletes, or as the arc lithosphere overthickens and impinges on the subducting slab and cuts off asthenospheric circulation in the mantle wedge (Chin et al., 2012; Karlstrom et al., 2014), or as the hot mantle wedge is pinched out by slab flattening, which Haschke et al. (2002) argues leads to magmatic gaps in the Andean arc record. These waning stages appear to coincide with lower and more homogeneous bulk rock Sr/Y and La/Yb(N) (Fig. 4) and with trenchward retreat of the arc (Fig. 9). We speculate that these features result from the periodic removal of dense arc residue, which pushes the arc back and aligns it with thinned and/or depleted arc lithosphere (Fig. 10B). During this time, the arc enters a lull phase, and the depth of magmatic differentiation decreases. A return to normal supra-slab mantle convection and landward migration of the arc into the more fertile lithosphere is required to start the next flare-up cycle.

New zircon Hf and O isotopic data from Jurassic to Eocene plutons of the southern Coast Mountains batholith reveal mantle-dominated melt sources across a broad swath of the arc and over ∼120 m.y. of episodic arc activity. This suggests that magmatic episodicity is not driven by the periodic input of evolved, supracrustal materials in sub-arc melt zones but rather must be associated with spatiotemporal changes in mantle melt productivity. Our results support recent work that documents the significance of mantle melt contributions to continental batholiths (Decker et al., 2017; Martinez Ardila, 2019; Attia et al., 2020; Klein et al., 2021). This has important implications for the generation of continental crust in arc systems and implies significant new addition, and limited recycling, of crust along plate margins. Geochemical indices such as whole-rock Sr/Y and La/Yb, and zircon Ti and Yb concentrations, change as a function of magmatic flux, which suggests that although the mantle remains the dominant melt source, the thickness of arc lithosphere increases during flare-ups and decreases during lulls. A closer look at spatial trends in the available geochronology for the southern Coast Mountains batholith reveals periods of inboard arc advance and retreat, leading us to suggest that arc migratory trends fundamentally control magmatic tempo. In our preferred model, the arc advances inboard, accessing previously (re)hydrated and (re)fertilized lithosphere, which triggers mantle-like magmatic flare-ups as proposed in Chapman and Ducea (2019). These inboard arc advances are accompanied by thickening of the arc via shortening, magmatic underplating, relamination, or a combination thereof. Eventually, dense arc residue formed at the base of the arc column becomes gravitationally unstable and is convectively removed, which forces the trenchward retreat of the arc into regions where the sub-arc mantle lithosphere is absent or depleted. Such a shift produces lower volumes of melt and magmas with geochemical signatures that are consistent with generation and emplacement in tectonically thinned lithosphere. Overall, the results of this study highlight the fact that the growth of large-volume Cordilleran batholiths does not require significant recycling of preexisting crust and that magmatic tempos can be highly episodic if the mantle is temporally and spatially heterogeneous in terms of its degree of hydration and melt fertility.

This research was supported by National Science Foundation awards EAR-1347219 and EAR-1655152 to M.R. Cecil, EAR-1347375 to G.E. Gehrels, EAR-1347212 to M.E. Rusmore, EAR-1347341 to H.H. Stowell, and EAR-1338583 for support of the Arizona LaserChron Center. WiscSIMS is supported by the U.S. National Science Foundation (EAR-1658823) and the University of Wisconsin–Madison. The authors are grateful for insightful reviews by Scott Paterson and Jonathan Miller, which greatly improved the paper. They also gratefully acknowledge Mark Pecha and Nicky Geisler for their help in acquiring the zircon Hf and trace element data and Kyle Oye and Christine Truong for their help with whole-rock geochemical analyses. Don Willson of Silver King Ventures and Mike King and Jim Henderson of White Saddle Air Services provided critical field logistical support.

1Supplemental Material. Includes sample location information, whole rock geochemical data, and individual zircon trace element, Lu-Hf isotope, and O isotope data. Please visit to access the supplemental material, and contact with any questions.
Science Editor: Andrea Hampel
Associate Editor: Julie Roberge
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.