Trace-element and Sr and Nd isotopic geochemistry of Cretaceous bentonites in Wyoming and South Dakota tracks magmatic processes during eastward migration of Farallon arc plutons

Bentonite beds, which are clay deposits produced by the submarine alteration of volcanic tephra, preserve millions of years of volcanic products linked to magmatic systems for which records are oth-erwise lost through erosion and alteration. Cretaceous strata from the Bighorn Basin, Wyoming, and southwestern South Dakota contain bentonites that originated from arc magmatism produced by subduction of the Farallon plate. We analyzed the bulk major- and trace-element geochemistry, and the 87 Sr/ 86 Sr ( n = 87) and 143 Nd/ 144 Nd ( n = 26) isotopic compositions of individual bentonite beds from these areas spanning 40 m.y. of volcanism to recover signals of magmatic processes and to attempt to trace bentonite geochemical and isotopic signatures to contemporaneous Cordilleran plutonic rocks. Using multiple immobile elements (e.g., Zr, TiO 2 , Nb, Ta, and rare earth elements), distinct temporal trends show variations in the effects of mineral fractionation and changes in crustal thickness. Bentonite Sr and Nd isotopic compositions allow ash beds to be correlated with specific batholithic complexes in Idaho and western Montana. With this data set, we observed the following: (1) The volcanic arc migrated across the 0.706 isopleth between 115 and 105 Ma; (2) between 105 and 95 Ma, magmatism stalled in central Idaho and was supported through significant MASH (mixing-assimilation-storage-homogenization) processing; (3) by 85 Ma, a shallowing subduction angle resulted in the eastward migration of the volcanic front into western Montana while volcanism in Idaho diminished; and (4) around 75 Ma, evidence of Idaho volcanism is lost. Montana plutonism with significant assimilation of radiogenic basement and regional centers of local magma emplacement (i.e., Pioneer batholith).


INTRODUCTION
Tracking the history and geochemical evolution of ancient magmatism, and the ways in which it is linked to plate-tectonic models, is often compromised by preservation. When magmatism persists in one region for long periods of time, younger plutons intrude and overprint older plutons, for example, in the Tuolumne suite of the Sierra Nevada (Coleman et al., 2004). Accreted arcs are typically metamorphosed, which affects their mineralogy, mineral chemistry, and even whole-rock composition (Jamieson et al., 1998). Ancient plutonic complexes lose their volcanic edifice during syn-to postorogenic exhumation and can be affected themselves by millions of years of erosion, which can remove kilometers of upper-and midcrustal material. These processes result in the limited preservation of plutonic roots, such as in the Idaho batholith (Gaschnig et al., 2017) and the Sierra Nevada batholith (Hill and Hill, 1975). Superimposed processes erase much of the history and chronology of regional arc magmatism, leaving behind an amalgamation of the deep roots of multiple generations of plutons.
However, airborne tephra deposits produced by explosive volcanism often exist in basins proximal to volcanic arcs, preserving an extensive record of magmatism within the stratigraphic record. These tephra deposits, known as bentonite, can preserve the original geochemistry of their source magmas (Hong et al., 2019;Huang et al., 2018;Jones et al., 2016;Knechtel and Patterson, 1956;Huff and Kolata, 1989;Christidis, 1998;Kiipli et al., 2008). This is particularly true for immobile elements (e.g., Ti, Zr, Th, and Nb) and rare earth elements (REEs; Zielinski, 1982). Preservation of the Sr and Nd radiogenic composition of tephra in whole-rock bentonite also occurs in some circumstances .
The aim of this study was to use the extensive tephragenic record of Cretaceous volcanism preserved from Wyoming to South Dakota in basins proximal to Cretaceous plutonic rocks to (1) document the geochemical makeup of ∼40 m.y. of Cretaceous tephra, now preserved as bentonite, and assess how it preserves source magmas; (2) attempt to trace bentonite geochemical and isotopic signatures to contemporaneous Cordilleran plutonic rocks, and (3) assess the viability of using these bentonites to track the emplacement, geochemical evolution, and migratory history of Cretaceous Cordilleran magmatism related to the eastward subduction of the Farallon plate beneath western Laurentia.

