Abstract Flood basalts, the largest volcanic events in Earth history, are thought to drive global environmental change because they can emit large volumes of CO 2 and SO 2 over short geologic time scales. Eruption of the Columbia River Basalt Group (CRBG) has been linked to elevated atmospheric CO 2 and global warming during the mid-Miocene climate optimum (MMCO) ~16 million years (Ma) ago. However, a causative relationship between volcanism and warming remains speculative, as the timing and tempo of CRBG eruptions is not well known. We use U-Pb geochronology on zircon-bearing volcanic ash beds intercalated within the basalt stratigraphy to build a high-resolution CRBG eruption record. Our data set shows that more than 95% of the CRBG erupted between 16.7 and 15.9 Ma, twice as fast as previous estimates. By suggesting a recalibration of the geomagnetic polarity time scale, these data indicate that the onset of flood volcanism is nearly contemporaneous with that of the MMCO.

INTRODUCTION The Columbia River Basalt Group (CRBG) is the youngest, smallest, and best-preserved continental flood basalt. It erupted ~210,000 km3 of lava in the Pacific Northwest, United States, between ~17 and 5 million years (Ma) ago. Forty-three distinct stratigraphic members with volume estimates have been defined using regional correlations based on detailed mapping, geochemistry, and paleomagnetic data (1). While many flood basalts have been implicated in mass extinction events (2), the CRBG is not associated with a mass extinction event and its total volume is up to an order of magnitude smaller than other flood basalts. However, single members yielded comparable volumes of lava (thousands of cubic kilometers) to single formations of large igneous provinces implicated in mass extinctions [for example, Deccan Traps (3)] and therefore may have had comparable short-term climate impacts (4). The ~17 to 15 Ma mid-Miocene climate optimum (MMCO) is marked by high-latitude sea surface temperatures 4° to 6°C above background temperature (5) and is associated with vertebrate migrations and increased species originations (6). In paleoclimate records, the MMCO is marked by a benthic δ18O minimum, a benthic δ13C maximum, an ice sheet extent minimum (7), and a variety of pco 2 (partial pressure of atmospheric CO 2 ) proxies, which indicate a possible doubling of atmospheric CO 2 levels to greater than 400 parts per million (ppm) (8–10). Given the apparent timing of both events, many studies have suggested that the environmental perturbations of the MMCO are connected to CRBG eruptions (1, 8, 11), but high-precision geochronologic data that link the two events are lacking (7). While most flood basalt provinces are thought to relate to mantle plume–driven processes, models for the origin of the CRBG include a mantle plume source (12) as well as subduction-related processes such as slab tear (13) or slab rollback (14). Central to this debate is the fact that the CRBG was erupted during a period of active regional volcanism, including subduction volcanism of the Cascade arc, rhyolitic volcanism of the Yellowstone–Snake River Plain hotspot track, the High Lava Plains of central Oregon, and bimodal volcanism related to Basin and Range extension in northern Nevada (Fig. 1). The CRBG began erupting from a north-trending linear fissure system in eastern Washington, eastern Oregon, western Idaho, and northern Nevada, in a back-arc setting between the Cascades and Rocky Mountains. Volcanism progressed from south to north, and the regional tectonic and geomorphic setting guided flows hundreds of kilometers from east to west (1). An improved understanding of the timing of CRBG eruptions, the volumetric rates at which they were emplaced, and the rate at which they propagated geographically through the province are essential constraints on geodynamic models for the CRBG origin. Fig. 1 Map of CRBG and regional volcanism. The map shows the areal extent of each formation of the CRBG, and the legend provides the volume contribution of each formation. Stars represent geochronology sample collection sites; dashed lines enclose areal extent of source dike swarms. The Prineville Basalt (PVB) and Picture Gorge Basalt (PGB) are coeval with the Grande Ronde Basalt, represent 1.4% of the total CRB volume, and are grouped with the Grande Ronde Basalt for all volume estimates presented here (1). Here, we aim to test the hypothesis that there is a temporal relationship between CRBG eruptions and the MMCO by establishing an accurate and precise age model for the eruption of the CRBG. Decades of study have produced a high-resolution stratigraphic framework for the CRBG. The 350 tholeiitic basalt to basaltic andesite flows of the CRBG are divided into five formations: the Steens (31,800 km3; 15.3% of the total volume), Imnaha (11,000 km3; 5.3%), Grande Ronde (150,100 km3; 72.3%), Wanapum (12,175 km3; 5.9%), and Saddle Mountains Basalts (2424 km3; 1.2%), which are composed of a total of 43 stratigraphic members, containing 1 to 20 lava flows each (1). Magnetic field reversals were ongoing during the eruption of the CRBG. While the magnetic stratigraphy first was developed in the field using a portable fluxgate magnetometer during mapping efforts (1), some detailed modern paleomagnetic data have since been published (15, 16). The Steens Basalt erupts during a reversal from reversed to a normal polarity interval [sometimes referred to as R 0 and N 0 (17)], which continues through the Imnaha Basalt. The Grande Ronde Basalt is marked by two couplets of magnetic field reversals (locally defined magnetostratigraphic units are known as R 1 , N 1 , R 2 , and N 2 ). The N 2 normal magnetozone continues through the majority of the Wanapum Basalt, which exhibits a final reversed interval (1, 18, 19). The majority of ages published on the CRBG have uncertainties that preclude the development of an unambiguous chronology for the timing and duration of CRBG volcanism. For example, in a review of K-Ar and 40Ar/39Ar geochronology for the CRBG, a preferred chronology was developed with eruption of the Steens Basalt at 16.9 to 16.7 Ma, the Imnaha Basalt at 16.7 to 16.0 Ma, the Grande Ronde Basalt at 16.0 to 15.6 Ma, the Wanapum Basalt at 15.6 to 15.0 Ma, and the Saddle Mountains Basalt in several distinct events between 15 and 6 Ma (4). This eruptive age model is based on geochronologic analyses with large uncertainties (>1 Ma) that make adherence to stratigraphic order difficult to address, and the result is that it is inconsistent with the geomagnetic polarity time scale (GPTS) (4, 20). For instance, this age model suggests a normally magnetized interval lasting ~600 thousand years (ka) through the Imnaha Basalt, which exceeds the duration of any normal chrons occurring around 16 Ma in different calibrations of the GPTS (20). Other inconsistencies are discussed in greater detail in Results. New high-precision 40Ar/39Ar dates derived from feldspar phenocrysts in silicic tuffs interbedded in the Steens Basalt have revised its eruptive duration to ~16.75 to 16.54 Ma and propose that the Steens magnetic field reversal occurred at 16.603 ± 0.028 Ma (21) [dates recalculated with Fish Canyon sanidine age of Kuiper et al. (22); see the Supplementary Materials]. These new data highlight the potential that additional precise geochronology has to resolve the timing and duration of CRBG volcanism and its correlation to the GPTS and MMCO.

