Tracking 40 Million Years of Migrating Magmatism across the Idaho Batholith Using Zircon U-Pb Ages and Hf Isotopes from Cretaceous Bentonites

: Cretaceous strata preserved in Wyoming contain numerous large bentonite deposits formed from the felsic ash of volcanic eruptions, mainly derived from Idaho batholith magmatism. These bentonites preserve a near-continuous 40 m.y. chronology of volcanism and their whole-rock and mineral chemistry has been used to document igneous processes and reconstruct the history of Idaho magmatism as emplacement migrated across the Laurentian margin. Using LA-ICP-MS, we analyzed the U-Pb ages and Hf isotopic compositions of nearly 700 zircon grains from 44 bentonite beds from the Bighorn Basin, Wyoming. Zircon populations contain magmatic autocrysts and antecrysts which can be linked to the main pulses of the Idaho batholith and xenocrysts ranging from approx. 250 Ma to 1.84 Ga from country rocks and basement source terranes. Initial ε Hf compositions of Phanerozoic zircons are diverse, with compositions ranging from − 26 to nearly +12. Based on temporal trends in zircon ages and geochemistry, four distinct periods of plutonic emplacement are recognized during the Mid- to Late Cretaceous that follow plutonic emplacement across the Laurentian suture zone in western Idaho and into western Montana with the onset of Farallon slab shallowing. Our data demonstrate the utility of using zircons in preserved tephra to track the regional-scale evolution of convergent margins related to terrane accretion and the spatial migration of magmatism related to changes in subduction dynamics.


Introduction
Bentonites, volcanic ash beds (airfall tuffs) devitrified and transformed into clay, are found globally and range in age from near-modern to Precambrian [1]. Airfall tuffs can be deposited over vast areas, including areas distant from the interiors of orogenic belts where the effects of thermal overprinting are absent. These distally deposited ash beds have the potential to serve as a robust source of geochronologic and geochemical data that track the evolution of their magmatic source terranes. Bentonite eruptive ages and whole-rock geochemistry have previously been used to identify the tectonic setting of coeval magmatic rocks in various settings across the globe (for example, [2][3][4]) and refine accretionary events at convergent plate boundaries [5] and more detailed study of bentonite phenocrysts has been used for stratigraphic correlation [6] and to infer magmatic processes [7]. Prior studies typically have sampled bentonites that formed during a limited time interval of only a few million years and so long-term changes in tectonic setting and the magmatic process that ultimately led to eruption of tephra have not been considered. In our study, we have analyzed a series of bentonites beds that span 40 million years of eruption.
The Cretaceous-Paleogene Idaho batholith, one of the largest in western North America and related plutons in western Montana that also have related volcanics [8], record magmatism from 125 to 55 Ma [9]. This volcanism is well-preserved over 500 km away in the Bighorn Basin, Wyoming and in western South Dakota by bentonites in Cretaceous strata of the Western Interior Seaway [10]. The whole-rock geochemistry and Sr and Nd isotopic compositions of these bentonites track changes in crustal thickness and the migration of the magmatic front across Idaho into Montana. By combining variations in crustal thickness and Nd depleted mantle model ages, magma source terranes have also been identified [11].
These same Cretaceous bentonites contain volcanic zircon, allowing analysis of their U-Pb ages and associated Hf isotopic compositions. While it is possible that each bentonite bed we sampled represents an amalgamation of more than one eruption of ash accumulated during periods of low sedimentation (e.g., [12]), our goal is to sample ash erupted over a time span of 10 s of millions of years to test models of large-scale shifts of magmatic activity across our larger study area, rather than analyze individual eruptions of ash. Within the zircon populations we analyzed, zircon age data in bentonites from all six formations we sampled are in accord with biostratigraphically-constrained ages [13,14]; these are autocrystic zircons that can be assigned to eruption of ash. There are also older zircons that are unlikely to be detrital in origin based on mineralogical characterization [15]: these can be xenocrysts, recording the ages of basement source rocks, assimilated country rocks, or older (unrelated) plutons through which younger magma was intruded and antecrysts that track early stages of magma crystallization related to the eruptions of ash that produced the bentonite beds.
Although it is possible that some grains may be detrital, few of them are rounded and previous analysis of the mineralogy of the bentonites we sampled showed the complete absence of detrital quartz and feldspar [15]. Zircon autocrysts, antecrysts and xenocrysts are also present throughout our different zircon populations. Most of these different types of zircon phenocrysts have maintained their isotopic signatures through the deposition of ash as tephra layers, any reworking during sedimentation and the transformation of tephra into bentonite. A small portion of the zircon grains we analyzed have anomalously young ages, attributed to hydrothermal or weathering processes.
Hf isotopic compositions can be used to constrain source rocks and in the case of the Cretaceous bentonites we sampled, test that magmatism migrated across Idaho into western Montana based on whole-rock Sr and Nd. We acquired 697 analysis points from 44 bentonite beds sampled from the Bighorn Basin, Wyoming and southwestern South Dakota (Figure 1), covering a forty-million-year span of eruption activity and report the ages of the different types of zircon phenocrysts and Hf isotopic compositions for the first time in these rocks. Our goals are to (1) determine a more precise chronology of magma migration associated with the formation of the Idaho batholith and (2) determine crustal inheritance using U-Pb geochronology of zircon xenocrysts combined with 176 Hf/ 177 Hf compositions and associated Hf depleted mantle model ages, the estimated age of crustal formation (source age) based on an assumed mantle source Lu-Hf composition.

