Evidence for Cyclical Fractional Crystallization, Recharge, and Assimilation in Basalts of the Kimama Drill Core, Central Snake River Plain, Idaho: 5.5-Million-Years of Petrogenesis in a Mid-crustal Sill Complex

Basalts erupted in the Snake River Plain of central Idaho and sampled in the Kimama drill core link eruptive processes to the construction of maﬁc intrusions over 5.5 Ma. Cyclic variations in basalt composition reveal temporal chemical heterogeneity related to fractional crystallization and the assimilation of previously-intruded maﬁc sills. A range of compositional types are identiﬁed within 1,912 m of continuous drill core: Snake River olivine tholeiite (SROT), low K SROT, high assimilation. These evolved lavas may be sourced from the Craters of the Moon/Great Rift system to the northeast. The Kimama drill core is the longest record of geochemical variation in the central Snake River Plain and reinforces the concept of magma processing in a layered complex.


INTRODUCTION
The relationship between volcanic rocks erupted on Earth's surface and their putative plutonic roots is a long-standing geologic problem (e.g., Bachmann et al., 2007;Annen et al., 2015;Cashman et al., 2017). The petrologic and chemical evolution of volcanic provinces and the formation of continental crustal plutonic complexes are typically studied in isolation from one another, yet each of these topics has implications for the other related to magmatic flux, intrusion geometries, fractionation mechanisms, and the effects of partial melting and assimilation (Annen et al., 2015). Previous efforts to understand linkages between pluton construction and eruptive processes have focused largely on silicic systems: the assembly of granitic batholiths and rhyolite ignimbrite eruptions (e.g., Lipman, 1984;Metcalf, 2004;Annen et al., 2006). Although large mafic layered intrusions have been studied extensively (e.g., Irvine, 1970) and the effect of basalt intrusions on crustal heat flux has been modeled (Annen and Sparks, 2002), there are few studies that infer mafic plutonic evolution from basalt compositions in long-lived mafic volcanic provinces (e.g., Shervais et al., 2006). In this paper we use the petrologic and geochemical evolution of basalt recovered in drill core to infer crustal magmatic processes over 5.5 Ma of mafic volcanism in continental crust.
The Snake River Plain (SRP) of central Idaho is a manifestation of Yellowstone-Snake River Plain (YSRP) hotspot volcanism and contains intriguing evidence for mantle hotspot impingement on continental crust (Pierce and Morgan, 1992;Hanan et al., 2008;Shervais and Hanan, 2008;Smith et al., 2009;Jean et al., 2014). The region contains a record of continuous bimodal volcanism extending over 17 Ma Coble and Mahood, 2012;Henry et al., 2017) to the present, and documents the migration of timetransgressive rhyolitic volcanism from the Bruneau-Jarbidge volcanic center (circa 12 Ma) to its present location beneath the Yellowstone Plateau (Pierce and Morgan, 1992;Anders et al., 2009). Interaction between the mantle hotspot and overlying continental lithosphere has resulted in large rhyolite calderaforming eruptions, followed by eruption of smaller basaltic shield volcanoes (Bonnichsen et al., 2008;Christiansen and McCurry, 2008;McCurry and Rodgers, 2009). Post-caldera basaltic flows form a veneer over the rhyolite ash flows and lavas and mask the complete volcanic record. Understanding the origin and evolution of the post-caldera basalts is a challenge because the lack of later uplift and erosion means that younger lava flows conceal older basalts. High-recovery continental scientific drilling is an effective method to address questions of geochemical and volcanic evolution through time, especially in relatively young volcanic provinces that lack exposure.
Project Hotspot drilled three deep boreholes in the Snake River Plain in order to provide a more complete understanding of the volcanic history of the SRP (Shervais et al., , 2012. The 1,912 m deep Kimama borehole, located in the Axial Volcanic Zone of the Snake River Plain (Figure 1), recovered over 1,900 m of continuous core, including 1,803 m of basalt and 110 m of intercalated sediment (including sidetrack core; Potter, 2014). We present here a detailed petrologic and geochemical investigation of basalts from the Kimama core to determine the nature and extent of chemical changes through time at a fixed location, where physical and chemical characteristics of the crust and mantle lithosphere are relatively set. These data document long-term (5.5 Ma) petrogenetic processes common to intracontinental basalts associated with mantle hotspots and large igneous provinces.  Shervais and Hanan (2008); location of Great Rift and Craters of the Moon from Hughes et al. (2002a).
Other deep boreholes on the SRP include the Sugar City borehole (0.7 km; Embree et al., 1978;Jean et al., in revision), the Idaho National Laboratory (INL) hole WO-2 (1.52 km; Shervais et al., 2006), and the Wendell-RASA borehole (0.3 km; Jean et al., 2013) (Figure 1). Geochemical and isotope analyses conducted by previous workers indicate that significant fractionation occurred at lower or mid-crustal depths in layered mafic magma reservoirs, and that interaction with continental lithosphere is the primary external influence on SRP basalt chemistry. The most important mechanism of magma differentiation is thought to be fractional crystallization accompanied by assimilation of genetically-related, previously-intruded mafic sills in mid-crustal magma chambers (Shervais et al., 2006;Jean et al., 2013). The range of primary magma compositions suggests the involvement of multiple, small magma batches (Leeman, 1982;Vetter and Shervais, 1992;Geist et al., 2002;Hughes et al., 2002a,b;Shervais et al., 2006;Jean et al., 2013). Shifts in the composition of primary magma sources through time and space have been documented by Vetter and Shervais (1992), Shervais et al. (2006), and Shervais and Vetter (2009). As the longest and most complete record of basalt volcanism and magma differentiation in the SRP, the Kimama drill core ties mafic volcanism to a transcrustal network of sills and sill complexes where magmas are processed and cycled to the surface.

Regional Setting
The central SRP is loosely defined as the portion of the SRP between the Owyhee Plateau, a highland in southwest Idaho, and the Great Rift, a north-northwest-trending fissure system that extends ∼50 km southward from Craters of the Moon lava flow field to the Wapi shield volcano (Kuntz et al., 1982(Kuntz et al., , 1992. Major rhyolitic volcanic fields of the central SRP include the 11-6 Ma Twin Falls eruptive center (Knott et al., 2016), and the 10.4-6.6 Ma Picabo eruptive center, for which the boundaries are poorly defined (Figure 1) (Bonnichsen et al., 2008;Drew et al., 2013). These rhyolitic centers are buried by younger basaltic lava flows. Beyond the SRP, to the north and south of the Kimama drill site, is the Miocene to recent Basin and Range Province.
During the late Pliocene through Pleistocene, the SRP was the locus of densely-spaced mafic volcanic centers along the axis of the YSRP hotspot track that erupted thick packages of basalt flows and formed the Axial Volcanic High (Hackett et al., 2004;Shervais et al., 2005). Mafic volcanism on the SRP began within 1-2 Ma of the cessation of YSRP silicic volcanism, and is primarily expressed by the eruption of monogenetic olivine tholeiite basalts (Hughes et al., 2002a). Although timeprogressive basalt eruptions are evident at the inception of postrhyolitic volcanism, Pliocene-Holocene basalt lavas on the SRP are distributed throughout the volcanic province and do not follow the progression of silicic volcanism to the NE (Hughes et al., 1999;Bonnichsen and Godchaux, 2002). Most SRP basalt volcanism is manifested as monogenetic, single-pulse lava fields, erupted as low volume (3-5 km 3 ) flows from shields during short-duration (tens to hundreds of years) eruptions (Kuntz et al., 1992).
The unique style of SRP volcanism has been recognized as a product of small, mid-crustal magma chambers feeding eruptions from coalesced low-relief shield volcanoes, over relatively short durations. Volcanism is described as "plainsstyle" by Greeley (1982). "Plains-style" volcanism is similar to modern flow processes on the island of Hawaii and transitional between quiescent Hawaiian-style and the large volume flood-style volcanism typified by the Columbia River Basalt Group (Greeley, 1982). Vent constructs for SRP volcanoes are typically low due to the efficient transport of lowviscosity lava away from the vent in lava tubes. Eruptive centers often blend in with the surrounding topography because relatively little lava accumulates near the vent (Self et al., 1998). Between periods of volcanism, lava flows in low areas were mantled by loess and by lacustrine and fluvial sediments.

Snake River Plain Basalt Stratigraphy
The high effusion rate of basalts on the SRP favors emplacement via inflated pahoehoe flows, whereby lava travels through the molten core of the flow and breaks out of the flow front as flow lobes (Walker, 1991;Chitwood, 1994;Hon et al., 1994;Self et al., 1998). The transport of pahoehoe basalt lava over large distances is facilitated via lava tubes, tube complexes, or surface channelization by inflation processes (Self et al., 1998;Hughes et al., 1999). On the Snake River Plain, lava flows may extend as much as 50-60 km from the vent of origin (Shervais et al., 2005). Self et al. (1998) recognized that the inflation mechanism of pahoehoe produces lava flows that display similar geometries at variable spatial scales. In order of increasing scale, basalt lavas are emplaced as flow lobes (flow units), flows, and flow groups (Walker, 1991;Chitwood, 1994;Self et al., 1998).
Over the duration of a "plains-style" eruption, sheet flows and flow units from a single magma reservoir interlayer to form a complex aggregate of flows termed flow group (Hughes et al., 2002a,b;Welhan et al., 2002). Lava flow groups have thicknesses of <1-106 m (flow groups 12 and 73, respectively). They are the two-dimensional equivalents of compound lava flow fields, such as the Wapi and Hell's Half Acre lava fields on the SRP, which have areal dimensions of kilometers to tens of kilometers. Each flow group or lava flow field comprises all of the volcanic activity from a monogenetic volcanic center over time scales of tens to hundreds of years (Greeley, 1982;Kuntz et al., 1992;Welhan et al., 2002;Champion, pers. commun. 2014). In core, flow groups are distinguished from one another by lithostratigraphic relationships, paleomagnetic polarity reversals or inclination changes, sediment interbeds, and by significant geochemical variation.
Monogenetic eruptive centers on the SRP are most likely fed by individual magma sources, as demonstrated by chemical variation between and within flow groups (Hughes et al., 2002b). However, in order to produce separate magma batches for every monogenetic center on the SRP, a stratified source region in which melts were produced over a range of depths and degrees of melting would be required (Hughes et al., 1997).

METHODS
We selected 265 whole rock samples from the Kimama core for analysis, representing 74 of 78 identified basalt flow groups. Samples were preferentially selected from the massive interiors of lava flows to avoid oxidation and alteration. Major elements were analyzed by fused bead X-ray fluorescence spectrometry (XRF), and trace elements were analyzed by inductively coupled plasma mass spectrometry (ICP-MS).
One-inch diameter mini cores were extracted from the whole round cores and a portion of each plug was broken in two to three fragments using a rock hammer. Samples were crushed using a Gyral Grinder shatterbox with a tungsten carbide vessel, and then ground again with an agate mortar and pestle. Samples for major element analysis were ignited at 800 • C for 24 h, after which 1 g of ignited sample was mixed with 5 g of a Claisse Liborate flux and 6 drops of LiI (added as a releasing agent). Sample mixes were melted in Pt-Au metal crucibles at 1200 • C in a muffle furnace, then poured into a heated Pt-Au metal disk mold and quenched into glass. Major element (SiO 2 , TiO 2 , Al 2 O 3 , MnO, Fe 2 O 3 , MgO, CaO, Na 2 O, K 2 O, P 2 O 5 Cr 2 O 3 ) were analyzed with a Philips 2400 X-Ray fluorescence (XRF) spectrometer at Utah State University or with a Rigaku ZSX Primus II XRF spectrometer at Brigham Young University.
Trace-element concentrations were measured at Centenary College (Shreveport, LA) using a PerkinElmer 600 inductively coupled plasma-mass spectrometer (ICP-MS). Approximately 60 mg of each sample was dissolved in 2 mL HF and 3 mL HNO 3 for 3 h, with watch glasses preventing evaporation. Watch glasses were then removed and samples dried, after which another 3 mL of HNO 3 was added to the samples. The sample solution was left at 50 • C overnight, and then dried at 90 • C. A further addition of 3 mL of HNO 3 to the sample preceded immediate drying. Finally, the sample was brought into solution with 2-3 mL of 50% HNO 3 , and brought to a total volume of 50 mL with 5% HNO 3 . This procedure is modified from Jenner et al. (1990) and Neal (2001). Five milliliters of 10 ppb In, Rh, and Ru were added to the sample solution as internal standards to calibrate measured concentrations. Plasma lab software was used to map out sampling order and record measurements over the duration of the experiment. Major and trace element analyses made at Utah State University (XRF), Brigham Young University (XRF), and Centenary College (ICP-MS) are distinguished by symbols in Table 1 in Supplementary Material. To compare Kimama major and trace element analyses with compositions from regional SRP basalt lavas, we used published values from Kuntz et al. (1992) and Hughes et al. (2002b).

Flow Groups and Flows in the Kimama Core
The Kimama area includes late Neogene to Quaternary basalts that were erupted from low-relief shield volcanoes, most of which are normal magnetic polarity lavas erupted during the (>781 ka) C1n Brunhes Normal Polarity Chron (Shervais et al., 2005). 21 paleomagnetic Chrons and subchrons were identified in the Kimama core (Potter, 2014). Basaltic lava flows overlie rhyolite in the Twin Falls and Picabo volcanic centers (Shervais et al., 2005(Shervais et al., , 2012Knott et al., 2016). Rhyolite was not encountered in the Kimama drill core Potter, 2014).
Potter (2014) identified 432 basalt flow units comprising 183 lava flows. In core, flows and flow groups are distinguished by lithostratigraphic relationships, paleomagnetic polarity reversals or inclination changes, sediment interbeds, and by significant geochemical variation. The lithostratigraphy, paleomagnetic characteristics, the paleomagnetic and Ar-Ar age model, and the accumulation rate of the Kimama drill core basalts are discussed in Potter (2014).
Using this stratigraphic framework, paleo-secular variations in magnetic stratigraphy, and geochemical indicators (e.g., FeO * = total iron as ferrous, TiO 2 , K 2 O, La/Lu), we recognize 78 distinct flow groups in the Kimama core, ranging from 0.274 to 105.6 m thick. (Figure 2, Table 2 in Supplementary Material). Flow groups 41 and 42, identified by paleomagnetic indicators, were too thin and rubbly to sample. Flow groups 57 and 75, identified by distinct stratigraphic and paleomagnetic indicators, were not sampled due to alteration. Based upon paleomagnetic and Ar-Ar data, as well as stratigraphic relationships, we document that this sequence represents ∼5.5 million years of volcanism (from 6 Ma to 475 ka) (Potter, 2014).
Loss on ignition (LOI) is largely negative above 1055 m depth, with an average of −0.56 (range 0.34 to −2.24) wt% as a result of weight gain as FeO is oxidized to Fe 2 O 3 during heating. Below Low K SROT, high Fe-Ti, high K-Fe, and evolved flow groups are represented by green, blue, and red layers (respectively). A green gradient on gray layers demonstrates the location of low K alteration within SROT flow groups. Sediment interbeds of loess (orange) in the upper section of the drill core, and sandstone and siltstone (orange and yellow) are present in the lower depths of the drill core. Hyaloclastite (purple) is present at the base of the drill core.
Frontiers in Earth Science | www.frontiersin.org 1,055 m depth, loss on ignition values are largely positive, with an average of 1.43 (range 9.41 to −0.83) wt%. The highest LOI (8.5-9.4 wt%) are found below 1,700 m depth, indicating extensive alteration by hydrothermal fluids.
These data reveal that the Kimama basalts define three geochemically-distinct groups: (1) Snake River Olivine Tholeiite (SROT), which includes most of the basalts from surface outcrops on the SRP, (2) tholeiitic basalts rich in iron and titanium ("Fe-Ti basalts"), and (3) a suite of rare high K 2 O-FeO * -rich basalts similar to mafic lavas from the Craters of the Moon. In addition, there is a fourth group of low-K basalts that includes members of Group 1 basalts and is characterized by low K 2 O (typically <0.25 wt%) and K 2 O/TiO2 ratios of ≤0.11 wt% (Figure 2).
Group 1 SROT dominate the sample set, with >90% of all analyzed samples from throughout the length of the drill core. They are similar in composition to ocean island basalts, and characterized by phenocrysts or glomerocrysts of olivine and plagioclase. Group 2 Fe-Ti basalts are characterized by high TiO 2 (3.6-4.5 wt%) and FeO * (15-18 wt%). They occur at five depth ranges: 343-359, 459-465, 522-548, 561-566, and 1,117-1,151 m depth. The Fe-Ti basalts fall on extensions of the SROT trends on element variation diagrams, and are inferred to be related to SROT by low-pressure fractionation. The Group 3 high K-Fe lavas are characterized by K 2 O >0.65 wt% and 17.5-18.5 wt% FeO * , with distinctively high K 2 O/TiO 2 ratios >0.5, and low MgO contents. They are found at only two depth ranges: near 318 m (flow group 12), 1,048 m (Flow group 49), and 1,077 m (flow group 51) (Figures 2, 5, Table 2 in Supplementary Material).
Group 4 includes 36 samples; all but one of these occur below 1,043 m depth. The low-K basalts appear to be normal SROT basalts that have lost K 2 O in response to hydrothermal alteration. This is supported by the correlation of low-K basalts with high K 2 O/TiO 2 ratios and high LOI (Figure 6, Table  1 in Supplementary Material). Alteration within the Kimama drill core is manifest as vesicle fillings at around 700 m depth, forming amygdules filled with calcite and quartz. By 950 m depth, amygdules filled with smectites and/or zeolites become common. Below 1,300 m, alteration becomes more pervasive, with iddingsitization of olivine phenocrysts and the alteration of interstitial glass to clay minerals (Sant, 2013).
Finally, there is one sample that cannot be easily fit into any of these groups, although it seems most comfortable in the high-K-Fe group: KA1B3901 (flow group 59). This sample is characterized by relatively high SiO 2 (53 wt%), high K 2 O (1.57 wt%), and moderate MgO (6.21 wt%) ( Table 1 in Supplementary Material).

Trace Elements
Trace element variations in the Kimama basalts generally mimic trends observed in correlative major elements. Thus, Zr, Nb, Y, and the REE correlate positively with FeO * , TiO 2 and P 2 O 5 , and Ba and Rb correlate with K 2 O, whereas compatible elements such as Ni and Cr correlate positively with MgO.
Kimama basalts are modestly enriched in the light rare earth elements (LREE) relative to the heavy rare earth elements (HREE), with chondrite-normalized La/Lu = 3-8 and smooth patterns that lack Eu anomalies (or have very modest negative anomalies); they are generally similar to those observed in surface basalts of the eastern and central SRP ( Figure 7A). This includes all SROT, Fe-Ti, and low-K SROT. Evolved high K-Fe basalts (Group 3) have similar slopes on chondrite-normalized plots but have higher concentrations, similar to or slightly higher than evolved basalts and trachybasalts from Craters of the Moon ( Figure 7A).
Multi-element (spider) diagrams document further similarities to published data for surface basalts and evolved Craters of the Moon-type lavas. Notable features include strong enrichments in Ba (similar to other basalts in the western USA), and depletions in Th, Sr and Y relative to primitive mantle values ( Figure 7B).

Chemical Suites and Stratigraphic Trends
Despite the overall chemical similarity in SRP Neogene basalts, flow groups in the Kimama drill core exhibit several compositional trends (Figure 2). Whole rock MgO, FeO * , TiO 2, P 2 O 5 , and K 2 O (all wt%), plotted against stratigraphic depth isolate three basalt compositional types (Figure 4). Oxides such as Al 2 O 3 and Na 2 O vary without discernible trends, perhaps owing to the crystallization and flotation of plagioclase within the magma, however CaO (not shown) also varies systematically, but Al 2 O 3 and Na 2 O (also not shown) have no discernable trends.
SROT and high K-Fe flow groups are thought to represent unrelated magma batches, whereas the Fe-Ti flow groups probably represent evolved SROT magmas. In the following paragraphs, we discuss whether Fe-Ti, high K-Fe, and SROT basalts are compositionally related by similar petrogenetic processes, and the significance of low-K basalts below 1200 m depth. We also discuss the significance of the anomalously evolved sample KA1B3901 in flow group 59.
SROT flows below ∼1,200 m depth are dominantly primitive, with low concentrations of FeO * , TiO 2, P 2 O 5 , and K 2 O, and relatively high MgO, CaO, Cr, and Ni. The latter are less susceptible to alteration than K and are indicators of recharge FIGURE 4 | Selected major and trace element variation diagrams with compositions for basalts from the Kimama drill core (colored symbols denote compositional type) compared to compositions of eastern Snake River Plain (ESRP) olivine tholeiites (crosses; Hughes et al., 2002b) and Craters of the Moon evolved lavas (circles; Kuntz et al., 1982). All diagrams show variation with respect to MgO wt%: (A), TiO 2 wt%; (B), FeO* wt%; (C), K 2 O wt%; (D), SiO 2 wt%; (E), P 2 O 5 wt%, and (F), Zr ppm.
processes. Many of these SROT flows are also low-K basalts, with abnormally low K 2 O due to hydrothermal alteration. Nonetheless, the low-K SROT flows below 1,200 m depth are also low in FeO * , TiO 2, and P 2 O 5 , with high MgO, CaO, Cr, and Ni. Above 1,200 m depth, most SROT flows display cyclic variations in composition that are correlated with fractional crystallization and magma recharge, as discussed below.
SROT flows also exhibit a generalized increase of incompatible elements up section (e.g., K 2 O and P 2 O 5 ) that may reflect fractionation enrichment and/or the assimilation of previously emplaced ferrodioritic melts. Flow groups deeper in the drill core exhibit more primitive major and trace element concentrations (e.g., Figure 5). This trend may result from the addition of more enriched primitive magmas or progressively higher amounts of assimilation, as discussed below. Concentrations of REE and other incompatible elements also vary with depth, and substantiate recharge and fractionation cycles observed in the major element data (Figure 7).
Three flow groups were classified as high K-Fe lavas based upon anomalously high values of FeO * (≥17.8 wt%), K 2 O (1.77 to 1.84 wt%), and TiO 2 (3.43-3.54 wt%), coupled with low MgO (3.40-3.94 wt%). Flow groups 12, 49, and 51 have compositions that are similar to those erupted within Craters of the Moon during the latest Pleistocene and Holocene. High K-Fe lavas are located at 318, 1,045, and 1,077 m depths, respectively ( Table  2 in Supplementary Material). Sample KA1B3901 (flow group 59) has low FeO * (10.2 wt%) and TiO 2 (1.91 wt%) relative to the high K-Fe lavas, but nevertheless displays a similar evolved geochemical signature, including high SiO 2 (52.9 wt%), high Na 2 O (2.70 wt%), and high K 2 O (1.57 wt%). The members of the high K-Fe and evolved suites are the most chemically distinct flows in the Kimama drill core.
The REE patterns of SROT, low-K SROT, and Fe-Ti suites are broadly similar, implying a genetic connection. High K-Fe flow groups are 10x more enriched in both LREE and HREE concentrations than other Kimama chemical groups and fall within the Craters of the Moon compositional array (Figure 7).

Hydrothermal Alteration with Depth
Normal SRP lavas contain little to no water and are characterized by negative loss on ignition values, which reflects the conversion of FeO to Fe 2 O 3 during ignition; this effect is especially noticeable because all SRP tholeiites have relatively high FeO/Fe 2 O 3 . Positive ignition loss reflects bound water in excess of that needed to offset the oxidation of ferrous iron. Since these basalts are basically anhydrous when fresh, this bound water represents hydrothermal alteration of groundmass glass to smectite, as documented by Sant (2013). Temperature logs from the Kimama bore hole record a geothermal gradient of 72-75 from about 980-1,440 m depth, with an average heat flow of 132 mW/m 2 (Lachmar et al., 2017). Plots of LOI vs. depth show that hydrothermally altered basalt is dominant below 1050 m depth, although the effects are variable and do not correlate in a simple fashion; only the deepest flows (>1,700 m) have consistent low K 2 O and high LOI (Figure 6).
The effect of hydrothermal alteration on major element concentrations is seen most clearly in the K 2 O/TiO 2 ratios: ratios > 0.11 follow normal igneous trends consistent with other SROT basalts (e.g., Hughes et al., 2002b). Samples with K 2 O/TiO 2 ratios ≤ 0.11 plot below other SROT analyses. During normal differentiation of SROT basalts, K 2 O and TiO 2 behave incompatibly. Since TiO 2 is generally immobile in hydrothermal fluids, the low K 2 O/TiO 2 ratios must represent leaching of K 2 O from the samples during the conversion of groundmass glass and mesostasis to smectite. On a plot of TiO 2 vs. K 2 O, the low-K group clearly lies below the main field of SROT basalts, consistent with this interpretation (Figure 3C). All samples with K 2 O/TiO 2 ratios ≤ 0.11 have positive LOI values, but not all samples with positive LOI have low K 2 O/TiO 2 ratios.
A plot of K 2 O/TiO 2 ratio vs. depth shows that the low-K samples all lie below 1,050 m depth, although many samples from >1,050 m depth plot on a trend consistent with the SROT at less than 1,050 m depth, suggesting that these samples still retain their primary signatures.

Fractionation/Enrichment and Recharge Cycles
The progressive enrichment of incompatible elements in stratigraphically upsection flow groups is commonly interpreted to represent ongoing eruptions from a fractionating magma chamber. Likewise, progressive depletion in incompatible elements, and concomitant enrichment of compatible elements upsection in a flow group, is interpreted to represent magma chamber recharge with primitive or parental melt compositions (Shervais et al., 2006). Upward fractionation cycles are indicated by a decrease in MgO (Cr, and Ni) and increases in FeO * , TiO 2 , P 2 O 5 , and K 2 O as a result of fractionation of Mg-rich olivine, Ca-Al-rich plagioclase, and Cr-spinel. Recharge cycles are indicated by increasing MgO and Ni, and decreases in FeO * , TiO 2 , K 2 O, and P 2 O 5 (Figure 5).
Stratigraphically defined flow groups based on lithology, paleosecular magnetic variations, and sediment intercalations, as defined by Potter (2014), appear to correspond to magma fractionation or recharge cycles (Figure 5). However, flow groups 7, 9, 18, 20, 23, 48, 53, 54, 57, 58, 66, 68, 70, 75, and 78 show relatively constant compositions with depth, suggesting the emptying of a magma reservoir over relatively short time spans, with little or no coeval fractional crystallization. This is consistent with relatively high eruption rates that could empty a magma reservoir quickly.
Other flow groups document progressive recharge of their magma chambers coeval with ongoing eruptions, such that lavas become progressively more primitive upsection, trending to low-K, high MgO compositions. This can be observed in six flow groups: 4, 16, 26, 27, 30, and 39. Flow groups 14, 38, 50, 52, and 60 also display brief intervals of recharge superimposed on the dominant fractionation trends.
The three high K-Fe flow groups and the evolved flow group (groups 12, 49, 51, and 59) are compositionally constant, but too thin to represent prolonged time intervals. Although all four sit atop underlying upward fractionation cycles, the large offset from REE patterns between SROT and high K-Fe or evolved flow groups suggests either a distinct parent magma or the assimilation of a crustal source, likely compositionally similar to Graveyard Point Intrusion gabbros or ferrodiorites (White, 2007). Therefore, connecting high K-Fe and evolved flows to underlying flow groups is problematic.
The occurrence of upward fractionation cycles and recharge cycles is consistent with the proposals of Shervais et al. (2006) and Brueseke et al. (2016) that magmas are processed through mid-crustal intrusions, in which individual sills form complex layers as the cumulate extract of crystal-melt fractionation. These layered mafic intrusions act as reactive filters which process the magma before it erupts. Individual sills may represent a single pulse of primitive magma which differentiates as it feeds a single monogenetic volcano on the surface.

Fractionation-Assimilation
Major processes that may control the evolution of continental basalts include fractional crystallization with no significant assimilation, fractional crystallization accompanied by assimilation of older crustal materials, and fractional crystallization accompanied by assimilation of previously intruded mafic magmas (ACMI-assimilation of consanguineous mafic intrusions; Shervais et al., 2006). These processes can be distinguished using plots of lithophile element ratios vs. a fractionation indicator such as Mg# (Shervais et al., 2006). In Figure 8 we plot Mg# (which decreases with fractional crystallization) against ratios of La/Lu and K 2 O/P 2 O 5 . The ratio of La/Lu varies little during fractional crystallization of olivine and plagioclase, but will increase progressively in response to fractional crystallization-assimilation of either older continental crust or older mafic intrusions. In contrast, K 2 O/P 2 O 5 ratios will increase strongly with assimilation of older crustal materials, but decrease in response to assimilation of mafic intrusions. In Kimama SROT lavas, La/Lu ratios increase progressively as Mg# decreases, whereas K 2 O/P 2 O 5 ratios decrease progressively,  (Kuntz et al., 1982;Hughes et al., 2002b). Trace elements are arranged in order of decreasing compatibility from right to left.
showing that AFC involving older mafic intrusive rocks was the dominant process (Figure 8). The high K-Fe flow group displays the same trends of increasing La/Lu and decreasing K 2 O/P 2 O 5 ratios with decreasing MgO, showing that this group also must form by large degrees of fractional crystallization and assimilation of a gabbroic intrusion. In contrast, sample KA1B3901 (flow group 57) displays an increase in K 2 O/P 2 O 5 ratio with decreasing MgO, and may be a result of fractional crystallization with assimilation of continental crust (Figure 8A). The variation in La/Lu at high mg# in SROT implies a range of parental magma compositions formed either by different degrees of partial melting or by small variations in source composition.

Origin of Highly Evolved Flow Groups
High K-Fe lavas do not appear to be simply related to the parental basalts (Figures 3-5, 7). Variation diagrams of trace elements (Zr/Y-Nb/Y), and MgO variations demonstrate a lack of congruity, and elevated concentrations of REE and other incompatible elements suggest either a different magma source or a unique magmatic differentiation scenario. These lavas were generated in a similar fashion as Craters of the Moon lavas, via extreme fractional crystallization accompanied by assimilation (e.g., Whitaker et al., 2006Whitaker et al., , 2008Putirka et al., 2009). However, the decrease in K 2 O/P 2 O 5 ratios (Figure 8) implies that older felsic crust was not involved forming in these flows (or Craters of the Moon), and that extreme fractional crystallization perhaps accompanied by assimilation of slightly older mafic sills was the dominant process.
Although flow group 59 (sample KA1B3901) has a similar major and trace element composition to Craters of the Moon evolved basalts, high K-Fe flows in the Kimama drill core have distinctly higher FeO, Zr, Sr, Sc, La, and lower Al 2 O 3 , at the same MgO compositions when compared with more "traditional" Craters of the Moon evolved lavas and may represent a different fractionation or assimilation process (Hughes et al., 2002b;Putirka et al., 2009).
The evolved and high K-Fe flow groups in the Kimama drill core may represent the distal portions of early Craters of the Moon lava flows that reached the Kimama site. The fact that these flows are thin is consistent with this idea. High K-Fe lavas range from 0.27 m (flow group 12) to 3.75 m (flow group 49) thick, and the evolved flow group (59) is 6.49 m thick. These thicknesses are well below the average flow group thickness of 23.6 m. However, the Kimama drill site is about 20 km from the exposed terminus of the exposed Craters of the Moon lava flows (Kuntz et al., 2007). This is less than half the maximum observed flow lengths in the SRP (50-60 km; Shervais et al., 2005), so long, thin flows may have reached the Kimama site. On the other hand, age projections from Potter (2014), suggest that high K-Fe flow groups 12, 49, and 51 (317.7-318.9, 1,045.0-1,048.8, and 1,077.1-1,078.5 m depths, and the evolved flow group 59 (1,185.0-1,191.5 m depth) erupted at ∼1.0, ∼3.3, ∼3.4, and ∼3.7 Ma (respectively), well beyond the 480 ka maximum recorded age of Craters of the Moon basalts (Putirka et al., 2009). Thus, either Craters of the Moon-type lavas erupted from other sources farther to the west, or magmatism at Craters of the Moon is much older than documented by surface flows, and the Kimama drill core samples the distal portions of early lava flows.

Melt Source Modeling
Chemical variations in parental basalts may be attributed either to differing degrees of partial melting, different depths of melting, or to different source compositions (e.g., Hughes et al., 2002b;Putirka et al., 2009;Shervais and Vetter, 2009). We have constructed a series of melt models to test these options, using three representative source compositions (primitive mantle, depleted MORB-source mantle, and enriched E-MORB source mantle. Each source region is modeled for a range of melt fractions (from 1 to 20%) using non-modal batch melting in two different depth regimes: shallow spinel-facies melting (<3 GPa or ∼100 km depth) and deeper garnet-facies melting (>3 GPa or ∼100 km depth). We compare the model results to five representative high-MgO (>9.28 wt%) basalts from the Kimama drill core on primitive mantle-normalized multi-element plots (Figure 9). Partial melting of an N-MORB composition source in the spinel lherzolite facies results in depleted LREE concentrations (relative to HREE) in model melts. At extremely low melt fractions (less than 2%), HREE concentrations in the model melts are too high and LREE/HREE ratios too low to match observed primitive basalt compositions (Figure 9A). High field strength elements (HFSE) and K 2 O also show suppressed values in the model melts. Within the garnet lherzolite facies, partial melting of an N-MORB composition source results in model melts that are depleted in LREE relative to HREE. An exception to this observation occurs at extremely low melt fractions (<2%), in which case HREE, K 2 O, and HFSE concentrations in the model melts are too low to fit the data (Figure 9B).
When primitive mantle source compositions are melted in the spinel lherzolite facies (Figure 9C), the model melts generated are too low in LREE and LREE/HREE ratios. Partial melting of a primitive mantle composition in the garnet lherzolite facies produces model melts with LREE concentrations that are too low at large melt fractions and HREE concentrations that are too high at low melt fractions. At all but the lowest melt fractions, the model melts have HFSE concentrations that are too low to match SRP tholeiites (Figure 9D).
The best fit to the trace element patterns of the primitive Kimama basalts was for model melts generated from the partial melting of an E-MORB source in the spinel lherzolite facies (Figure 9E). Coherence is observed for all elements in the 7-15% melting range with the exception of Sr, which is consistently high in all of the model melts. Using the same source composition in a garnet lherzolite mode resulted in fits for the LREE within the ∼10-20% melting range, but poor fits for the HREE (Figure 9F).
The depth of melting inferred from the trace element models (see also Jean et al., 2013, in revision) is also consistent with estimates using phase equilibria and silica geobarometry for Kimama basalts which give melt segregation depths between 80 and 110 km and within the spinel lherzolite facies (Bradshaw, 2012). These relatively low pressures imply that melting must occur within the sublithospheric conduit that has been imaged seismically (Schutt and Dueker, 2008;Stachnik et al., 2008;Obrebski et al., 2010) beneath the eastern Snake River Plain as inferred by Hanan et al. (2008).

Implications for Magma Evolution in Continental Volcanic Settings
Although SRP olivine tholeiites (SROT) are relatively homogenous, chemical variation between separate magmatic events (flow groups) is evident. Previous workers have demonstrated that chemical variation between SRP basalt flows is a result of fractional crystallization, crustal contamination, and partial melting occurring within the mid-crustal sill at pressures of ∼0.6 GPa and temperatures of 1,205 ± 27 • C (Kuntz et al., 1992;Shervais et al., 2006;Christiansen and McCurry, 2008;McCurry and Rodgers, 2009;Miller and Hughes, 2009;Putirka et al., 2009). For Kimama drill core basalts specifically, Bradshaw (2012) estimated pre-eruptive crystallization occurred at 0.2-0.5 GPa (7-17 km) and between 1,155 and 1,255 • C. Putirka et al. (2009) argue for a three-stage process to explain the entire range of SRP lava compositions. For typical SRP olivine tholeiites (SROT), mantle-derived picrites ascend to the middle crust at depths of 20-10 km, where they undergo partial crystallization of olivine ± clinopyroxene. Storage of olivine tholeiite magmas in the middle crust (20-10 km) causes magma compositions to evolve to moderate MgO wt% (10%), at which point positive buoyancy is reached and migration through the middle crust occurs. The lack of clinopyroxene phenocrysts in SRP basalts requires that SRP magmas continued to evolve in shallow subvolcanic reservoirs prior to eruption, forming the low-pressure phenocryst assemblage olivine-plagioclase.
More evolved magmas, such as those associated with Craters of the Moon, require additional processing at depths of 15-0 km, where differentiation and resulting volatile content increases (1-2 wt% H 2 O) causing the final eruption and ascent of magma through the middle and upper crust (e.g., Annen et al., 2006;Putirka et al., 2009).
The data presented here also document the cyclic nature of magma intrusion, fractional crystallization, and magma recharge, consistent with magma evolution in a series of stacked layered intrusions in the lower to middle crust, which has also been documented in other deep drill cores from the SRP (e.g., Shervais et al., 2006;Jean et al., 2013). These data support the SRP magma petrogenetic model of Shervais et al. (2006), in which magmas are shown to evolve through complex pathways of fractional crystallization, assimilation, and mixing at multiple crustal levels before eruption. A particular aspect of this model also supported here is the process referred to as "assimilation of consanguineous mafic intrusions" (ACMI; Shervais et al., 2006). This refers to the process in which the assimilant is not preexisting continental crust (e.g., Shirley, 2013), but rather slightly older intrusions of genetically-related mafic magmas that form the layered mafic sill complex in the middle crust. This process is documented by the fractionation-related increase in LREE/HREE ratios accompanied by decreases in lithophile/high-field strength element ratios such as K 2 O/P 2 O 5 ( Figure 8A).
Assimilation of consanguineaous mafic intrusions is facilitated by two processes. First, newly intruded magma is less dense than previously intruded sills that are now mostly crystallized. This means that new magma will tend to pond above older intrusions, and subsequent intrusions will continue to stack above older ones. Second, the newly intruded magmas will interact with more evolved portions of the older sills, which are typically located in the upper part of the crystallized sill. Overall, the evolved diorites of an older sill will have lower melting points than the sill, thereby facilitating assimilation.
In some cases, if the older sill is not completely crystallized, the freshly recharged magma will preferentially inject into the evolved residual magma (mush) of the older sill, mixing with the evolved residual magma (mush) to form a hybrid magma. This may result in wholesale replacement of the older magma mush, or may proceed gradually, resulting in progressively more primitive magma over time. Our data present evidence for both scenarios: in some cases, fractionation cycles succeed one another with no observed gradual recharge; in others, gradual recharge of the magma chamber takes place, as documented by the recharge cycles. These processes are observed in other deep core suites from the central and eastern SRP (Shervais et al., 2006;Jean et al., 2013).
The occurrence of recharge cycles in particular suggests that magma evolution in the middle crust takes place in relatively large (possibly km-scale?) layered intrusions, not 100m scale sills (Shervais et al., 2006), and that layered mafic intrusions such as the Muskox (Irvine, 1970) may provide more robust analogies for this evolution than smaller sills like the Graveyard Point Intrusion (White, 2007). The continuous increase in concentrations of incompatible elements upsection in the drill core (e.g., La, Figure 5) is more compatible with longer-lived large scale intrusions than with many small, ephemeral intrusions, which are unlikely to result in continuous enrichment of the magma system. These sill-shaped layered intrusions are stacked in the upper middle crust by density contrasts with the underlying lower crust and overlying upper crust (Annen et al., 2015). Each subsequent intrusion will tend to pond above previous intrusions, which are denser than the newly intruded magma, or may tend to localize within the zone of residual melt of a previously intruded sill.
One consequence of this recharge-mixing-assimilation model is that successive magma bodies become progressively more enriched in incompatible elements over time. This is reflected in the progressive increase in K 2 O, P 2 O 5 , and La/Lu upsection (Figure 5). This process is depicted schematically in Figure 10, which shows how mixing of newly recharged magma with residual melt (or assimilated fractional melt) from prior intrusions leads to enrichment of the mixed melt. Over time, successive mixtures become progressively enriched even though the recharge magma remains the same. Similar enrichments may also occur due to in situ crystallization in a large magma chamber (e.g., Langmuir, 1989).
The physical manifestation of our stacked sill model is depicted in Figure 11. Mantle-derived mafic magmas move upwards through the lower crust (not shown) and form large sill-like intrusions in the upper middle crust. The larger silllike plutons (which form during episodes of high magma flux) evolve by crystal fractionation to form layered mafic intrusions (e.g., Irvine, 1970) that are periodically refreshed by new magma. Magma escaping from these layered plutons may erupt directly to the surface, or may be processed through small subvolcanic magma bodies at shallow depths (2-4 km) that are typically associated with a single volcano. As each pluton freezes, new silllike intrusions pond above the older plutons to form a stacked sill complex. This process leads to assimilation of consanguineous mafic intrusions, as later intrusions cannibalize older ones and evolve by AFC processes that generally do not include preexisting crustal rocks (e.g., Shervais et al., 2006).
We suggest that the physical and chemical connections that we observe in the Kimama drill core are common to intra-continental basalts and their crustal plutons, and may characterize mafic magma systems in continental crust, especially when magma flux is relatively high.

CONCLUSIONS
Geochemical data, along with previously-determined drill core stratigraphy and paleomagnetic data (Potter, 2014), show that basalts from the Kimama drill core on the axial high of the Snake River Plain vary between two main compositional arrays throughout 5.5 Ma of continuous eruption. Snake River olivine tholeiites (SROT) and their high Fe-Ti differentiates make up >99% of sampled flows. The remaining flows (<1%) comprise minor volumes of evolved high K-Fe lavas that share a compositional affinity to evolved flows erupted at Craters of the Moon.
The Kimama borehole samples lava flows presumably erupted from multiple shields in the region, but the similarity of mantle melting conditions, transport paths, or depths of evolution has resulted in lavas with generally similar compositions. The compositional variance within SROT and high Fe-Ti lavas observed in the Kimama drill core must reflect differing processes of petrogenesis from magmas of the same source. Low K flows below 1,050 m depth are correlated to increasing LOI, and we invoke secondary alteration of SROT lavas to explain these composition variants. Upward fractionation, or upsection increases in incompatible elements, K 2 O and FeO * suggest fractional crystallization cycles as observed in layered mafic intrusions. Chemical reversals, or upsection increases in compatible elements and MgO concentrations, suggest episodes of progressive recharge of the system by more primitive melt.
We interpret the multiple chemical cycles evident in stratigraphic comparisons of Kimama drill core geochemistry to represent the fractionation of individual magma batches and the progressive recharge of crustal magma reservoirs, with each cycle representing eruptions from a nearby vent or vent system. FIGURE 10 | In large mafic igneous provinces with high magma flux such as the SRP, pre-existing mafic magma intrusions are assimilated by genetically-related newer magma intrusions in a process we call "assimilation of consanguineous mafic intrusions" process (ACMI; Shervais et al., 2006). This conceptual diagram shows the time-progressive effects of ACMI. Ci and Cj are the concentrations of two different incompatible trace elements for which the bulk distribution coefficient (D) of C i is less than that of C j (e.g., D ci <D cj <1). The curved red arrows demonstrate assimilation and fractional crystallization (AFC) pathways that lead to progressively more evolved compositions and where the assimilant is solid intrusive rock of the same magma system. The blue dashed lines indicate recharge by more primitive magma. The mixed magma then experiences another cycle of evolution beginning with AFC. The compositional trend (black line) of Kimama basalts shows the "stair-step" effect of multiple cycles of AFC and recharge leading to more evolved compositions and echoing patterns visible in Figures 5, 6B.
FIGURE 11 | Conceptual diagram of a layered mafic sill complex illustrating the process of "assimilation of consanguineous mafic intrusions" (ACMI). In this model, mantle-derived mafic magmas move upwards through the lower crust and form large sill-like plutons within the middle crust (10-20 km depth) during periods of high magma flux. The large sill-like plutons evolve through fractional crystallization and are periodically refreshed by newer magmas, eventually forming a layered mafic intrusion comprising cumulates of dunite, olivine gabbro, troctolite, and late-stage differentiates including ferrogabbro and ferrodiorite (e.g., Irvine, 1970). New sill-like intrusions pond above cooled sill-like intrusions to form a stacked sill complex. Within the stacked sill complex, new sills cannibalize older sills and undergo ACMI processes without the influence of pre-existing continental crust. Magmas may be transported to the surface through dikes, or they may be stalled and processed in small, subvolcanic intrusions at shallow depths (2-4 km).
Fractional crystallization is accompanied by assimilation, as shown by increases in La/Lu ratios with decreasing MgO, but the concomitant decrease in K 2 O/P 2 O 5 ratios shows that the assimilant is not pre-existing continental crust, but older mafic intrusions related to the continuing SROT volcanism. This is consistent with our FC-AFC models that use the Graveyard Point Intrusion diorite (White, 2007) as an analog for the mid-crustal sill complex. These fractionation and reversal trends provide evidence of a dynamic magma processing system over relatively short to long periods of time.
Our data show the importance of "assimilation of consanguineous mafic intrusions" (Shervais et al., 2006). In this model, large-volume mafic igneous provinces evolve by AFC processes in which the assimilant is dominantly older intrusions related to the same magma suite, not pre-existing continental crust. We suggest that model may be generally applicable to all mafic igneous provinces with high magma flux. This likely require large, km-scale layered intrusions in the middle to lower crust in order to produce the gradual replenishment trends seen throughout the drill core (Figure 11).
The stratigraphically thin (<1.5 m) high K-Fe lavas found in the Kimama drill core are compositionally unrelated to SROT. We propose that high K-Fe flows erupted from Craters of the Moon and flowed southwest into the Kimama area, and that evolved magmas have been generated and episodically erupted from the Craters of the Moon region since at least 3.4 Ma.
Finally, we note that continental scientific drilling of volcanic provinces with >99% drill core recovery is a highly effective way to study both the overall chemical evolution of these provinces and their evolution through time (e.g., Shervais et al., 2014). It is particularly effective for understanding magma evolution in the crust and for inferring magma chamber processes in relatively young volcanic provinces.

AUTHOR CONTRIBUTIONS
JS provided guidance during the preparation of the manuscript, specifically discussions of assimilation processes in the mafic mid-crustal sill complex of the Snake River Plain (central Idaho). JS oversaw the progress and completion on the dissertation that inspired this manuscript as my primary academic advisor at Utah State University (USU). JS also provided lab facilities and materials at USU where the majority of major element analyses on the Kimama drill core samples were conducted. EC provided guidance in the interpretation of major and trace element geochemical trends in Snake River Plain basalts, including identifying the relationship between some major element oxide percentages and loss on ignition (LOI) values. EC also oversaw the preparation and analysis of several Kimama drill core samples via X-ray fluorescence at Brigham Young University. SV provided lab facilities at Centenary College for the preparation of samples and the analysis of trace element concentrations via inductively coupled plasma mass spectrometry (ICPMS). KP logged and sampled the Kimama drill core while a Ph.D student at Utah State University. While a graduate student, KP prepared samples and conducted geochemical analyses, obtaining the data to identify compositional trends and make the larger-scale magma petrogenesis interpretations upon which this manuscript is based. The manuscript is based, in part, on KP Utah State University dissertation research (Potter, 2014), but later geochemical analyses were conducted by KP during her post-doc at Utah State University. The majority of the manuscript was written by KP with minor additions contributed by JS.