RESULTS Here, we use U-Pb zircon geochronology by chemical abrasion–isotope dilution–thermal ionization mass spectrometry (CA-ID-TIMS), which can achieve the precision (~0.1%) and accuracy required to address the issues with previous age models outlined above. Because basalts are generally too low in Si and Zr to saturate zircon, we follow the sampling strategy of Schoene et al. (23) by collecting silicic volcanic ash from between basalt flows and dating them by U-Pb geochronology on single zircon crystals. These ash beds were sourced from contemporaneous regional silicic volcanism, and deposits include ash-bearing paleosols (referred to as “redboles” in the remainder of this paper), interflow volcaniclastic sediments, and pumice-bearing airfall tuff (Fig. 2 and fig. S1). Zircons separated from these samples (fig. S2) were dated by CA-ID-TIMS at Princeton University. Fig. 2 CRBG geochronology samples. Certain interflow lithologies were targeted for sampling and found to contain magmatic zircons, allowing high-precision dating techniques to be applied: (A) redboles, analogous to those found in the Deccan Traps (similar to CRB1519, CRB1556, CRB1624, CRB1625, and CRB1634); (B) ashes deposited from Cascade subduction volcanism or other silicic centers (CRB1506 and CRB1586); (C) ash-bearing sediments trapped in brecciated basalt flow tops; and (D) pumice clasts in the Vantage sedimentary interbed (CRB1533); the horizon with these centimeter-scale clasts is outlined in white (photo credit: Jennifer Kasbohm, Princeton University). Data from eight horizons within the CRBG, reported as 206Pb/238U dates from single zircons taken from a population of dated grains with 95% confidence intervals, are shown in Fig. 3. Dates on single crystals within each sample spread beyond analytical uncertainty, a common observation in volcanic ash beds, because zircon can retain radiogenic Pb at temperatures of >900°C. A spread in crystallization dates beyond analytical uncertainty can reflect growth of zircon in the magmatic system before eruption, incorporation of pre-eruptive zircon from the volcanic edifice during eruption, or inheritance from the rock hosting the magmatic system (24, 25). In light of these potential sources for complex zircon populations deposited within the CRBG ashes, we use the youngest single concordant analysis, excluding low-precision analyses, as an estimate for the age of the ash (Fig. 3) (26). This approach is tested by dating a large number of zircons from each sample to obtain multiple analyses that overlap with the youngest grain and further requiring that samples from stratigraphically higher positions in the CRBG are in chronological order. A growing database of zircon geochronology from a variety of tectonic settings, using a variety of U-Pb analytical methods, shows that the youngest zircon analysis from continuous age population (excluding outliers that can be explained by residual Pb loss) overlaps with independent eruptive ages (27, 28). We also show in the Supplementary Materials that our use of single youngest zircon ages overlaps with weighted mean ages calculated from multiple grains for each sample, although we prefer the youngest zircon age as the more conservative interpretation. Our dates from the Steens Basalt show excellent agreement with ages derived from high-precision 40Ar/39Ar sanidine geochronology (Fig. 3) (21). Fig. 3 U-Pb zircon CA-ID-TIMS geochronological data. Rank order plot of geochronological data presented in this study, with sample names and the youngest, most precise zircon age starred and labeled. Error bars represent 2σ uncertainty, and stratigraphic younging is from left to right. We present a timeline, shown with background shading and labels, for the eruption of each formation based on the stratigraphic position of these samples and consistent with our geochronology results. Steens Basalt geochronology is compared to the most recent 40Ar/39Ar geochronology results (21) and is found to overlap U-Pb data when recalculated with the Fish Canyon Sanidine age of Kuiper et al. (22). Our zircon ages improve estimates for the timing and duration of each formation of the CRBG stratigraphy (Fig. 3) and indicate that 95% of its total volume erupted in 758 ± 66 ka (fig. S3). The upper 72% of the Steens Basalt volume erupted between 16.653 ± 0.063 Ma and 16.589 ± 0.031 Ma. The latest Steens eruptions occurred concurrently with Imnaha Basalt eruptions, which we have dated at 16.572 ± 0.018 Ma. While we currently do not have an estimate for the onset of the overlying Grande Ronde Basalt, four samples from its upper half fall in stratigraphic order and show its termination by 16.066 ± 0.040 Ma. Two samples from the bottom and top of the Wapshilla Ridge Member of the Grande Ronde Basalt, which comprises 20% of the total CRBG volume, gave dates overlapping within uncertainty of 16.288 ± 0.039 Ma and 16.254 ± 0.034 Ma. These ages come from a single stratigraphic section, and the age of the upper bound is indistinguishable at the 95% confidence interval from the age from a redbole that is in the same stratigraphic position but is found 44 km away (16.210 ± 0.043 Ma). Finally, the lower 77% of the Wanapum Basalt finished erupting before 15.895 ± 0.019 Ma (Fig. 3). Our dated samples bracket 95% of the eruptive history of the CRBG, and these data show that the CRBG erupted 2.4 times faster than the previous estimates (4). Our ages agree with the relative chronology of the existing stratigraphic framework (1) and bolster regional correlation through geochemical and paleomagnetic data. Combined with detailed volume estimates for each member of the CRBG (1), our new zircon ages yield effusion rates throughout the eruptive history of the CRB and provide a timeline for the release of CO 2 , SO 2 , and other gases into the atmosphere. For the main phase of eruption (Steens, Imnaha, and Grande Ronde Basalts), we calculate an average effusion rate of 0.334 ± 0.042 km3/year, which is double the previous estimate of 0.178 km3/year for the same interval (12). During the eruption of the Wanapum Basalt, the rate slows to 0.055 ± 0.014 km3/year, whereas during the Wapshilla Ridge Member, the average rate calculated is 1.18 km3/year, with a minimum rate of 0.376 km3/year. However, because our ages for the top and bottom of this voluminous member in a single stratigraphic section overlap, effusion rates were likely far greater than the average values. A potential constraint for the maximum rate can be derived from previous estimates for the minimum duration of the 1300 km3 Roza Member to be 14 years (29), and the ca. 100 years or more needed to develop thin redbole horizons (30). Because there are at least 18 (and as many as 28) lava flow units in the Wapshilla Ridge Member (31) and 8 more redboles between the samples we dated, we estimate that peak eruption rates could be as high as 20 to 40 km3/year for a duration of thousands of years. The eruptive pulse for the Wapshilla Ridge Member is an order of magnitude higher than average effusion rates calculated for other large igneous provinces associated with mass extinction events: 1 to 2 km3/year for the Deccan Traps (23), 3 to 5 km3/year for the Central Atlantic Magmatic Province (32), and 1 to 4 km3/year for the Siberian Traps (33). However, records from other flood basalts do not have the stratigraphic resolution to calculate effusive rates during high flux pulses compared to long-term average rates. More generalized models for the total climatic impact during flood basalt eruption will be dependent on delineating the timing of eruption pulses, such as the Wapshilla Ridge Member, from background fluxes.