Geologic Background
In the Cordillera of western North America, several extensively studied Cretaceous plutonic complexes, including the Sierra Nevada, Peninsular Range and Idaho batholith, record magmatism from the eastward subduction of the Farallon plate beneath North America. Emplacement of these plutonic centers began following a change in subduction polarity, resulting in an Andean-style continental margin stretching from Mexico to Alaska [16]. Widespread crustal shortening of the margin from dextral transpression [17] (largely synchronous with convergence [9]), magma intruding the upper plate and significant subduction erosion on the continental margin [18], resulted in the retrograde (landward) movement of the volcanic front through time. This migration has been documented through the abrupt eastward change of juvenile to evolved Sr, Nd, Hf and O isotopic compositions of plutons as emplaced magma encroached upon Laurentia basement [19][20][21]. In each of the major Cretaceous plutonic complexes, magmatism persisted for tens of millions of years producing plutons that overlap in time and space [22]. This long history of melt generation resulted in a complex amalgamation of plutonic lobes and associated local metamorphism [23]. Later periods of uplift and subsequent erosion have removed kilometers of upper crust, exposing the deep plutonic roots formed during subduction magmatism [24], similarly to the southern exposure of the Sierra Nevada batholith [25]. in north-central Wyoming and near the Angostura Reservoir (ANG) in southwestern South Dakota (see [11] for more detailed locations). Larger map shows simplified locations of the Idaho batholith, including the Southern Atlanta lobe (SAL), Northern Atlanta lobe (NAL) and the Bitterroot lobe (BR), in addition to the Boulder batholith (BB); the 0.706 Sr isopleth, the isotopic boundary separating Laurentia from exotic and accreted terranes now exposed in the Blue Mountains. Generalized boundaries drawn with dashed lines for various basement terranes. Other abbreviations: BMA = Blue Mountain Arcs; GFTZ, Great Falls Tectonic Zone; PB = Pioneer Batholith; SRSZ = Salmon River Suture Zone.

Geologic Background
In the Cordillera of western North America, several extensively studied Cretaceous plutonic complexes, including the Sierra Nevada, Peninsular Range and Idaho batholith, record magmatism from the eastward subduction of the Farallon plate beneath North America. Emplacement of these plutonic centers began following a change in subduction polarity, resulting in an Andean-style continental margin stretching from Mexico to Alaska [16]. Widespread crustal shortening of the margin from dextral transpression [17] (largely synchronous with convergence [9]), magma intruding the upper plate and significant subduction erosion on the continental margin [18], resulted in the retrograde (landward) movement of the volcanic front through time. This migration has been documented through the abrupt eastward change of juvenile to evolved Sr, Nd, Hf and O isotopic compositions of plutons as emplaced magma encroached upon Laurentia basement [19][20][21]. In each of the major Cretaceous plutonic complexes, magmatism persisted for tens of millions of years producing plutons that overlap in time and space [22]. This long history of melt generation resulted in a complex amalgamation of plutonic lobes and associated local metamorphism [23]. Later periods of uplift and subsequent erosion have removed kilometers of upper crust, exposing the deep plutonic roots formed during subduction magmatism [24], similarly to the southern exposure of the Sierra Nevada batholith [25].
The Idaho batholith and related plutons extending into western Montana represent the most inboard occurrence of Farallon-related plutonism, which lasted from the Middle Cretaceous to the Eocene [9]. Idaho plutonism began in western Idaho primarily in ac- Location maps for this study. Inset map shows location of bentonite beds collected from the Bighorn Basin (BHB) in north-central Wyoming and near the Angostura Reservoir (ANG) in southwestern South Dakota (see [11] for more detailed locations). Larger map shows simplified locations of the Idaho batholith, including the Southern Atlanta lobe (SAL), Northern Atlanta lobe (NAL) and the Bitterroot lobe (BR), in addition to the Boulder batholith (BB); the 0.706 Sr isopleth, the isotopic boundary separating Laurentia from exotic and accreted terranes now exposed in the Blue Mountains. Generalized boundaries drawn with dashed lines for various basement terranes. Other abbreviations: BMA = Blue Mountain Arcs; GFTZ, Great Falls Tectonic Zone; PB = Pioneer Batholith; SRSZ = Salmon River Suture Zone.
The Idaho batholith and related plutons extending into western Montana represent the most inboard occurrence of Farallon-related plutonism, which lasted from the Middle Cretaceous to the Eocene [9]. Idaho plutonism began in western Idaho primarily in accreted exotic terranes [20] west of the 0.706 Sr isopleth (the isotopic boundary between the North American craton and exotic terranes) and the Salmon River Suture Zone (SRSZ), the boundary between oceanic terranes to the west and Laurentian crust to the east, which now coincides with the Western Idaho Shear Zone (WISZ). The magmatic front migrated into Laurentia by the Cenomanian (100.5 to 93.9 Ma) where it developed the main body of the batholith, the Early Metaluminous/Border Zone Suite and the Atlanta Peraluminous Suite in central Idaho [26].
Relationships between individual plutons and original fabrics associated with Idaho batholith magmatism underwent multiple periods of overprinting and destruction. The SRSZ is an example of convergence-related deformation and movement along this boundary, producing significant alteration of the original magmatic fabrics [27,28]. Plutons along this boundary initially spanned nearly 100 km but experienced extensive right-lateral transpression along the WISZ in the Late Cretaceous, shortening the region to <10 km [28]. This ductile faulting combined with regional tectonic shortening caused magmatism to migrate into central Idaho, where the main phase of batholith activity was initiated during the Cenomanian [8] and persisted through most of the Late Cretaceous.
During the Campanian (~80 Ma), subduction of the Shatsky Rise Large Igneous Province caused the angle of subduction to change to sub-horizontal [29] resulting in an eastward shift of the magmatic front (e.g., [30]). A broadened region of pluton emplacement occurred throughout the Campanian, stretching across Idaho and into western Montana [31]. By the Maastrichtian, the magmatic front was volcanically active in Montana, forming plutonic complexes including the Boulder batholith (Elkhorn equivalent; [31,32]), the Pioneer batholith [33] and the Tobacco Root batholith [34,35].