Construction of the Cordillera
The Mesozoic was a time of significant plate convergence as Pangea fragmented, resulting in extensive belts of arc magmatism by Early Jurassic time (Dickinson, 2004). The westward migration of North America initiated a period of island-arc obduction onto the passive western margin beginning in the late Permian and lasting through the Cretaceous (Jones, 1990). During the Middle Jurassic, the accretion of many suture zone, a Cordilleran/Laurentian boundary in western Idaho (McClelland et al., 2000). All North American arc accretion had finished by 160 Ma, marking the end of Cordilleran terrane accumulation. During the Early Cretaceous, Farallon subduction transitioned to a single-trench, Andean-style, eastward-dipping margin (Jones, 1990;LaMaskin et al., 2015). This was a reversal in subduction and resulted in a period of magmatic quiescence, marked by a distinct absence of plutonic rocks and preserved volcanic ejecta eastward of the plate boundary (Armstrong and Ward, 1993;Dorsey and LaMaskin, 2008;Mc-Clelland et al., 2000).
At ca. 125 Ma, magmatism across the western margin of the continent began a 40 m.y. interval of significant batholith formation from Mexico to Alaska (Armstrong and Ward, 1993;Dickinson, 2004). These large plutonic systems include the Omineca-Selwyn belt and Coast plutonic belt in Canada, the Sierra Nevada and Peninsular Ranges from Mexico to California, the Idaho batholith in Idaho and Montana, and isolated plutons across Washington and Nevada. All of these continental arcs record a general eastward migration of plutonic emplacement, with the oldest dated rocks on their western margins and plutons with younger ages to the east (Chen and Moore, 1982;Gaschnig et al., 2010Gaschnig et al., , 2017Krummenacher and Doupont, 1975). Geochemical trends, from sodic-rich to normal calc-alkaline compositions, follow a similar west-to-east transition (Bateman and Roddick, 1983). The emplacement of these batholiths was generally within accreted Cordilleran crust or along the transition zone between oceanic and continental crust. The Idaho batholith is a major exception. Its isotopic composition reflects a more evolved crustal signature, indicative of emplacement within Laurentian crust (Gaschnig et al., 2011).
During the Late Cretaceous, magmatism broadened eastwards following shallowing of the Farallon plate, possibly due to a buoyant response from the subduction of the Shatsky Rise large igneous province (Liu et al., 2010;Livaccari et al., 1981). The shallowing resulted in the termination of many magmatic provinces across the continent while pushing magmatism eastward into western Montana (Robinson et al., 1968) by 80 Ma. Beginning in the Late Cretaceous, magmatic centers characterized by high-K compositions formed in southwestern Montana (for example, the Absaroka Mountains; Meen and Eggler, 1987), and by the early Eocene, small intrusive bodies extended into central Montana, but by this time, magmatism had changed character, represented by mafic alkalic compositions (including lamproites and minettes; O'Brien et al., 1995). Concurrent with eastward Farallon subduction, the development of the Western Interior Basin, extending from central Canada southward to the Gulf of Mexico cratonward of the orogenic belt, continued through the combined forces of crustal loading and depression associated with tectonic convergence (Bertog et al., 2007;Miall, 2009). The basin was further enhanced by dynamic loading through mantle convection (Liu et al., 2014), and, combined with a global eustatic rise in sea level, these factors resulted in the Western Interior Seaway extending across much of North America by middle Cretaceous time. The Western Interior Seaway proved to be an effective host for the vast quantity of volcanic ejecta produced through the end of the Cretaceous. Deposits of windblown ash are commonly found in the shale deposits of the Late Cretaceous and can be regionally traced for hundreds of kilometers (Elder, 1988).

Preservation of Whole-Rock Geochemical Signatures in Bentonite
Volcanic ash alters rapidly when exposed to marine conditions, resulting in the restructuring of amorphous volcanic glass into clay minerals, dominantly montmorillonite, a 2:1 (T-O-T) swelling clay with Na + and Ca 2+ as the interlayer charge balance (Cuadros et al., 1999). Na and Ca are replaced by K + when devitrified ash undergoes burial and diagenesis (and anchimetamorphism; Huff, 2016), and montmorillonite recrystallizes into a mixed-layer illite-smectite clay. Elemental mobility during devitrification and diagenesis strongly influences most other major cations as well. The structure of the alteration products (e.g., montmorillonite, kaolinite, and zeolites) repurposes many major elements, and changes to bulk-rock chemistry are thus variable depending on the environment of ash deposition. Large ion lithophile alkali and alkaline group elements only fit within the interlayer space of the 2:1 clay unit cell, resulting in significant concentration changes of Ca, Na, and K, and, to a lesser degree, Rb, Sr, Cs, and Ba during the formation of montmorillonite, for example. The major oxide(s) composition of smectite clays (e.g., montmorillonite), when normalized to the parent volcanic glass composition, shows immobility of Si, Ti, and Al and significant mobility of Fe, Mg, Ca, Na, and K (Christidis, 1998).
Many studies have described the behavior of minor and trace elements during the alteration of volcanic glass to clay, documenting both the overall mobility and durability of the volcanic geochemical signature, and the way in which the original chemistry is retained in the final product (Christidis, 1998;Kiipli et al., 2008Kiipli et al., , 2017. Still, caution regarding the use of classic discrimination diagrams is necessary (e.g., Winchester and Floyd, 1977) because of the loss of certain trace elements through alteration. However, the retention of many immobile elements with minor variations is similar among sets of bentonite samples, making comparisons among them viable, and, in some circumstances, making direct comparisons with the source material possible (Jones et al., 2016).
In a study of Paleozoic bentonite from Baltoscandia, incompatible high field strength elements (HFSEs), including Al, Ti, Nb, Zr, Sn, Pt, Hf, Ta, and Th, were mostly immobile (Kiipli et al., 2017). In this same study, other elements commonly used for bentonite interpretation, e.g., Sc, V, Ga, Y, and REEs, showed only slight mobility. In an example of recent bentonite formation from the Aegean islands, Zr, Nb, V, and Ni were immobile, and Y (and thus the heavy [H]REEs) was similarly immobile, though overall REE mobility increased as devitrification progressed (Christidis, 1998). In a study of Miocene-aged rhyolitic bentonite from Colorado, no noticeable difference between the REE composition of glass and that of the montmorillonite was observed, but there was a slight enrichment in Th, Ta, Hf, and Al, and additionally in Ca, Mg, Sr, Sc, P, and Cr in the clay, which resulted from structural incorporation or adsorption to the clay surface (Zielinski, 1982). The light (L)REEs, Ti, Zr, Hf, and Ta were found to be immobile in alteration of bentonite preserved within alkaline lake sediments (Summa and Verosub, 1992), whereas greater mobility was observed in the alkalis and HREEs.
Despite these documented variations in element mobility, the transformation of tephra to bentonite did not significantly alter the original magmatic radiogenic isotopic compositions of the sequence of Cretaceous-aged bentonite that we analyzed. Sr isotopic compositions are distinctive between time periods of eruption and can be distinguished from the contemporaneous oceanic Sr composition associated with fossil material . Oceanic Sr concentration is ∼7 ppm, i.e., considerably less than bentonite (see Supplementary Table S1 1 ), and thus foreign Sr is not a significant component of bentonite composition. However, certain geologic conditions may influence the Sr composition. For example, Sr may be exchanged with seawater if porosity increases due to physical disturbance of the ash bed (Jones et al., 2014). Hydrothermal processes have also been documented as a significant influence on the Sr composition, such as that seen in midocean-ridge ophiolite sequences (Bickle and Teagle, 1992). Nd isotopic compositions show similar durability, bolstered by the immobility of Sm and Nd during ash alteration.  also demonstrated that the Nd isotopic compositions of Cretaceous bentonite preserved in South Dakota from the same series analyzed here were unique in both Cenomanian and Campanian beds, and when plotted as 87 Sr/ 86 Sr versus ε Nd , the Cenomanian samples fit within the compositional regime of the Idaho batholith during the same time interval (Gaschnig et al., 2011). The Sr and Nd isotopic composition of the Campanian ash showed mixing between mantle values and the Boulder (Montana) batholith composition (du Bray et al., 2012).

Sample Collection
Eighty-seven bentonite samples were collected from Cretaceous strata on the eastern and southern margins of the Bighorn Basin, from Billings, Montana, in the north to Thermopolis, Wyoming, in the south (Fig. 1). Bentonite samples from the Pierre Shale were collected south of the Black Hills, on and around the margin of the Angostura Reservoir in southwest South Dakota ( Fig. 1, inset). Many Bighorn Basin samples were collected from open-pit mines, prospective mining localities, and sets of outcrop sections using published stratigraphy (all courtesy of Wyo-Ben, Inc.; Fig. 2). The mines often exposed entire bentonite layers, revealing deposits up to 2-3 m thick. Only bentonite beds thicker than 10 cm were sampled. To avoid weathering rinds, material to a depth of ∼20 cm was removed, and sample material was collected from the center of the deposit when possible. To assure that the material was not subjected to secondary diagenesis (conversion of smectite to illite), X-ray diffraction was performed to identify the clay mineralogy. All samples analyzed were composed of unaltered smectitic clay, indicative of deposition within a normal marine environment .
The samples are from the Late Cretaceous stratigraphic record, spanning ∼40 m.y. of volcanism ( Fig. 1). From oldest to youngest, the sampled bentonite beds are from following formations and named sections: the Cloverly Formation, the Thermopolis Shale, the Mowry Formation, the Frontier Formation, the Cody Shale, the Pierre Shale, the Mesa Verde For-mation, and the Meeteetse Formation. The Cloverly Formation is the basal Cretaceous unit in the Bighorn Basin and represents a fluvial/estuarine/tidal flat sequence with bentonitic mudstones (Finn, 2010). The Thermopolis Shale is a 70-m-thick Albian-aged marine deposit composed of black shales and siltstones with thick interbedded bentonite (Eicher, 1960). The Mowry Formation is late Albian to Cenomanian and is composed of a 120-m-thick siliceous shale unit with many bentonites ranging from 1 cm to several meters thick (Davis, 1963;Finn, 2010). The Frontier Formation is composed of marginal marine deposits of sandstone/siltstone interbedded with shale, containing abundant bentonite up to several meters thick (Kirschbaum and Roberts, 2005). The Cody Shale is a 1200-m-thick marine shale unit (Hagen and Surdam, 1984) that extends into the Campanian and contains several prominent bentonites; outcrop sections rarely expose the bentonite, but mining has revealed several prominent units. The Campanian-aged Pierre Shale is a marine unit, consisting of over 300 m of homogeneous black shale (Crandell, 1958). The Pierre Shale bentonites are up to 40 cm thick. The Mesaverde Formation is contemporaneous with the Pierre Shale but is characterized by a transitional marine depositional environment (Finn, 2010;Severnz, 1961) and is composed of marginal marine deposits, including sandstone, siltstone, mudstones, and coal, with several large bentonites found within the transgressional sequences. Last, the Maastrichtian-aged Meeteetse Formation consists of sandstone, siltstone, shale, and coal, and it was deposited in a coastal setting as the seaway began to diminish (Finn, 2010;Keefer, 1961). Several prominent bentonite units, 100-200 cm thick, within the Meeteetse Formation are present within deposits of carbonaceous shale.  presented a more detailed explanation of the formations sampled.

Bulk Geochemistry
The bulk sample was first homogenized using a pulverizer and agate mortar and pestle and then fused using a lithium metaborate/tetraborate fusion technique. Geochemical analysis was completed at Actlabs, Ancaster, Canada, using a combination of inductively coupled plasmaoptical emission spectrometry (ICP-OES) and inductively coupled plasma-mass spectrometry (ICP-MS). Standards NIST 694, DNC-1, GBW 07113, W-2, SY-4, BIR-1a, and NCS DC73304 were used for analysis calibration. See Hall and Plant (1992) for an overview of methodology precision. All geochemical data are presented in Supplementary Table S1.

Sr Isotopes
The bulk bentonite (n = 87; 0.5 g) was digested using 8 M nitric acid for ∼12 h to ensure complete clay dissolution. The minor phenocryst fraction was left behind for Sr analysis; there was no significant observed difference between the Sr composition of the bulk bentonite versus that of the isolated clay fraction . The Sr from the solution was then isolated using column separation with Eichrom Sr-spec resin. All samples were checked for Sr concentration using a single-collector ICP-MS, and the isotopic ratio was measured with a multicollector ICP-MS at the University of Illinois-Urbana-Champaign. The international standard NBS 987, "Coral," and "E&A" were used for correction and have certified 87 Sr/ 86 Sr ratios of 0.710255, 0.70918, and 0.70804, respectively. Mass bias was mitigated by dilution of strontium concentrations with 2%-3% nitric acid to maintain an 88 Sr intensity of 8.5-11.5 V. A Sr standard using SRM 987 was used with sample repeats to check data consistency. The 87 Sr/ 86 Sr composition was then corrected for Rb-decay using the approximate age of the bentonite stratigraphic position. All Sr isotopic data are presented in Supplementary Table S2 (see footnote 1).

Nd Isotopes
Twenty-six of the 87 bentonites were analyzed for their Nd and Sm isotopic compositions. A multicollector Triton Plus thermal ionization mass spectrometer was used for Nd analysis. The 147 Sm/ 144 Nd value was determined using isotope dilution with a 149 Sm-150 Nd standard. Approximately 100 mg aliquots of sample were dissolved using HNO 3 and exposed to microwaves. The REEs were separated by cation exchange using Re resin before the isolation of Nd and Sm. Ln resin was used to isolate Nd and Sm following the methods of Míková and Denková (2007). The isotopic ratio was normalized to 146 Nd/ 144 Nd = 0.7219 (O'Nions et al., 1977). Measured values in the international standard sample JNdi-1 (Tanaka et al., 2000;147 Sm/ 144 Nd = 0.512115) were 143 Nd/ 144 Nd = 0.512095 ± 10 (2SD, n = 5). For calculation of initial ε Nd values, 143 Nd/ 144 Nd CHUR = 0.512630 and 147 Sm/ 144 Nd CHUR = 0.1960 (Bouvier et al., 2008) were used, where CHUR is chondritic uniform reservoir. All Nd isotopic data are presented in Supplementary Table S2.

Depleted Mantle Model Age
A depleted mantle model (DMM) age was determined using the calculated 143 Nd/ 144 Nd isotopic composition for each sample. The model used a two-stage melt separation history: The first stage was a crustal component derived from the depleted mantle, followed by a second (much later) stage of melting to produce the measured volcanic material. The second stage (from the Cretaceous melting event to the present) used the measured Sm/ Nd ratio, whereas the first stage (from the depleted mantle melt to the Cretaceous melting event) used the average continental crustal (ACC) value of 0.12. Although the crust involved in Cretaceous-aged melting was likely far from homogeneous, and there is uncertainty with respect to the timing when it separated from the mantle, our assumptions allowed us test broad relations between the bentonite Nd isotopic compositions and those of the various basement provinces found across western North America during the Cretaceous. This two-stage approach has been used for felsic igneous rocks to account for changes in Sm/ Nd ratios due to magmatic processes, such as magma mixing, partial melting, and fractional crystallization (e.g., Liew and McCulloch, 1985;Champion, 2013). This two-stage model follows the depleted mantle (DM) model equation of Liew and Hofmann (1988, their

Volcanic Classification and Bulk Major-Element Geochemistry
Use of the trace-element ratios to classify the bentonites (Figs. 3A and 3B) resulted in a range of rock types, including trachyandesite and trachyte, and a highly variable alkalinity index Nb/Y (Winchester and Floyd, 1977), likely the result of Y mobility (Christidis, 1998). The analyzed samples showed a grouping with respect to age in the differentiation index Zr/TiO 2 , where the Aptian and Maastrichtian samples had overlapping, lower values compared to Albian and Cenomanian samples. For all samples, the rare earth element (REE + Y) concentration exceeded 100 ppm, but it rarely exceeded 400 ppm (with one sample in the Aptian [>800 ppm] as the exception), LREEs ranged from 20× to 500× chondrite values, a variable negative Eu anomaly was observed, and HREEs ranged from 5× to 100× chondrite values. Taken together, these Cretaceous bentonites show a REE signature typical of evolved magmas. When comparing the REE signature of the bentonites within a particular unit or across ages, they can show discrete compositional differences.
All analyzed samples showed major-oxide concentrations typical of differentiated felsic magmas (dacite to rhyolite), even though the transformation of volcanic ash to bentonite through the formation of clays has likely altered the relative weight percentages of the major oxides (Fig. 4A). Cations that occupy the interlayer site in clays showed the greatest variability compared to cations in the octahedral and tetrahedral sites. Octahedral and tetrahedral Al and tetrahedral Si showed the least variation and had mean values relative to continental crust of nearly one. Bulk sample SiO 2 ranged from ∼60 to 80 wt% with no clear trend with age. The concentration of Al 2 O 3 ranged from 10 wt% to as high as 30 wt%. A negative correlation was observed between SiO 2 and Al 2 O 3 (R 2 = 0.77) and, to a much lesser degree, MgO (R 2 = 0.28), an artifact of aluminum and magnesium substitution for silicon in the octahedral layer of the clay minerals. A negative correlation was observed between CaO and SiO 2 , and a positive correlation was observed between CaO and Al 2 O 3 , another artifact of clay formation. All samples also had total [Na + K + Ca] concentrations ranging from ∼0.5 to 10 wt%, resulting in a nominally highly peraluminous composition. However, the progressive formation of clay minerals removes [Na + K + Ca], and the low concentrations of the alkalis (K 2 O ranged from 0.1 to <3 wt%, Na 2 O ranged from 0.1 to <6 wt%) are significantly less than an average evolved granitic melt. Though only slightly correlated, the concentrations of Na

Trace Elements
In detail, trace element ( Fig. 4) and REE (Fig. 5) ratios show distinctive signatures with age, although not in a simple age progression. The Eu anomaly is a key to interpreting the role of plagioclase fractionation within magma and provides an approximation of the magnitude of differentiation (e.g., Belousova et al., 2002). All bentonites showed variable magnitudes of the Eu* anomaly, Eu/Eu* Eu / Sm Gd / n n n = + [( ) ] 2 , which ranged from 0.1 to 0.8 (Fig. 4B). The Albian and Cenomanian bentonites had the lowest values of Eu/Eu*, ∼0.14, ranging as high as 0.60, and with mean Eu/Eu* = 0.33 and 0.35, respectively. Values for Aptian, Campanian, and Maastrichtian bentonites were, with a few exceptions, higher than 0.6 and reached values as high as 0.78. The mean values of Eu* of Aptian, Campanian, and Maastrichtian samples were = 0.63, 0.66, and 0.71, respectively, i.e., over double the mean values of Albian and Cenomanian samples. When Eu* is plotted against the differentiation index Zr/TiO 2 , there is a negative linear relationship (R 2 = 0.72), in which the most differentiated samples (more felsic bentonite with higher Zr/ TiO 2 ) have lower Eu* values, whereas less differentiated bentonites (more dacitic with lower Zr/TiO 2 ) have higher Eu* values (Fig. 4D). This relationship appears to be primarily the result of Ti distribution rather than Zr: Eu* plotted against Zr is nonlinear, whereas the Ti relationship is highly linear (R 2 = 0.85).
The Dy* anomaly, Dy/Dy* = Dy/(La4/13 * Yb9/13) (Davidson et al., 2013), varied based on the existence of amphibole in the restite (Davidson et al., 2013), effectively approximating the depth of magma generation. There is a variation through time in the total range of Dy/ Dy* (Fig. 4F). Albian and Cenomanian samples showed the greatest range of values, between 0.41 and 0.94, with mean Dy/Dy* = 0.67 and 0.62, respectively. Most values of Dy/Dy* for Aptian, Campanian, and Maastrichtian samples were between 0.5 and 0.7, with mean Dy/Dy* = 0.58, 0.57, and 0.50, respectively. The ratio Dy/Yb (Fig. 4G) varied in parallel with Dy/Dy*, with Albian and Cenomanian samples having the largest range, 1.1 to as high as 4, and values for Aptian, Campanian, and Maastrichtian samples at the lower end of this range, typically 2 and less. Thus, Dy/Dy* and Dy/Yb showed a strong positive correlation, with Albian and Cenomanian samples dominating the trend toward the highest values. Albian and Cenomanian bentonites had the most significant HREE depletion, whereas Aptian, Campanian, and Maastrichtian samples had no noticeable HREE depletion (Fig. 4G).
The ratio of Nb/Ta is an important indicator of fractionation (Li et al., 2017), and both Nb and Ta are relatively immobile during glass alteration (Christidis, 1998), resulting in a useful proxy for bentonite provenance (Fig. 4C). The mean Nb/ Ta value was highest in the Aptian and Maastrichtian samples, with ratios of 14 and 14.3, and maximum values of 23.3 and 20, respectively. Albian and Cenomanian samples had mean values of 8.4 and 7.1 and maximum values of 12.2 and 12.5, respectively. Campanian samples had the largest range of values, with a maximum value of 20, and a mean value of 10.2.

Sr/ 86 Sr Composition
The 87 Sr/ 86 Sr composition of all samples (Fig. 6A), with values ranging from 0.70461 to 0.71101, indicated a significant crustal component. The Aptian bentonites contained the lowest measured value, 0.70461, the only measured value lower than the isotopic boundary value of 0.706. The mean isotopic value of the Aptian samples was near 0.706, followed by an increase to 0.70725 from the Aptian to the Albian bentonites. The Albian samples were consistently more radiogenic, with values ranging from 0.70708 to 0.71021 and a mean of 0.7081. The Cenomanian samples were similarly radiogenic but had a narrower range, 0.70743-0.70858, with a comparable mean of 0.708.
The Campanian samples, which included the Cody Shale, Pierre Shale (South Dakota), and Mesaverde Formation, showed lower 87 Sr/ 86 Sr values, ranging between 0.70625 and 0.70992, with a mean value of 0.70761. However, since bentonites of Campanian age extend over a larger geographic area and several stratigraphic units, their Sr isotopic values showed some variation within the stage. Bentonite from the Cody Shale ranged from 0.70625 to 0.70992, with a mean value of 0.70826; Pierre Shale bentonites showed a downward shift of the mean to 0.70701, the lowest of all Campanian units, with values ranging from 0.70642 to 0.70769. The Mesaverde Formation bentonites had a very limited range of isotopic composition, ranging from 0.70823 to 0.70853, with a mean of 0.70844.
Finally, the youngest samples from the Maastrichtian Stage had a range of 0.70967-0.71101, with a mean of 0.71057, the highest 87 Sr/ 86 Sr compositions. Taken together, the Campanian-Maastrichtian bentonites followed a similar trend with age as the Aptian-Cenomanian samples, with lower 87 Sr/ 86 Sr compositions in older samples, ranging to higher mean values through time.

ε Nd Composition and Depleted Mantle Model Ages
The 26 analyzed samples showed a range of Nd isotopic values (Fig. 6B). The Albian-aged bentonites had a ε Nd composition ranging from −10.8 to −8.9, with a mean value of −9.9. The Cenomanian bentonites displayed a range of ε Nd values extending from −10.4 to −4.0, with an average of −8.4. The Campanian bentonites were all from the Pierre Shale; these samples had ε Nd values ranging from −10.4 to −2.6, with an average of −6.0. The Campanian samples had the largest isotopic diversity and were the most radiogenic of the three stages analyzed.
Using the two-stage depleted mantle model of Liew and Hofmann (1988), model ages for each ε Nd composition were calculated assuming that (1) the basement of the Laurentian craton represents an initial separation from depleted mantle and (2) Cretaceous plutonism melted this lithosphere. However, in addition to Laurentian crust being composed of a heterogeneous amalgamation of plutons that represent a range of mantle separation ages, the terranes in question are poorly exposed, resulting in considerable uncertainty  (Hammarstrom, 1982).
in the extent of specific plutons throughout the region. Due to this variability, the calculated DMM age is best used by fitting the samples within age bins of candidate terranes. The Sm/Nd ratio for each bentonite that was used in the DMM calculations ranged from 0.14 to 0.25, with an average of 0.2. The model age for the Albian-aged bentonites was the most restricted and fell within the mid-Paleoproterozoic, extending from 1961 Ma to 1791 Ma, while the Cenomanian-aged samples gave model Paleoproterozoic ages as old as 1892 Ma and extended to a younger crustal age in the Mesoproterozoic of 1362 Ma. Three Campanian-aged bentonites yielded a range of ages from Paleoproterozoic to Mesoproterozoic, extending from 1859 Ma to 1209 Ma. The mean DMM age of each set of analyzed bentonites decreased through time, from 1878 Ma in the Albian, to 1739 Ma in the Cenomanian, to 1494 Ma by the Campanian.

DISCUSSION
The major-element chemistry of all the bentonites we analyzed demonstrates that the transformation of ash into bentonite and the formation of clays control changes to major-element abundances and ratios, and so they cannot be used for unequivocal volcanic classification. To classify bentonites, trace elements are required, but even Nb/Y, for example, can be affected by the formation of clays during devitrification. Nevertheless, the trace-element composition of the bentonites we analyzed is a more accurate representation of the original magmatic signature. In our sample set, when various diagnostic trace-element ratios are plotted against stratigraphic age, distinctive shifts in the magmatic signature are evident. The REE and trace-element compositions of the analyzed bentonites provide evidence that these erupted magmas record the variation in crustal thickness of the leading Laurentian crustal margin during the Cretaceous. This is supported by overlaying plutonic data from Idaho and Montana onto the range of REE values of each stage sampled. The Aptian (Fig. 5A), Albian (Fig. 5B), and Cenomanian (Fig. 5C) bentonites show a mixture between arc crust and Proterozoic crust, a result of the incorporation of diverse crustal sources across the Western Idaho suture zone (Kuntz and Snee, 2007). The Campanian bentonites (Fig. 5D) mirror the REE data from the Boulder batholith (du Bray et al., 2012), while the Maastrichtian bentonites (Fig. 5E) are similar to both the Pioneer (Hammarstrom, 1982) and Tobacco Root batholith(s) (Sarkar et al., 2009), plutonic sources proximal to the Boulder batholith in southwestern Montana.
The thickness of evolving suprasubduction arc lithosphere is thought to be directly related to the progression of differentiation in silicic magmas, and thus it correlates with indices of differentiation (e.g., Zr/TiO 2 ) and associated signatures of mineral fractionation (e.g., Eu/Eu*). The connection between crustal thickness and magma evolution has been well documented using an analysis of a global data set (Farner and Lee, 2017), showing that as lithospheric thickness increases, magma becomes, on average, more silicic and enriched in incompatible elements. The data set of Farner and Lee (2017) was potentially biased, however, due to regional oversampling, but this was corrected for by averaging compositional data over discrete regions. This phenomenon is associated with increased transit time from the mantle wedge to the upper crust and longer cooling periods, allowing for enhanced fractionation of plagioclase, pyroxene, and amphibole, which dominate equilibrium assemblages in evolving intermediate magmas (Grove and Brown, 2018). However, modeling studies of the mass transfer of elemental components that undergo chemical reaction with partial molten rock during melt migration offer another mechanism for chemical differentiation in magmas (Solano et al., 2012(Solano et al., , 2014. Here, too, crustal thickness plays a role, since differentiation in these models is proposed to largely occur at the base of thick, hot crust.

Evidence of Mineral Fractionation
The Albian-and Cenomanian-aged tephra deposits contain the most pronounced negative Eu anomalies (Fig. 4B). In both the older Aptian samples and in Campanian and Maastrichtian samples, negative Eu anomalies are much less pronounced, and in some samples, they are nearly absent. The shift in Eu depletion could indicate a shift in the relative amount of plagioclase fractionation through time, thus correlating with differentiation (e.g., Zr/TiO 2 ; Fig. 4D), but B A

Figure 6. Radiogenic Sr and Nd compositions of bentonite samples through the Cretaceous. (A) 87 Sr/ 86 Sr vs. age (n = 87). Bentonites of Aptian age shifted from low to high values as volcanism moved across the 0.706 isopleth. Note stability of values through the Cenomanian and renewed downward shift and spread of values in the Campanian. Bentonites of Maastrichtian age have the highest values. (B) ε Nd(i) vs. age (n = 26). Radiogenic composition increases from Albian-to Cenomanian-aged bentonites, and Cenomanian and Campanian bentonites range from low to high values. The highest Campanian-aged value coincides with the initiation of the Elkhorn volcanics in Montana, which are enhanced in mafic trace-elemental signatures. See text for discussion.
other evidence of plagioclase crystallization is lost from large ion lithophile element mobility (e.g., Ca, Sr) and cannot be correlated with Eu/ Eu*. Titanite fractionation, with partition coefficients of REEs ranging from 1 to 10 crystal/ melt (Tiepolo et al., 2002) to nearly 1000 crystal/melt in the middle (M)REEs (Glazner et al., 2008), could account for the lessening negative Eu anomaly from Cenomanian to Maastrichtian time, counteracting Eu depletion associated with plagioclase removal. Fractionation models show that a small increase in titanite crystallization (99.69 wt% plagioclase/0.31 wt% titanite) can cause the Eu/Eu* to go from negative to positive (Loader et al., 2017) due to high MREEs and low Eu sequestration in titanite, regardless of the significantly larger proportion of plagioclase. The positive correlation in our analyzed bentonites between Eu/Eu* and TiO 2 may be evidence of the role of titanite in the sequestration of Eu from the magma during the eruption of Cretaceous tephra. In addition, there is an indication that crustal thickness (pressure) has played a role in suppressing plagioclase saturation: Campanian and Maastrichtian bentonites consistently have the highest Eu/Eu* values and low Dy/Yb (Fig. 7); Albian-Cenomanian bentonites have much greater variability. Taken together, variation of the Eu/Eu* signal, Eu/Eu* versus TiO 2 correlation, and changing MREE depletion likely mean that the fractionation of plagioclase played a role in controlling REE distribution, likely buffered by variations in the crystallization of both amphibole and titanite in the magma.
Changes in the degree of melt differentiation are indicated by bentonite Nb/Ta ratios. There is a separation of Nb/Ta ratios through the bentonite chronology, with the Nb/Ta of ∼12 acting as a dividing point. Most of the bentonites had Nb/ Ta ratios lower than the mid-ocean-ridge basalt (MORB) value of 15.4 (Gale et al., 2013), but the values within Aptian, Campanian, and Maastrichtian ash were close to the MORB value, and sometimes exceeded it (Fig. 4C). This seems to demonstrate that, in the many of the samples, very little Nb has been fractionated from Ta, showing either a lack of a Nb fractionating phase (such as rutile) or less crustal differentiation. It should be noted that rutile can have a significant influence on Nb/Ta while not readily accepting REEs into its crystal lattice (Foley et al., 2000). With progressive differentiation, a decrease in magma Nb/Ta occurs during rutile crystallization (Tang et al., 2019). While garnet, amphibole, and pyroxene also fractionate Nb over Ta, their partition coefficients are lower than 1 (Klimm et al., 2008;Skora and Blundy, 2010) and do not have a significant effect on the overall Nb/Ta composition, though this may change in highly silicic melts. Additionally, Nb is generally concentrated in within-plate granites (Pearce et al., 1984), and the Atlanta lobe of the Idaho batholith shows unusually high Nb concentrations (Gaschnig et al., 2011). Bentonites erupted through the Albian/ Cenomanian transition are characterized by low (<3) La/Nb (Fig. 4E).
Bentonite MREE and HREE compositions retain evidence for changes in mineralogical stability associated with variations in crustal thickness and associated changes in the pressure/temperature regime. A geochemical proxy for the presence of crustal amphibole is the Dy anomaly (Dy/Dy*; Davidson et al., 2013), and the ratio of MREE/HREE (Gd to Lu) is particularly useful for observing amphibole and garnet behavior. In a plot of Dy/Dy* versus Dy/Yb, HREE depletion decreases with a correlative decrease in Dy/Dy* (Fig. 4H). Restitic garnet, stable at greater depths (O'Hara et al., 1971;Robinson and Wood, 1998) and preferentially accommodating HREEs (Wood et al., 2013;Kay and Kay, 1986), can be linked to Albian and Cenomanian bentonites with greater HREE depletion. The positive correlation between Dy/Dy* and Dy/Yb (Fig. 4H) could be explained by anatexis of lower-crustal amphibolite resulting in a garnet-bearing restite (at 10 kbar; Wolf and Wyllie, 1993). Albian and Cenomanian bentonites indicate that the magmas at this time were generated in the thickest crust, but significant overlap in their MREE/ HREE ratios suggests that these magmas were not generated in a uniformly thick crust. It is during this time that the zone of active magmatism broadened while shifting eastward with the shallowing of the Farallon subduction angle.

Magma Migration in the Bentonite 87 Sr/ 86 Sr Record
The 87 Sr/ 86 Sr ratios of our analyzed bentonites can be used to link them to plutonic centers, considering, too, overall proximity. For example, the Sierra Nevada arc was a productive source of tephra throughout the middle to Late Cretaceous (Armstrong and Ward, 1993;Dickinson, 2004), but the Sr isotopic composition of plutonic rocks there is less radiogenic than the bentonites of this study. Sierra Nevada plutons range as low as 0.7041 in Aptian-aged igneous material as old as 121 Ma and maintain ratios between 0.7052 and 0.7058 in the Turonian (DePaolo, 1981); this resulted from much of the Sierra Nevada being intruded across the oceanic-continental crust boundary. These values are significantly lower than age-equivalent bentonites sampled for this study. The Sierra Nevada batholith is over 1000 km away from the Bighorn Basin, and given the general thickness of the bentonites sampled (up to 3 m), the Sierra Nevada batholith seems unlikely to have been the magmatic source, based on known occurrences of significant tephra dispersion (e.g., Taupo; Wilson et al., 2007) and models of supereruptive activity (e.g., Mastin et al., 2014), pointing toward a likely proximity of <100 km. Other potential sources include the Coastal Mountain complex along the coast of British Columbia, but the isotopic composition there is generally too juvenile to be the source (Girardi et al., 2012). Additionally, air currents at this latitude transported ash in an east-northeastern direction (Elder, 1988), placing the Bighorn Basin too far south and distant for this to be likely.
The Idaho batholith magmatic epicenter, one of the few plutonic systems emplaced almost entirely within Laurentian crust (Gaschnig et al., 2011), is a likely source for our analyzed samples based on the fact that its Sr isotopic chronology mirrors the bentonite Sr isotopic chronologic sequence through the Cretaceous (Figs. 8A and 8B). Assimilation of radiogenic Laurentian basement into partial melts generated from Farallon subduction would result in an isotopic mixture (e.g., Ferrara et al., 1986), and the Idaho batholith reflects such mixed signatures (e.g., Gaschnig et al., 2011). The Western Idaho shear zone (also known as the Salmon

Figure 7. Plot demonstrating the correlation between Eu/ Eu* and Dy/Yb of all bentonite data. Less Eu depletion is correlated with a decrease in heavy rare earth element (HREE) sequestration, indicating a thinner (lower-pressure) crustal source. Enhanced Eu depletion is broadly correlated with an increase in HREE depletion, likely indicating that the Eu anomaly is influenced by both crustal thickness and subsurface homogenization.
River suture zone), which marks a 10-km-wide suture between accreted, exotic terranes to the west (e.g., the Wallowa in Oregon) and the Precambrian craton to the east (Kuntz and Snee, 2007), exhibits a sharp isotopic gradient that changes from values less than 0.706 in the west to values greater than 0.706 in the east. From west to east across the Salmon River suture zone, the intermediate to granitic plutons of the Hazard Creek Complex have juvenile signatures less than 0.705, the Little Goose Creek Complex ranges from juvenile to an evolved composition, and the Payette River Complex generally has an 87 Sr/ 86 Sr composition equal to or higher than 0.708 (Manduca et al., 1992). Similar isotopic gradients can be found in the Sierra Nevada system (Chen and Moore, 1982) and many other cratonic margin subduction zones, such as the Andes (Folguera and Ramos, 2011;Ramos et al., 2002). The cratonward migration of the volcanic front as subduction is established is common within continental arc systems. The bentonites show a significant increase from low to high 87 Sr/ 86 Sr through the Aptian into the Albian, mirroring values in plutonic rocks in western Idaho during the eastward migration of pluton emplacement. The late Albian and Cenomanian bentonite Sr compositions mark a stabilization of 87 Sr/ 86 Sr, averaging 0.70803 with only minor deviation for nearly 7 m.y. This time of stabilization coincides with the main phase of Idaho batholith emplacement, classified as part of the "early metaluminous phase" of Idaho activity, and most likely was the result of a well-established mixing-assimilation-storagehomogenization (MASH) zone that reached Sr isotopic equilibrium between basement hostrock assimilation and mantle input. The range of Sr ratios of the bentonites is similar to Idaho plutonic data of this age, though many of the early metaluminous granitoids have been lost (Gaschnig et al., 2017).
Major magmatic activity in Idaho east of the Western Idaho shear zone began in the Cenomanian, continued during the Turonian, and was sustained into the Late Cretaceous (e.g., the Atlanta peraluminous suite), and by 80 Ma, magmatism extended eastward to include a narrow belt of activity in western Montana. Thus, during the Late Cretaceous, magmatism extended across a broad zone (as opposed to the emplacement of plutons migrating eastward as a narrow band). Volcanism in western Montana began with the pre-batholith Elkhorn Volcanic Complex, a precursory system that signaled the emplacement of the Boulder batholith and associated plutons across western Montana (Robinson et al., 1968). The Elkhorn volcanics and subsequent early stages of the Boulder batholith have average 87 Sr/ 86 Sr compositions lower than those of the Idaho batholith (du Bray et al., 2012), and they are more mafic in composition (Tilling, 1973). The isotopic composition of the Atlanta lobe of the Idaho batholith, a later phase of plutonism, ranges from 0.708 to over 0.71; this is in direct contrast to the Campanian-aged Montana-associated plutons, which range from under 0.707 to just over 0.709 at their highest.
Campanian-aged bentonites from the Pierre and Cody Shales and the Mesaverde Formation are also characterized by a wide diversity in their 87 Sr/ 86 Sr compositions, ranging from 0.706 to over 0.71, with distinct clusters of compositional similarity (Fig. 6A). Pierre bentonites are the least radiogenic of all samples (except for the oldest Aptian bentonite), with one Cody bentonite plotting closely. The Cody bentonites show a wide isotopic range, while the Mesaverde bentonites plot tightly near the median Cody bentonites. The Sr data show there may have been at least three source regions of Campanian-aged ash: a less radiogenic source, similar to the Boulder batholith; a source that had a Sr composition around 0.708, similar to the Idaho batholith; and one that was much more radiogenic, exceeding 0.71. This radiogenic source appears to have continued into the Maastrichtian-aged Meeteetse bentonites, where the radiogenic composition drastically narrows from the wide Campanian range. The Meeteetse samples contain the highest radiogenic Sr values of any bentonite, closely matching the composition of the Pioneer batholith in Montana (Zen et al., 1975;du Bray et al., 2012).
Based on their Sr ratios, our analyzed bentonites appear to record two phases of volcanism. The first phase followed the migration of pluton emplacement across the Western Idaho shear zone into central Idaho, where the main phase of Idaho batholith magmatism began during the Cenomanian and was marked by relative constancy of 87 Sr/ 86 Sr at 0.708. The second followed the diversification of Sr composition at the start of the Campanian, as magmatism extended into the area of the Elkhorn volcanics and Boulder batholith in Montana, but most likely continued to produce some volcanism in Idaho. By the Maastrichtian, the bentonites show an enhanced crustal influence with a mean 87 Sr/ 86 Sr value close to 0.71, a compositional range similar to the southwestern Montana Pioneer batholith (Arth et al., 1986). This two-phase sequence follows the geochemical trends observed by  using linear discriminant analysis; they observed two distinct periods of differentiation based on statistical groupings of elemental signatures.

ε Nd Compositions and Model Ages
The distribution of analyzed ε Nd compositions is limited in scope (n = 26) compared to 87 Sr/ 86 Sr measurements (n = 87) but reveals a similarly distinctive isotopic transition from the older (ca. 100 Ma) to younger (ca. 80 Ma) bentonites. Samples show an increase in mantle influence and/ or a decrease in crustal assimilation from Cenomanian to Campanian time. Albian-aged benton-ites have a restricted Nd composition (−11 to −9) that diversifies by the Cenomanian, raising the average ε Nd value to more positive values (−11 to −5). Melting and assimilation of crust likely brought the ε Nd composition of magmas closer to Precambrian basement signatures, explaining the low ε Nd composition of Albian bentonites (near −10) and the coeval Idaho batholith ε Nd (Unruh et al., 2008). As mantle input persisted, the average value slightly increased through the Cenomanian. Strontium compositions support this, but Sr isotopic ratios remain more homogeneous than ε Nd . Through Albian to Cenomanian time, it appears that a mantle component was outpacing the assimilation of crustal ε Nd while maintaining equilibrium with crustal Sr. This discrepancy between Nd and Sr evolution was likely a combination of mantle input outpacing crustal assimilation and the higher Sr/Nd ratio of the continental lithospheric mantle (∼15; Mc-Donough, 1990), meaning that the Nd compositions of magmas were more easily influenced by a mantle contribution. Campanian bentonite shows a diversification of Nd isotopic values, ranging from the low ε Nd values of the Cenomanian to near chondritic reservoir (ε Nd = 0). This range supports the idea that the bentoniteproducing zone saw a significant broadening that extended from Idaho to Montana during the Late Cretaceous. The more positive radiogenic Nd compositions coincide with the initiation of the Elkhorn Volcanic Complex, with low 87 Sr/ 86 Sr, high ε Nd , and elevated concentrations of compatible elements, which resulted from an enhanced mantle contribution (du Bray et al., 2017) and the correlative Nd increase.
Depleted mantle model ages support this Late Cretaceous migratory history and help to isolate the probable source localities by matching basement terranes (Fig. 9A) with DMM Nd ages (Fig. 9B). The Idaho batholith and the Montana plutonic system intruded through two primary, tectonically related terranes: the Selway terrane, a poorly exposed, Proterozoic-aged crustal block that covers much of north-central Idaho (Foster et al., 2006), and the Great Falls tectonic zone, a northeast-trending Proterozoic suture zone of the northern Medicine Hat and southern Wyoming blocks, coeval with the Selway terrane (Foster et al., 2006(Foster et al., , 2012Gifford et al., 2018). The Great Falls tectonic zone and parts of the Selway terrane were later subjected to rifting in B A the Mesoproterozoic, resulting in the formation of the Helena Embayment and deposition of the Belt Supergroup sediments (Foster et al., 2012). The DMM ages of the Albian bentonites are narrowly restricted to the late Paleoproterozoic, and though this Paleoproterozoic signature could have been sourced from the Great Falls tectonic zone, the absence of large DMM age variability that would arise from a shift in provenance suggests the Selway terrane as the probable source. DMM ages diversify in the Cenomanian, ranging from Paleoproterozoic to Mesoproterozoic. The Campanian shows a similar range but extends further into the late Mesoproterozoic. The only source terrane comparable in age to these late Mesoproterozoic signatures is the rift-related material found within the Great Falls tectonic zone, suggesting that this complex zone of rifting and extensive NW-trending basement faults (Gu et al., 2018) channeled plutonic emplacement here as opposed to further south or north within Archean basement terranes. This likely means that some magmas were intruding near the Great Falls tectonic zone as far back as the Cenomanian, and that the eastward migration began broadening early on in the history of the Idaho batholith.

Two Cycles of Idaho Magmatism
A suite of Aptian-aged tephra deposits found in the Cloverly Formation are the oldest ash beds preserved in the Bighorn Basin, and they are part of the renewed magmatism that followed volcanic and tectonic quiescence that began in the Jurassic. This first cycle (Fig. 10A) of renewed magmatism recorded by these bentonites began ∼10 m.y. after the start of the mid-Cretaceous magmatic culmination of the Cordilleran system. Aptian bentonite 87 Sr/ 86 Sr compositions follow a significant increase, mirroring plutonic Sr signatures as magmatism migrated from exotic terranes in the west eastward into Precambrian Laurentian crust in western Idaho. In the bentonite data, there is a small Eu anomaly, indicating differentiation generally had not progressed past rhyodacite, and the Nb/Ta ratio is high. As plutonic emplacement advanced across the complex suture zone, magmas most likely had a direct route toward the upper lithosphere, decreasing crustal residence time, and depressing differentiation. Subsequent crustal shortening during this magmatic flux (Smith et al., 1988) resulted in rapid eastward migration across the crustal boundary and into central Idaho.
The second stage of this first cycle (Fig. 10B) began with bentonites of Albian age (collected from the Thermopolis Shale) and continued through bentonites in the Cenomanian-aged Mowry and Frontier Formations. The 87 Sr/ 86 Sr compositions in the Albian are wide-ranging, but they stabilize by the Cenomanian, maintaining a signature near 0.708 through the rest of the stage. There is a concurrent increase in Eu depletion coupled with enhanced differentiation. Plagioclase and titanite fractionation was an important factor in the generation of these magmas and their subsequent evolution. Low Nb/Ta ratios match the enhanced differentiation index typical of highly evolved granitic magma (Green, 1995). Homogeneous 87 Sr/ 86 Sr and ε Nd compositions are most likely due to lower-crustal MASH processing (Hildreth and Moorbath, 1988), similar to long-lived plutonic systems in the Sierra Nevada (e.g., Lackey et al., 2008) and Andes (e.g., Hildreth and Moorbath, 1988), and this mechanism was responsible for geochemical and isotopic averaging between mantle-and crustal-derived melts. The age-equivalent rocks in the Idaho batholith correspond to the early metaluminous phase plutons, which were largely lost during later Cretaceous magmatism. The Nd DMM ages of the Cenomanian bentonites coincide with the Selway basement terrane of central Idaho. This Precambrian crust would have been conducive to significant lithospheric underplating, allowing for the well-developed MASH processing buffered by a thickened crust. In the Cenomanian bentonites, their large Eu anomaly and enhanced crustal differentiation were likely a result of significant plagioclase extraction during MASH generation (e.g., Sierra Valle Fértil batholith; Walker et al., 2015) during a long residence time from increased crustal thickness. Additionally, this cycle of bentonites has the most pronounced HREE depletion, demonstrating that the crust was thick enough (>30 km) to stabilize garnet, affirming the connection between an increase in crustal thickness and the signatures associated with progressive differentiation.
Flattening of Farallon subduction occurred during the Cenomanian-a key tectonic event in the Cordilleran that is often associated with the convergence of the Shatsky conjugate (Liu et al., 2010). Slab flattening caused the eastward migration of magmatic activity, pushing plutonic emplacement into western to central Montana. The Campanian marked the beginning of a significant change in bentonite composition and the beginning of a second interval when the centers of plutonism migrated eastward (Fig. 10C), signaled by 87 Sr/ 86 Sr ratios varying widely between more primitive and evolved signatures. The ε Nd compositions followed the same diversification, and 87 Sr/ 86 Sr versus ε Nd plots show a trend toward chondritic uniform reservoir. The initial pulse of magmatism in western Montana began with the Elkhorn Volcanic Complex, and bentonite compositions reverted to a reduced Eu anomaly linked to higher Dy/Yb, less differen-tiation, and a larger Nb/Ta ratio, all indicators of diminished crustal residence. Bentonite chemistry coincides with geochemical observations from the Boulder batholith, which include minimal internal differentiation, limited Eu anomaly, and modest plagioclase fractionation (du Bray et al., 2012). The intrusion of plutons through the Great Falls tectonic zone may have been the reason for this change. In this complex crustal boundary characterized by an ancient suture zone and later extension, MASH zone processing is less likely to have been sustainable, resulting in a more direct pathway for magma to reach the upper crust. This may be the reason that the Campanian Elkhorn volcanics (Tilling, 1973) and Boulder batholith (du Bray et al., 2012) have mafic geochemical signatures.
The last period of volcanism (Fig. 10D) is found entirely within the Meeteetse Formation bentonites of Maastrichtian age. It is the last evidence of volcanism preserved well enough within the Bighorn Basin (with a minimal terrigenous input) to reliably trace the geochemistry back to the source. Their Sr isotopic compositions are the highest of all sequences and closely match those of the Pioneer batholith, part of the Montana plutonic belt located southwest of the Boulder batholith. The emplacement of the Pioneer batholith marginally overlapped the Boulder batholith but extended later, matching the depositional age and isotopic composition of the Meeteetse bentonites (Arth et al., 1986). The trace elements follow the low differentiation progression of Campanian bentonites, but their high 87 Sr/ 86 Sr composition shows significant crustal assimilation of a Precambrian source with high isotopic values. This fits the Vipond Park gneiss with an 87 Sr/ 86 Sr of ∼0.72 (the Vipond Park gneiss is a Proterozoic unit underlying the Pioneer batholith; Arth et al., 1986). The mixing of this radiogenic source with mantle input could have produced the Sr composition of the batholith and the resulting ash. The most voluminous magmas of the Pioneer batholith have an initial Sr ratio of 0.7112, a composition isotopically similar to the Meeteetse bentonites.

CONCLUSIONS
The trace-element and Sr and Nd geochemical signatures preserved within 40 m.y. of Cretaceous-aged bentonite sampled from available Cretaceous strata in the Bighorn Basin and South Dakota record evidence of magmatic processes that appear to follow trends associated with differences in crustal thickness and allow these ash beds to be correlated with plutonic rocks in Idaho and Montana. While glass alteration results in some bentonite source-information loss, immobile elements provide key information, including REE patterns, the Eu anomaly, and Zr/TiO 2 and Nb/Ta values. The role of garnet is preserved in bentonite LREE/HREE ratios, and their Dy anomalies show the role of amphibole in crustal processing. Bentonite 87 Sr/ 86 Sr and ε Nd compositions are also valuable in determining the reconstructed age of the assimilated basement.
Bentonite geochemistry can be divided into two distinct periods with different signatures, one beginning in the Aptian and ending in the Turonian, and the other beginning in the Campanian and continuing through the end of the Cretaceous. Bentonite geochemistry can be correlated to plutonic rocks that transgressed across the Cordillera-Laurentian boundary into Idaho, stalled in central Idaho until the Campanian, and then pushed further eastward into Montana with the onset of Farallon slab shallowing. As the active zone of plutonic emplacement migrated eastward from accreted terranes into Laurentian crust, much of the ash preserved in the Cloverly Formation and early Thermopolis Shale bentonites of Aptian age was produced. Once the plutons reached the thick Precambrian crust of Idaho, magmatism was likely buffered by MASH processing, and the stable composition and position of plutons in central Idaho until ca. 85 Ma are recorded in the bentonite deposits of the late Thermopolis Shale, Mowry Shale, and Frontier Formation, all characterized by low Nb/ Ta, enhanced Eu anomalies, and stable 87 Sr/ 86 Sr of 0.708. By 80 Ma, the influence of slab shallowing is seen by the wide broadening of the zone of plutonism, spreading from central Idaho to western Montana, and concurrent lowering of Dy/Dy* and Dy/Yb, diminished Eu anomaly, and increased range of Sr and Nd isotopic values, including lower 87 Sr/ 86 Sr ratios and ε Nd values approaching MORB. Magmatism continued to evolve through a system of upper-crustal plutons, producing isolated, small plutonic bodies in western Montana that continued to produce ash through the Maastrichtian, resulting in the Meeteetse Formation bentonites, which have the most radiogenic Sr, the smallest Eu anomalies, and low Dy/Yb values.
The mobility of key elemental signatures through devitrification, including alkali metals and Y, limits the complete one-to-one comparison of ash to bentonite; however, when used to reconstruct the large-scale magmatic history, Cretaceous bentonites in the Bighorn Basin and South Dakota record distinct changes in geochemistry that reflect temporal patterns in crustal inheritance, magmatic differentiation, residence time, and associated lithospheric thickness. Thus, bentonite data can be a unique proxy for tracking magmatic processes, including the timing of magmatic migration in regions where plutons are preserved poorly, and future efforts could enhance this ability with well-constrained phenocrystic data (i.e., zircon).

ACKNOWLEDGMENTS
This work was completed as the culmination of Hannon's Ph.D. dissertation with the University of Cincinnati's Department of Geology. We would like to thank Richard Brown of Wyo-ben, Inc., for allowing us access to sampling localities, expertise, and hospitality over two field seasons. We would also like to thank Kristen Hannon for assistance in the field. This work was partially funded by the Clay Minerals Society, the Geological Society of America, and the University of Cincinnati Graduate Student Governance. We would like to thank two anonymous reviewers of an earlier version of this manuscript, whose comments and suggestions significantly improved this paper. The comments and questions of Associate Editor Jean Bédard also sharpened our thinking and helped to improve the manuscript, and for this, we are grateful.