DISCUSSION Our age model provides quantitative constraints that must be satisfied by any geologic or geodynamic model for CRBG volcanism. In particular, the mechanism must be consistent with (i) eruption duration of ~750 ka from 16.65 to 15.90 Ma; (ii) an average effusion rate of 0.334 ± 0.042 km3/year, with pulses of >>1 km3/year; (iii) simultaneous eruptions at Steens Mountain and in the Imnaha Basalt vents 300 km away; and (iv) an average linear geographic propagation rate of eruption of 0.37 ± 0.08 m/year to the north, given distances between the vents sourcing Imnaha through Wanapum eruptions. These criteria alone may be currently insufficient to fingerprint a mantle plume or subduction-related origin of the CRBG, both of which allow eruptions to occur in this time frame (13, 34). The geographic propagation rate of 0.37 ± 0.08 m/year for CRB volcanism is also compatible with either model: A small plume head has been modeled to spread at 0.2 to 0.3 m/year (35), while the proposed slab tear is modeled to propagate at 0.45 m/year (13). This northward propagation rate is about three times faster than those calculated for the McDermitt and High Rock dike swarms that propagate to the south of Steens Mountain (0.12 and 0.14 m/year, respectively) (36), showing that the potential radial propagation from Steens Mountain did not occur at the same velocity radially. Further modeling with our new quantitative constraints is required to better understand the process that allowed CRBG eruption. Determining the relative timing of CRBG volcanism and the MMCO requires independent chronologies that are equally precise. However, the early- and mid-Miocene is one of the most problematic periods in the Neogene for establishing precise independent chronologies in marine sediments due to the difficulty in obtaining undisturbed stratigraphic sections that yield reliable magnetostratigraphy, biostratigraphy, astronomical tuning, and radiometric ages (20). All time scales proposed for the mid-Miocene depend directly or indirectly on correlation with the GPTS, for which there are currently several proposals, the most recent being the Geologic Time Scale (GTS) 2012 (20). GTS 2012 was derived from the seafloor anomaly profiles of the Antarctic and Australian plates and assuming a relatively constant spreading rate tuned to give a 23.03-Ma age for the Oligocene-Miocene boundary (20). GTS 2012 rejected an astronomically tuned record of mid-Miocene δ18O and magnetostratigraphy from ODP (Ocean Drilling Program) site 1090 in the subantarctic south Atlantic, whose record extends from the Oligocene-Miocene boundary to ~15.9 Ma, because the tuned record yields ages for chron boundaries that do not meet the assumption of constant seafloor spreading rates in the Pacific (20). The most recent age model for the CRBG (4) attempts to reconcile 40Ar/39Ar geochronology with GTS 2012. However, the resulting age model is inconsistent with the existing GPTS and is in need of refinement (Fig. 4). Fig. 4 Revised CRBG eruptive timeline, magnetostratigraphy, and GPTS correlation. U-Pb geochronology suggests a timeline of the eruption for each formation (from Fig. 3) as well as a revised GPTS consistent with CRBG magnetostratigraphy [(1) and references therein]. These are compared to the eruptive chronology derived from 40Ar/39Ar geochronology (4) and different GPTS calibrations (20, 37, 46). Given the magnetic polarity of different stratigraphic members, U-Pb geochronology constrains the age of four different chron boundaries (straight lines), identified with arrows and ages with internal and decay constant uncertainties. Estimated chron boundaries are shown with zigzag lines and are not yet constrained by geochronology. Lighter shades of color in the stratigraphic column represent reversed polarity intervals in CRBG magnetostratigraphy, also shown by adjacent reversal stratigraphy to the right of each CRBG age model. Stars indicate the youngest zircon ages obtained for each sample in the study, and letters label each formation (S, Steens Basalt; I, Imnaha Basalt; GR, Grande Ronde Basalt; W, Wanapum Basalt). The blue diamond represents the age of the Steens reversal obtained from Mahood and Benson (21), recalculated with the Fish Canyon Sanidine age of Kuiper et al. (22), to be 16.603 ± 0.028/0.36 Ma, which is consistent with our results. Our new age model for the CRBG permits a more robust correlation of CRBG magnetostratigraphy with existing proposals for the GPTS (Fig. 4). However, this exercise also indicates that some previous proposals for the GPTS, including GTS 2012, are in error. For example, the most recent age model for the CRBG has the Imnaha Basalt, which is entirely normally polarized, erupting through several magnetic reversals and thus is not permissible. Similarly, the existing age model places the Grande Ronde Basalt, which records two reversed and two normal intervals, within a single normal chron. By comparison, in our proposed correlation illustrated in Fig. 4, the Imnaha Basalt erupted entirely during chron C5Cn.3n, while the Grand Ronde Basalt erupted during C5Cn.2r-C5Cn.1n, consistent with observed magnetostratigraphy in the basalts. Using this baseline correlation with the GPTS, we can refine four proposed reversal ages (Fig. 4). Our ages in the Upper and Lower Steens bracket the “Steens Reversal” (between magnetozones R 0 and N 0 , and chrons C5Cr and C5Cn.3n), which can be conservatively constrained to 16.637 ± 0.079/0.089 Ma (95% confidence intervals given for internal uncertainty/decay constant uncertainty). This estimate compares favorably with the estimate of 16.603 ± 0.028/0.36 Ma obtained through recent 40Ar/39Ar sanidine geochronology (21). Our samples from the base and top of the Wapshilla Ridge Member constrain the timing and provide a minimum duration for C5Cn.1r, to begin no later than 16.288 ± 0.039/0.046 Ma, and to end no earlier than 16.210 ± 0.043/0.047 Ma, because the Wapshilla Ridge Member comprises the majority of volume of the second reversed magnetostratigraphic unit of the Grande Ronde Basalt (R 2 ) (31). The end of C5Cn.1n (N 2 ) is well constrained by our age of 15.895 ± 0.019/0.026 Ma for the top of the transitionally magnetized Roza Member, which immediately overlies the normally magnetized Frenchman Springs Member, especially given previous estimates that the Roza Member erupted in as little as 14 years (29). Our initial data do not identify any significant hiatuses in eruptions—no more than ~200 ka elapse between any two of our samples, during which volcanism is known to be ongoing, although we do not present zircon data from these intervals (fig. S3). Therefore, high-precision geochronology can be used to bound the ages of magnetically characterized CRB flows and to further refine the record of mid-Miocene magnetic field reversals. Our proposed GPTS is also consistent with the astronomically derived age model for the magnetic reversal stratigraphy at IODP (Integrated Ocean Drilling Program) site U1335 in the equatorial Pacific (Fig. 4) (37), indicating an independent verification for our proposed age model for the GPTS. Given the inconsistencies described above for the GPTS, demonstrating a link between the eruption of the CRBG and the MMCO requires a careful assessment of the age models used to develop proxy records across the MMCO. For example, the δ11B proxy record for pco 2 at ODP site 761 indicates that atmospheric CO 2 increases at 16.5 Ma (8), which agrees well for our suggested timing of the onset of voluminous Grande Ronde Basalt volcanism. However, the age model for site 761 (38) depends on biostratigraphic (39) or isotopic events (40) tied to calibrations of the GPTS (41) that we have shown to be inaccurate. Recent work describing the δ13C and δ18O records from IODP site U1337 identifies the onset of the MMCO at 16.9 Ma (42), which precedes our timing for all CRBG eruptions. This site has an age model derived from an astronomical solution (43) without radiometric age control or a magnetostratigraphy, thus adding subjectivity to the chosen isotopic tie points used to calibrate the tuning (44) and making correlation with our eruptive record difficult. One way forward is to use proxy records from sites that contain reliable magnetostratigraphy (37, 45). Benthic δ18O values—a proxy for deep-ocean temperature—from sites 1090 (46) and U1335 (37) (Fig. 5) indicate that the decline in δ18O values began during what is interpreted as C5Cr, reaching a nadir (the MMCO) during C5Cn.3n-C5Cn.1r. While it is currently difficult to validate the identification of C5Cr from site 1090 given the potential for hiatuses in the record, it is corroborated by the astronomical model of U1335 and is interpreted to be the same chron in which CRB eruptions began with the Lower Steens Basalt. Although the absolute timing of the onset of the MMCO cannot be confirmed by our data, the astronomically tuned record from U1335 compared to our geochronology shows that the decrease in δ18O preceded the eruption of Steens Basalt lava flows by 100 to 200 ka. CRBG volcanism may have played a role in eliciting global warming through cryptic degassing of CO 2 as magma migrated through dike swarms before surface eruptions (7). Alternatively, the apparent mismatch between the onset of the CRBG and the MMCO may indicate that the two events are unrelated. Regardless, the δ18O minimum appears coeval with the eruption of the Grande Ronde Basalt, suggesting that a link may exist. Further work that refines age models for climate proxy records across the MMCO and investigates rates of eruptions for the CRBG is required before it can be determined whether the CRBG caused the MMCO. While our work constrains the age of the top of chron C5Cr, there are not yet absolute age constraints either for the bottom of C5Cr or for the onset of the MMCO. Because the record for U1335 has an astronomical age model constrained by isotopic correlations with site U1337 (37, 42) and does not have an absolute age defining the base of chron C5Cr, a tighter correlation of the CRBG and MMCO could be observed, particularly if the chron began later in time than is currently proposed. These uncertainties in the timing of this magnetic field reversal and the beginning of the MMCO must be resolved to better assess whether the CRBG played a causative role in the MMCO. Fig. 5 Correlation of the CRBG with the MMCO. (A) A compilation of proxy records exhibiting the MMCO (47), with age constraints as reported in each study. Although ages are susceptible to uncertainties in the mid-Miocene time scale, the magnitude of the isotopic signals is not. (B) To compare zircon geochronology results for CRBG eruptions to paleoclimate proxy records of the MMCO, it is necessary to bypass age models tied to outdated calibrations of the GPTS. The robust magnetostratigraphy of sites 1090 (45, 46) and U1335 (37) allows correlation of these isotopic records to our CRBG eruption chronology and refined GPTS. The area of each colored rectangle corresponds to the volume of each formation (1) (S, Steens Basalt; I, Imnaha Basalt; GR, Grande Ronde Basalt; W, Wanapum Basalt), with width constrained by zircon ages (slanted boundary indicates that the onset of Steens Basalt volcanism is not yet constrained); polarity of the basalt flows is taken from Reidel (1) and references therein. Yellow shading compares global proxy data at 17 to 16 Ma (lacking an age model based on absolute geochronology) with volcanic events occurring 17 to 16 Ma, while the light blue shading highlights the onset of the MMCO in both records with the drop in δ18O. Despite uncertainties present in mid-Miocene age models, global proxy data (47) indicate that the MMCO continued for >1 Ma after the cessation of the majority of CRBG volcanism (Fig. 5). The time lag between the cessation of volcanism and a return to cooler climatic conditions could be understood as a consequence of the long response time of negative feedbacks within the global carbon cycle that regulate atmospheric CO 2 and Earth’s temperature on geologic time scales. These feedbacks include interactions between temperature, the chemical weathering of continental silicate minerals, and the burial of CO 2 in marine carbonate sediments (48). While the sensitivity of the silicate weathering feedback remains poorly understood, recent estimates for response times vary from ~200 to 500 ka (49) and are consistent with the stabilization of atmospheric CO 2 (that is, return to baseline conditions) on ~1 million year time scales. Our age model of CRBG emplacement shortens the duration of volcanism from 1.9 Ma (4) to 750 ka and correlates the onset of CRBG volcanism and the onset of the MMCO to within ~100 ka. A shorter duration of CRBG volcanism implies higher average CO 2 emissions and higher peak CO 2 concentrations during volcanism, to be compared with marine proxy records. However, current proxy records for atmospheric CO 2 during the MMCO are too coarse for a close comparison to the eruptive history of the CRBG, further inhibiting the ability to assess whether or not the CRBG caused the MMCO. Furthermore, establishing a quantitative link between CRBG volcanism and changes in the global carbon cycle and atmospheric CO 2 is hampered by uncertainties in the amount of CO 2 emitted by flood basalts from dissolved mantle carbon in addition to “cryptic” sources, such as organic or inorganic sediments volatilized through contact with basaltic flows or sills (7). Armstrong McKay et al. (7), using a main phase CRBG eruptive duration of 900 ka, model that 4090 to 5670 Pg of emitted carbon can yield the observed changes in benthic δ13C and atmospheric CO 2 , although this amount includes a substantial component of cryptic degassing beyond the expected volatile release of subaerial basalt flows. Future studies should focus on further revision of the mid-Miocene time scale and a high-resolution climate proxy record spanning the 700-ka duration of CRBG volcanism to explore the extent to which the timing of CRBG volcanism agrees with changes in atmospheric CO 2 . Such studies will lead to an improved understanding of the MMCO and more general models linking volcanism to climate change and could be crucial for understanding why some flood basalts apparently result in mass extinctions and others do not.