Materials and Methods
Only one bentonite bed was sampled from each of the Cloverly, Mesaverde and Meeteetse formations from which we collected bentonites. In the other three formations, the Mowry, Frontier, and Pierre, we collected multiple bentonite beds and merged zircons from each bed into combined groups. Sampled zircons (see [11] for detailed bentonite information; stratigraphy and sampling locations in Figure 2) were analyzed via laser ablation witha Photon Machines Analyte G2 Excimer laser (λ = 193 nm; pulse width = 8 ns) for 10 s ablation bursts, using a He carrier gas on a 40 µm target ( Figure 3). The U-Pb data were acquired using a Thermo Element2 HR ICPMS with a dual mode secondary electron multiplier at the Arizona LaserChron Center. SLM zircon was used as the primary standard for mass normalization [36] (553.3 ± 13.5 Ma) and FC-1 (1096.2 ± 19.2 Ma)/R33 (417.9 ± 26 Ma) as secondary standards. Correction for common Pb was made using measured 206 Pb/ 204 Pb with the assumed composition of common Pb [37]. Many of our zircon ages (up to a third or so in sample populations) are concordant, but for analysis targets that are not, we report 206 Pb/ 238 U ages for those that are less than 1.1 Ga and 206 Pb/ 207 Pb ages for those that are older than 1.1 Ga [36]. Ages within each group of bentonite zircon (grouped by stage) were calculated using IsoplotR [38] as the weighted mean average of the 206 Pb/ 238 U ratio for grains younger than 1 Ga and the 207 Pb/ 206 Pb ratio for grains older than 1 Ga and rejected outliers using Chauvenet's criterion to assure normal distribution. We also calculated the weighted mean average age of each zircon population from each sampled bentonite bed using IsoplotR; these plots are in the Appendix. From these plots, we assessed which zircon grains are autocrysts (falling with the calculated 1-sigma error), which are likely antecrysts (falling above the calculated 1-sigma error), which are xenocrysts (rejected in the calculation of the weighted mean average age) and which grains record post-eruption effects (falling below the calculated 1-sigma error). In addition, the weighted mean average age of small sets of concordant zircons within each population were also calculated using IsoplotR to simplify reporting the ages of different zircon sub-populations.
Using a two-stage depleted mantle model age (Hf DMM ) for each analysis (from [11]), an age range for mantle extraction of the source rock can be determined that is dependent on the 177 Lu/ 176 Hf ratio of the continental crust used. While this measurement is not fully quantitative and is based on several assumptions (e.g., 176 Lu/ 177 Hf), it can provide firstorder constraints on minimum and maximum separation ages, as well as the magnitude of mantle versus basement contribution to the overall Hf signature. Table S1: Zircon U-Pb and Hf database; Table S2: Sample information; Figure S1: Weighted Mean Age (WMA) histograms.
Zircons were separated from a single bentonite bed within the Cloverly Formation whose biostratigraphic age is late Aptian-early Albian ( [53]; 115-110 Ma); only a small population of grains (n = 36) was recovered. Though the Cloverly Formation is a nonmarine fluvial succession, the zircons analyzed were not rounded, suggesting limited terrestrial transport. Within this small population, spot analyses on individual grains yielded a wide range of ages, from Cretaceous to Paleoproterozoic and more than half of the analyzed grains (23) Figure 3A, grain iv). One analyzed grain had a 206 Pb/ 238 U age that fell within the stratigraphic eruptive age range.
Zircons were separated from 25 bentonite beds, ranging in thickness from 20 to 120 cm, within the late Albian to early Cenomanian Mowry Formation ( [54,55]; 105-96 Ma), yielding a large population (n = 336) with a large proportion of grains with Cretaceous ages ( Figure 4B). One hundred and five grains (31% of the total) have concordant ages and every sampled bed has concordant zircons whose ages fall within stratigraphic eruptive age range of approximately 105-96 Ma. Zircons from the Mowry common show concentrically zoned cores and unzoned rims in CL; cores yield older ages than rims (for example, Figure 3B, grains ii, vi and viii). The population also yielded two slightly discordant grains with Paleoproterozoic 206 Pb/ 207 Pb ages of 1.7 and 1.82 Ga, one slightly discordant grain with a 206 Pb/ 238 U age of 369 Ma and one discordant grain with a 206 Pb/ 238 U age of 252 Ma ( Figure 4B). The discordant Mowry zircons lie on a chord (1.84 Ga and 252 Ma) whose lower intercepts fall within the stratigraphic eruptive age. Excluding the concordant grains, approximately 35% of the population has 206 Pb/ 238 U ages ranging from two to 35 million years younger than the stratigraphic eruptive age of the Mowry (for example, Figure 3B, grains iii, iv, v and vi).
Seventeen bentonite beds, ranging in thickness from 90 to 300 cm, were sampled from the late Cenomanian-Turonian Frontier Formation and equivalent Belle Fourche Shale whose stratigraphic age is approximately 89-96 Ma [14]. A population of 206 zircon grains were separated from Frontier bentonite beds: of these, 28 concordant ages fall within the stratigraphic age range (18% of the total population) and the 206 Pb/ 238 U ages of an additional 65 analyses also fall within the stratigraphic age range (together making up 59% of the total population; Figure 5). The weighted mean average age of the 28 concordant grains is 95 ± 0.3 Ma. Like zircons in bentonites from the Mowry Formation, zircons from Frontier bentonites have cores and rims seen in CL with older ages in cores, but some grains have more complex distribution of ages (for example, Figure 3C, grains ii, iv and vii).  Zircons were separated from 25 bentonite beds, ranging in thickness from 20 to 120 cm, within the late Albian to early Cenomanian Mowry Formation ( [54,55]; 105-96 Ma), yielding a large population (n = 336) with a large proportion of grains with Cretaceous ages ( Figure 4B). One hundred and five grains (31% of the total) have concordant ages and every sampled bed has concordant zircons whose ages fall within stratigraphic eruptive age range of approximately 105-96 Ma. Zircons from the Mowry common show concentrically zoned cores and unzoned rims in CL; cores yield older ages than rims (for example, Figure 3B, grains ii, vi and viii). The population also yielded two slightly discordant grains with Paleoproterozoic 206 Pb/ 207 Pb ages of 1.7 and 1.82 Ga, one slightly discordant grain with a 206 Pb/ 238 U age of 369 Ma and one discordant grain with a 206 Pb/ 238 U age of 252 Ma ( Figure 4B). The discordant Mowry zircons lie on a chord (1.84 Ga and 252 Ma) whose lower intercepts fall within the stratigraphic eruptive age. Excluding the concordant grains, approximately 35% of the population has 206 Pb/ 238 U ages ranging from two to 35 million years younger than the stratigraphic eruptive age of the Mowry (for example, Figure 3B, grains iii, iv, v and vi).  [51]) and 2.0-3.6 Ga (Wyoming basement [52]). Pie diagrams show the abundance of grains from each bin. Example inheritance data (bottom left) from the Idaho batholith [46]. Plot and age peaks generated using detritalPy for Python 3.x [56].
Seventeen bentonite beds, ranging in thickness from 90 to 300 cm, were sampled from the late Cenomanian-Turonian Frontier Formation and equivalent Belle Fourche Shale whose stratigraphic age is approximately 89-96 Ma [14]. A population of 206 zircon grains were separated from Frontier bentonite beds: of these, 28 concordant ages fall within the stratigraphic age range (18% of the total population) and the 206 Pb/ 238 U ages of an additional 65 analyses also fall within the stratigraphic age range (together making up 59% of the total population; Figure 5). The weighted mean average age of the 28 concordant grains is 95 ± 0.3 Ma. Like zircons in bentonites from the Mowry Formation, zircons from Frontier bentonites have cores and rims seen in CL with older ages in cores, but some grains have more complex distribution of ages (for example, Figure 3C, grains ii, iv and vii).
Twenty-two other analyses were concordant ages marginally older than the eruptive age, ranging from 148 to 97 Ma. Nearly 20% of the analyses have Cretaceous 206 Pb/ 238 U ages older than the stratigraphic eruptive age, including 7 grains within the Cenomanian,  [51]) and 2.0-3.6 Ga (Wyoming basement [52]). Pie diagrams show the abundance of grains from each bin. Example inheritance data (bottom left) from the Idaho batholith [46]. Plot and age peaks generated using detritalPy for Python 3.x [56].
Twenty-two other analyses were concordant ages marginally older than the eruptive age, ranging from 148 to 97 Ma. Nearly 20% of the analyses have Cretaceous 206 Pb/ 238 U ages older than the stratigraphic eruptive age, including 7 grains within the Cenomanian, 17 within the Albian and 11 within the pre-Albian Early Cretaceous (for example, Figure 3C, grains i and ii). The population includes five older discordant analyses, four of which have 206 Pb/ 238 U ages of 152, 157, 252 and 377 Ma and one with a 206 Pb/ 207 Pb age of 1.3 Ga; this Mesoproterozoic grain and the discordant Devonian grain define a chord whose lower intercept falls within the stratigraphic age range ( Figure 4C). Twenty-three analyses have younger 206 Pb/ 238 U ages than the stratigraphic age, ranging in age from the Upper Cretaceous to as young as the Lower Eocene.
Zircons were separated from four bentonite beds within the Campanian (84-72 Ma) Pierre Shale collected south of the Black Hills in South Dakota [15] and a total of 73 grains were analyzed at single spots in grain interiors excluding grain cores seen in cathodoluminescence. Nearly all these analyses have Campanian ages ( Figure 4D). Each sampled bed has Campanian concordant zircons, totaling 28 grains (38% of the total population); together these concordant zircons have a weighted mean average age of 81 ± 0.2 Ma. One other concordant grain has an older age of 86 Ma, the same as the 206 Pb/ 238 U ages of two other grains. All the other grains within this population have Campanian 206 Pb/ 238 U ages with two exceptions: an older Cretaceous age of 108 Ma and a younger Paleocene age of 67 Ma.
One bentonite bed from the late Campanian (76-72 Ma) Mesaverde Formation was sampled and single analyses of 32 zircon grains yielded much greater age variability compared to zircons separated from bentonites from the Pierre Shale. Eight zircons have concordant ages which define the stratigraphic eruptive age: these grains have a weighted mean average age of 73 ± 1 Ma. Seven other grains have older concordant ages: the weighted mean average age of two of them is 80 Ma, for another two, 96 Ma, both with a standard error of 2 m.y. and for a set of three grains, 1.75 Ga with a standard error of 25 m.y. The remaining analyses have a range of 206 Pb/ 238 U ages: seven grains are Campanian, three have mid-Cretaceous ages extending to the Aptian and two grains are Triassic (252 Ma) and Devonian (373 Ma); one grain is slightly younger than the stratigraphic eruptive age (71 Ma; Figure 3D, grain i). Two other grains have 207 Pb/ 206 Pb ages of 1.48 and 1.77 Ga ( Figure 3D, grain ii) and one zircon core has an age of 3.3 Ga ( Figure 3D, grain iii.).
The youngest bentonite collected is from the Maastrichtian-aged Meeteetse Formation; one bed was sampled which yielded a small zircon population (n = 29) and each zircon was analyzed at a single spot ( Figure 4F). Seven zircon grains are concordant and all of these have Cretaceous ages. Two grains give a weighted mean average age of 67 ± 2.5 Ma in accord with the stratigraphic eruptive age; six other grains have the same 206

Hf Isotopes
The initial εHf compositions of Phanerozoic zircons we analyzed are diverse, ranging from −26 to nearly +12 ( Figure 6A-C). The Cloverly zircon population primarily consists of grains ranging from 800-400 Ma. The εHf composition of these Paleozoic to Neoproterozoic zircons form two clusters: a radiogenic subset, ranging from approximately −4 to +10 and a non-radiogenic subset, ranging from −12 to −26. The radiogenic subset has a range of Hf DMM from Early Neoproterozoic to Mesoproterozoic (1353-677 Ma), whereas the non-radiogenic subset is primarily Paleoproterozoic to Neoarchean (2582-1779 Ma). The youngest Cloverly zircon (110 Ma), the approximate stratigraphic age, has a highly radiogenic εHf of +12, grouping with two data points from the Mowry. This cluster represents a Hf DMM age of approx. 300 Ma.
Zircons from the Mowry and Frontier bentonite populations contain a wide, continuous range of εHf compositions. Mowry zircons have an overall larger range of εHf than Frontier zircons, ranging from −23 to +12, compared to −19 to −3. The Hf DMM age of the Mowry populations varies from 2205-360 Ma, mid-Paleozoic to the Paleoproterozoic (Rhyacian), whereas the Frontier population ranges from 2047-1107 Ma, Mesoproterozoic to Paleoproterozoic (Orosirian). The zircon with the highest εHf is also the oldest Mowry grain, with a 207 Pb/ 206 Pb age of approximately 1800 Ma. Bentonite zircons from the Campanian-aged Mesaverde Formation and Pierre Shale also maintain a wide range of εHf values, from −26 to +8.6. Zircons from this population have a Hf DMM age range of 1623-2718 Ma, the oldest ranging population in this study. The Maastrichtian-aged Meeteetse Formation has a smaller εHf range than the Mesaverde, spanning from −16 to +2 with a range of Hf DMM ages of 2066-1565 Ma.