MATERIALS AND METHODS Methods are as described by Schoene et al. (23) and Samperton et al. (50). Zircon separation and preparation Zircons were separated from their host rock through standard methods of crushing as well as gravimetric and magnetic separation techniques using a Bico Braun “Chipmunk” jaw crusher, disc mill, hand pan, hand magnet, Frantz isodynamic separator, and methylene iodide. Zircons from the least magnetic and most dense mineral separate were transferred in bulk to quartz crucibles and annealed in a muffle furnace at 900°C for 48 hours after Mattinson (51). After annealing, 20 to 40 zircon grains from each sample were photographed (fig. S2) and picked in reagent-grade ethanol for analysis. Given the low radiogenic Pb content of the samples, cathodoluminescence images were not obtained. Euhedral grains with a range of morphologies were selected, while those with visible cracks, inclusions, and cores were avoided. Individual grains were transferred using stainless steel picking tools to separate 3-ml Savillex Hex beakers containing distilled acetone and taken to the clean laboratory for analysis. U-Pb zircon ID-TIMS analysis Single zircon grains were loaded into 200-μl Savillex “micro”-capsules with 100 μl of 29 M HF + 15 μl of 3N HNO 3 for a single leaching step in high-pressure Parr bombs at 185°C for 12 hours to remove crystal domains affected by Pb loss (51). Grains were rinsed after leaching in 6 N HCl, MQ H 2 O, 3N HNO 3 , and 29 M HF before spiking with EARTHTIME (202Pb-)205Pb-233U-235U tracer and addition of 100 μl of 29 M HF + 15 μl of 3N HNO 3 (52, 53). Zircons were then dissolved to completion in Parr bombs at 210°C for 48 hours. Dissolved zircon solutions were subsequently dried down, dissolved in 100 μl of 6N HCl, and converted to chlorides in Parr bombs at 185°C for 12 hours, after which solutions were dried again and brought up in 50 μl of 3N HCl. The U-Pb and trace element aliquots were then separated by anion exchange chromatography using 50-μl columns and AG-1 X8 resin (200 to 400 mesh, chloride from Eichrom) (54) and dried down with a microdrop of 0.015 M H 3 PO 4 . The dried U and Pb aliquot was loaded in a silica gel emitter (55) to an outgassed zone-refined Re filament. Isotopic determinations were performed using an IsotopX PhoeniX-62 TIMS at Princeton University, with Pb analysis performed in peak-hopping mode on a Daly photomultiplier ion-counting detector. A correction for mass-dependent Pb fractionation was applied in one of two ways. For double-Pb spiked analyses (202Pb-205Pb, ET2535), a cycle-by-cycle fractionation correction was calculated from the deviation of measured 202Pb/205Pb from the known tracer 202Pb/205Pb [0.99924 ± 0.00027 (1σ)]. For single-Pb spiked analyses (205Pb, ET535), a Pb fractionation of 0.182 ± 0.041% per atomic mass unit (amu) was used, as determined by repeat measurements of NBS982 at Princeton. A Daly photomultiplier dead time of 28.8 ns was used, as determined by repeat measurements of National Bureau of Standards (NBS) standards. Corrections for interfering isotopes under masses 202, 204, and 205 were made cycle by cycle by measuring masses 201 and 203 and assuming that they represent 201BaPO 4 and 203Tl and using natural isotopic abundances to correct for 202BaPO 4 , 204BaPO 4 , 205BaPO 4 , and 205Tl. UO 2 measurements were performed in static mode on Faraday cups with a bulk U fractionation correction calculated from the deviation of measured 233U/235U from the known tracer 233U/235U [0.995062 ± 0.000054 (1σ)], and an oxide composition of 18O/16O of 0.00205 was used (56). Data reduction was performed using the programs Tripoli and U-Pb Redux (57, 58) and the decay constants of Jaffey et al. (59). All Pbc was attributed to laboratory blank with a mean isotopic composition determined by total procedural blank measurements (see table S1 for values). Two different blank models were generated to assess data collected before [Outliers Culled (OC)] and after [Side Filaments (SF)] January 2017, when the laboratory began heating side filaments before collecting data on the mass spectrometer, which was found to reduce interferences. Uncertainties in reported U-Pb zircon dates are at the 95% confidence level and exclude tracer calibration and decay constant uncertainties. Correction for initial 230Th disequilibrium in the 206Pb/238U system was made on a fraction-by-fraction basis by estimating (Th/U) magma using (Th/U) zircon determined by TIMS and a mean (Th/U) zircon-magma partition coefficient ratio of 0.19 ± 0.11, which encompasses the range of values for (Th/U) zircon-magma partition coefficients obtained from glasses from a variety of volcanic settings (60). Uncertainties for the resulting (Th/U) magma were also calculated on a fraction-by-fraction basis, propagating the uncertainty in the (Th/U) zircon-magma partition coefficient. Overall, these corrections for 230Th disequilibrium affected our results by no more than ±10 ka compared to an alternative approach using a constant (Th/U) magma of 3.5 ± 1.0 (see the Supplementary Materials).