Zircon Autocrysts, Antecrysts and Xenocrysts and Migrating Magmatism
In the 49 bentonite beds we sampled, zircon populations contain autocrysts grown in the melts that were erupted to form the ash beds and the ages of these zircons fall within the stratigraphically constrained ages of the ash layers (with the exception of the one bentonite bed collected within the Cloverly Formation). The zircon populations also contain antecrysts [39] recycled from older igneous material of the same magmatic system of the bentonites and xenocrysts derived from the country rocks through which the erupted magma moved and crustal source rocks. Single zircon grains can have xenocrystic or antecrystic cores and antecrystic or autocrystic rims (for example, Figure 3C, grains ii, iii and vi) as well as autocrystic cores and rims with ages that are younger than the stratigraphic eruptive age (for example, Figure 3B, grains ii, vi and ix). The proportion of ages that are in accord with the stratigraphic eruptive age and that are older in each bentonite zircon population changes through time, recording differences in the tectonic setting of each set of ash eruptions as Cretaceous magmatism migrated spatially [4] across the western margin of Laurentia during Farallon plate subduction. These differences also reflect changes in magmatic processes related to crustal thickness and the depth of magma generation [4]: the late Albian and early Cenomanian, when Sr isotopes in the bentonites were homogeneous and their whole-rock trace element compositions record magma generation in a MASH zone in the thickest crust [11], is precisely when bentonite beds contain the largest proportion (40%) of autocrystic zircons.
The late Aptian-early Albian age of the Cloverly Formation and associated bentonites overlaps with the early stages of Idaho magmatism west of the 0.706 Sr isopleth, a position supported by bentonite whole-rock data [11] where a complex amalgamation of arc terranes and adjoining basins began to dock with Laurentia [7]. The Cloverly bentonite zircon population ( Figure 5) includes xenocrysts whose Mesoproterozoic inherited ages match those in the Northern Atlanta Lobe (NAL) of the Idaho batholith [46] whereas younger Early Cretaceous to Middle Jurassic inherited ages match parts of the Blue Mountain Province and associated island arcs [57]. Ordovician inheritance is consistent with magma being intruded through plutons produced following the Late Proterozoic breakup of Rodinia [13,58], a signal which is common in zircons from the Southern Atlanta Lobe (SAL) of the Idaho batholith [46]. The diversity of Devonian and Silurian zircon ages from both euhedral and sub-euhedral grains ( Figure 3) indicates a different source, likely consisting of significant volcaniclastics mixed with recycled detrital zircons from supracrustal sources. The fact that a detrital zircon population from a Cloverly Formation sandstone has similar Middle Paleozoic ages [13] supports some terrestrial input into the sampled bentonite bed, despite the fact that the zircons we analyzed are not rounded.
During late Albian-early Cenomanian time, magmatism migrated eastwards across the Salmon River Suture Zone (SRSZ), the boundary between accreted terranes and the western Laurentian margin, co-eval with transpression along the western Idaho shear zone (WISZ); [24,28] which now forms the western boundary of the Idaho batholith. This migration and deformation were coeval with eruption of ash now found as bentonite beds in the Mowry Formation. Magmatism during this time took place along the eastern region of the SRSZ, including the Idaho Border Zone Suite [17]. Xenocrystic zircons with Paleoproterozoic inheritance and the absence of Early Paleozoic and Mesoproterozoic grains in the bentonites from the Mowry Formation which we analyzed is evidence that the magmatic front had migrated to lie largely within Laurentian crust during the eruption of Mowry ash beds.
The migration of magmatism across the suture zone and into Laurentian crust led to magmatism centered within thick Precambrian crust where a sub-lithospheric crustal reservoir developed, homogenizing whole rock 87 Sr/ 86 Sr and developing enhanced levels of differentiation (high Eu anomaly and Zr/TiO 2 ; [11,15,59]). This was likely a mixingassimilation-homogenization-storage (MASH) zone, as magma pooled at the Moho and began a long period of supra-subduction calc-alkaline pluton formation. Much of the inheritance preserved in the small Cloverly bentonite zircon population is absent in the significantly larger Mowry population ( Figure 5), providing evidence that magmatism had largely migrated eastwards away from the suture zone by approximately 100 Ma. The large proportion of autocrystic zircons may reflect the dynamics of magma evolution where lower SiO 2 , Zr-richer mantle-derived melts mix with assimilated higher-SiO 2 , Zr-poorer melts in a system where magmas ascend steadily and those that erupt ash crystallize zircons at relatively low temperature in the shallow crust [60].
Zircons in the Mowry population with late Albian-early Cenomanian ages coincide with the emplacement age of the Early Metaluminous Suite (EMS) of the Idaho batholith [26]. Zircons with older Early Cretaceous ages between 115 and 125 Ma could be antecrysts from plutons now located along the WISZ [61] (Suture Zone Suite of [8]). The remaining zircons are anywhere from <1 to 18 million years younger than the eruptive age and plot discordantly. This fraction of the population is interpreted as a result of low temperature hydrothermal alteration and related Pb-loss in metamict outer zones [62] during burial to a probable depth of 2 to 3 km prior to Laramide deformation or due to weathering [63]. These metamict zones can be seen in CL images that show chaotic zoning, large secondary outer growth bands and cross-cutting textures (for example, Figure 3B, grain v).
The zircon population from the bentonite beds of the late Cenomanian-Turonian Frontier Formation is characterized by a large proportion (nearly half of the total population) of autocrystic grains. The zircon population also includes a significant proportion (36%) of antecrystic grains, including grains with concordant ages. These antecrysts have Cenomanian, Albian and pre-Albian Early Cretaceous ages. Cenomanian ages as old as 100 Ma coincide with the main emplacement phase of the Early Metaluminous Suite of the Idaho batholith in central Idaho [26]. Zircons with older Early Cretaceous ages, including sets of concordant grains with ages of 111 and 124 Ma, in Frontier bentonites beds closely match those in the Mowry population. Older grains include concordant and discordant Late Jurassic ages which link the Frontier population to arc terranes now juxtaposed with the westernmost plutons of the Idaho batholith. Inherited Precambrian zircons are absent from the Frontier with one Paleoproterozoic exception.
The crustal reservoir that began to form in Cenomanian time continued to produce volcanism through the Turonian. We interpret the increase of autocrystic zircons in the Frontier population as evidence that continuing mantle recharge of crustal magma reservoirs and further development of a MASH phase of magmatism resulted in an overall decrease in melt interaction with basement and country rocks, resulting in a decrease in xenocrystic material (e.g., [64]). Frontier bentonites also record homogenized 87 Sr/ 86 Sr, enhanced Eu anomalies and high Zr/TiO 2 [11]. During Cenomanian to Turonian time, the main magmatic bodies of the EMS in central Idaho were emplaced [26] as transpression along the WISZ continued to push magmatism eastwards into thick Precambrian crust.
The Campanian marks the beginning of significant change in the location and dynamics of magmatism across Idaho and Montana [32] attributed to flattening of Farallon plate subduction [65], the result of subduction of the Shatsky Rise Large Igneous Province [66,67]. This pushed the magmatic front further eastward and produced a broadened region of magmatism across Idaho and western Montana, subsequently causing crustal melting in the Southern and Northern Atlanta lobes and forming the Atlanta Peraluminous Suite of the Idaho batholith [24,26]. In contrast to the thick crust in central Idaho, western Montana is marked by thinner crust, likely resulting in decreased differentiation and an increase in mafic rock production, exemplified by the Elkhorn volcanic complex and a significant mantle component of Sr and Nd isotopes in intermediate plutonic rocks of the Boulder batholith [31].
Eruption of ash during the Campanian is first recorded by bentonite beds in the Pierre Shale with an eruption age of 81 Ma and approximately 8 million years later, in the late Campanian Mesaverde Formation ( Figure 4D,E). The zircon population in four bentonite beds of the Pierre Shale is composed nearly entirely of grains of Campanian age, with only a very minor proportion of grains with older ages, 86 and 108 Ma. Zircons from one bentonite bed in the Mesaverde define a more diverse population: 80 and 96 Ma, older Cretaceous extending to the Aptian, Triassic (252 Ma) and Devonian (373 Ma). Three xenocrystic zircons in the Mesaverde are concordant at 1.75 Ga and single grains have ages of 1.43 and 1.77 Ga; the core of one zircon yielded the oldest age of our entire study, 3.3 Ga. The eruptive ages of the bentonites in both formations coincide with ages of granite plutons of the Boulder batholith in western Montana [31] and the eruptive age of Pierre Shale bentonites matches the initiation of the Elkhorn Mountains volcanics [32]. By the time of ash eruptions leading to bentonites in the Mesaverde Formation, the volcanic front appears to have been restricted to Montana, forming the Boulder, Pioneer and Tobacco Root batholiths and the contemporaneous Elkhorn volcanic complex, all of which were primarily emplaced within the northwest-southeast trending Great Falls Tectonic Zone, the Paleoproterozoic suture between the Archean-aged Wyoming and Medicine Hat terranes. and a significant mantle component of Sr and Nd isotopes in intermediate plutonic rocks of the Boulder batholith [31].
Eruption of ash during the Campanian is first recorded by bentonite beds in the Pierre Shale with an eruption age of 81 Ma and approximately 8 million years later, in the late Campanian Mesaverde Formation ( Figure 4D,E). The zircon population in four bentonite beds of the Pierre Shale is composed nearly entirely of grains of Campanian age, with only a very minor proportion of grains with older ages, 86 and 108 Ma. Zircons from one bentonite bed in the Mesaverde define a more diverse population: 80 and 96 Ma, older Cretaceous extending to the Aptian, Triassic (252 Ma) and Devonian (373 Ma). Three xenocrystic zircons in the Mesaverde are concordant at 1.75 Ga and single grains have ages of 1.43 and 1.77 Ga; the core of one zircon yielded the oldest age of our entire study, 3.3 Ga. The eruptive ages of the bentonites in both formations coincide with ages of granite plutons of the Boulder batholith in western Montana [31] and the eruptive age of Pierre Shale bentonites matches the initiation of the Elkhorn Mountains volcanics [32]. By the time of ash eruptions leading to bentonites in the Mesaverde Formation, the volcanic front appears to have been restricted to Montana, forming the Boulder, Pioneer and Tobacco Root batholiths and the contemporaneous Elkhorn volcanic complex, all of which were primarily emplaced within the northwest-southeast trending Great Falls Tectonic Zone, the Paleoproterozoic suture between the Archean-aged Wyoming and Medicine Hat terranes. Zircons with older Cretaceous ages in Pierre and Mesaverde bentonites could be antecrysts matching suites of the entire Idaho batholith. Those in Mesaverde bentonites with Meso-and Paleoproterozoic ages could be xenocrysts indicating a new, rejuvenated magmatic source that assimilated fresh basement rock, for example the basement Wallace terrane (Selway, [51]) or the Great Falls Tectonic zone with which their ages match. How-ever, an alternative explanation is that Paleozoic and older grains in the Mesaverde at least are detrital since their ages match detrital zircons from Mesaverde sandstone units [13] which have prominent peaks in age distribution spectra at 1.4 and 1.7 Ga, considered to be derived from Belt Supergroup sediments.
The final phase of ash eruption included in our study is recorded by bentonites from the Maastrichtian-aged Meeteetse Formation (n = 29; Figure 5) whose zircon population is similar to that in the Mesaverde, but with a smaller number and more limited range of inherited grains with a similar and minor Paleoproterozoic inheritance. Nearly 30% of the zircon population has an age of 73 Ma, in accord with the stratigraphic eruptive age and matching the end-Cretaceous magmatic pulse in western Montana. Zircons with older Cretaceous ages between 96-84 and age of 110 and 132 Ma could signal inheritance from near the WISZ, but best match the restricted age range of the Bitterroot Lobe [46] and of the northeastern Idaho batholith [70]. During the Campanian and Maastrichtian, the reduced angle of Farallon plate subduction resulted in a trail of magmatism across Idaho that required several million years to sequentially dissipate (e.g., [30]). By the Maastrichtian, the magmatic front was primarily in western Montana.

Hf Compositions and Source Terranes
The inferences drawn from the ages of inherited zircons can be tested using zircon 177 Hf/ 176 Hf compositions ( Figure 6). In Figure 6A, sets of evolution lines demonstrate the range of potential crustal sources of the zircon grains. The upper, middle and lower 176 Lu/ 177 Hf value for each potential evolution line is 0.0037, 0.0115 and 0.0193, respectively ( 176 Lu/ 177 Hf mean and 1-sigma variation, [71]).
The relatively small Cloverly zircon population spans a wide range of potential crustal sources but closely matches both the age and Hf composition of the western radiogenic exotic terranes of the Blue Mountains [68] and the depleted composition associated with the Idaho batholith, specifically the slightly younger Early Metaluminous Suite [26] (Figure 6C). This late Aptian-early Albian period of volcanism coincides with the plutonic migration across the SRSZ, as magmas intruded a large swath of western Idaho before settling in central Idaho by the Cenomanian. Most of the Cloverly zircons are Early Paleozoic and Neoproterozoic in age, likely originating from volcanism associated with magmatism during the rifting of Rodinia beginning in Neoproterozoic time [46,58]. The radiogenic Hf composition of the youngest (110 Ma) zircons from the Cloverly confirms the involvement of exotic terranes with a Hf DMM of approx. 300 Ma. The wide range of Hf composition in the Cloverly zircon population is indicative of a diverse number of source areas, explained by the rapid migration of the volcanic front from juvenile exotic terranes (i.e., the Wallowa) to the marginal ensialic arc (i.e., the Olds Ferry; [48]) and then finally stabilizing in Laurentian crust east of the SRSZ.
Mowry and Frontier zircons contain a continuous range of εHf with a large variation, spanning the regional ranges of the Southern Atlanta lobe, the Northern Atlanta lobe and the Bitterroot lobe. This εHf continuum results in difficulty assigning a specific Hf DMM age and is likely the result of significant mantle input involved with the production of Idaho magmas. Zircons formed during the Early Mesozoic have εHf compositions that follow a general trend into the magmatic age, 105-90 Ma, of Mowry/Frontier volcanism ( Figure 6B). This trend, which spans from approx. 150-85 Ma, matches a Hf DMM of approx. 1.65 Ga, evidence that basement rocks of the Wallace/Selway basement were contributing to the main phase of Idaho magmatism. However, with the onset of major Idaho activity stabilized in Laurentian crust and development of a MASH-zone reservoir in Idaho, εHf compositions diverge widely away from inherited Hf DMM towards the radiogenic mantle input, explaining the wide range in values. In concordant zircons from the Mowry population that fall within the stratigraphic eruption age, εHf values define a shift from −12.5 to −8 from 105 to 95 Ma; εHf values in autocrystic concordant Frontier zircons which define an eruptive age of 95 Ma range between −18 and −4.5.
The zircon population from the single bentonite bed of the Mesaverde Formation defines a wide range of εHf, similar to the range seen in the Cloverly zircon population. A few grains lie within the 1.65 Ga trend and match the εHf composition of the Blue Mountains, Early Metaluminous Suite and Atlanta Peraluminous Suite, suggesting the continuation of Idaho volcanism into the Campanian, likely in the Southern Atlanta lobe. Zircons with the most negative εHf values follow a Lu/Hf evolution line towards an age of approximately 2.2 Ga, intersecting the population's subset of Paleoproterozoic zircon xenocrysts. This transition to the incorporation of an older source matches the geographic broadening of magmatic emplacement from the Paleoproterozoic Wallace/Selway terrane to the Great Falls Tectonic Zone suture in southwestern Montana. This suture formed during a Paleoproterozoic-aged collision between two major Archean-aged basement domains, the Medicine Hat terrane to the north and the Wyoming terrane to the south. The Meeteetse population is restricted to the 1.65 Ga trend, indicating that volcanism in Idaho persisted through the end of the Cretaceous, while volcanism in Montana narrowed away from the bordering Archean domains.