SUPPLEMENTARY MATERIALS Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/9/eaat8223/DC1 Supplementary Materials and Methods Fig. S1. Geochronology sample photos. Fig. S2. Zircon photos. Fig. S3. Thickness and volume versus age plots. Fig. S4. Concordia plots for U-Pb ID-TIMS geochronological data. Fig. S5. Alternate age interpretations. Table S1. U-Pb isotopic data. Table S2. Alternate age interpretations. References (61–66)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

Acknowledgments: We thank S. Reidel and B. Martin for assisting with stratigraphic classification of the samples, and J. Murray, S. Bartusek, and K. Duffey for assisting with fieldwork and sample preparation. M. Eddy, J. Higgins, M. Long, and S. Reidel provided thoughtful feedback on an earlier version of this manuscript. We thank S. Burgess and two anonymous reviewers for suggestions that greatly improved the manuscript. Funding: This material is based upon work supported by the NSF Graduate Research Fellowship under grant no. DGE-1656466, by Princeton Environmental Institute at Princeton University through the Walbridge Fund, and by the Princeton University Department of Geosciences Scott Vertebrate Fund. Author contributions: J.K. and B.S. devised the study, conducted fieldwork, interpreted results, and prepared the manuscript; J.K. performed most zircon analyses. Competing interests: The authors declare no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.