Conclusions
U-Pb ages of zircons from Cretaceous bentonite beds initially deposited as ash in the Bighorn Basin far to the east of the active magmatic front can be linked to the emplacement of plutons that formed the Idaho batholith and other plutons in Idaho and western Montana. Most bentonite zircon ages record eruptive ages, consistent with the stratigraphy, but many grains are inherited, recycled from older plutons (antecrysts) and basement source rocks (xenocrysts) and are important fractions of the zircon populations which are critical in interpreting provenance. Bentonite zircon εHf compositions provide another record which can be used to determine magmatic provenance. The model presented here is similar to the plutonic history presented previously [11], proving the robust nature of the original trace element and isotopic signatures of whole rock data and the correlative zircon geochemical data presented here. The significant difference between the two methods is the accuracy in which basement terranes can be identified. Trace element signatures can be reliably traced to relative variations in crustal thickness and juvenility, whereas zircons provide a robust age distribution to track specific basement associations. The absence of detrital quartz or feldspar grains found in the bentonite beds sampled [15] is key to the magmatic origin of the grains we analyzed. Even in the case of analysis of only one bentonite bed within a particular formation (for example, one bed from the Cloverly) containing only a relatively small zircon population (35 grains), useful data bearing on magmatic provenance can be obtained.
By integrating zircon U-Pb ages and Hf isotopic compositions-data from the same analysis target in individual zircon crystals-the history of volcanism within the Idaho batholith complex has been reconstructed as follows: 1.
The Aptian/Albian phase consisted of voluminous ash production across the Salmon River Suture Zone into the region of the Southern Atlanta Lobe of the Idaho batholith, pushed eastward by contemporaneous dextral transpression along the 0.706 Sr isopleth ( Figure 7A).

2.
After the magmatic front crossed into Laurentian crust in the Cenomanian, the primary phase of Idaho batholith development commenced. Magmatism spanned the entire N-S transect of the batholith across both the Northern and Southern Atlanta Lobes and retained an active geochemical equilibrium due to the development of a sub-lithospheric MASH zone. Eastward migration ceased as SRSZ transpression and associated uplift stabilized ( Figure 7B).

3.
By the Campanian, the Shatsky Rise Large Igneous Province subducted beneath Laurentia, resulting in the significant shallowing of Farallon subduction. This caused a large broadening of the active magmatic zone, pushing magmatism into western Montana with the onset of the Elkhorn Volcanic complex and the Boulder batholith and continued magmatism in Idaho ( Figure 7C).

4.
By the Maastrichtian, magmatism associated with Farallon subduction waned and Montana volcanism narrowed with the continued emplacement of the Boulder batholith and the Pioneer batholith. Minor activity also persisted in the South Atlanta lobe region of Idaho. This prolific, long-lived belt of magma production came to a halt with the onset of the Laramide orogeny and remained quiescent until the Eocene magmatic culmination ( Figure 7D).
sub-lithospheric MASH zone. Eastward migration ceased as SRSZ transpression and associated uplift stabilized ( Figure 7B). 3. By the Campanian, the Shatsky Rise Large Igneous Province subducted beneath Laurentia, resulting in the significant shallowing of Farallon subduction. This caused a large broadening of the active magmatic zone, pushing magmatism into western Montana with the onset of the Elkhorn Volcanic complex and the Boulder batholith and continued magmatism in Idaho ( Figure 7C). 4. By the Maastrichtian, magmatism associated with Farallon subduction waned and Montana volcanism narrowed with the continued emplacement of the Boulder batholith and the Pioneer batholith. Minor activity also persisted in the South Atlanta lobe region of Idaho. This prolific, long-lived belt of magma production came to a halt with the onset of the Laramide orogeny and remained quiescent until the Eocene magmatic culmination ( Figure 7D).  [59,60]. Simplified base map modified from [41].