The history of calcite diagenesis and origin of exceptionally negative oxygen isotope values in chalks of the Niobrara Formation, Denver Basin, USA

The Niobrara Formation of north‐east Colorado, USA, has anomalously negative δ18O values compared to all other Cretaceous chalks. These unique δ18O values have been attributed to elevated heat flow and/or freshening of the Cretaceous Western Interior Seaway. This work utilises clumped isotopes of calcite (Δ47), peak burial temperatures estimated from pyrolysis data, and strontium and neodymium isotopes of carbonate to re‐evaluate the origin of the calcite's 18O‐depletion. Peak temperatures indicate lateral variability in geothermal gradients of ca 20°C/km at the tens of kilometre scale, and corroborate prior studies proposing locally elevated palaeotemperatures. Greater insight is provided by numerical models of calcite recrystallisation and oxygen isotope evolution that are constrained by measured Δ47‐derived temperatures, calcite δ18O values and inferences from the 87Sr/86Sr and εNd values. The models indicate that (1) sea water in the seaway had normal marine δ18O values of −1 (VSMOW) except on the eastern margin of the basin where some freshwater dilution yielded −2 to −3‰ (VSMOW) water, and (2) the main driver of the anonymously negative calcite δ18O values was a semi‐open hydrologic system that provided a few percent by pore volume of meteoric groundwater derived from post‐Laramide recharge into the basin. Minor contributions were a Laramide‐aged heat pulse related to the underlying Colorado Mineral Belt, the thermal insulating effects of now eroded coals, and a small flux of compaction‐driven Cretaceous sea water evolved by smectite dehydration. However, those three factors alone were insufficient drivers of the calcites' 18O depletion. High burial temperatures are interpreted to have caused clumped isotope reordering in at least one well, but those temperatures cannot yield the observed calcite δ18O values. The study illustrates the unique attributes of the Niobrara's diagenetic system that results in its anomalous δ18O values, and reaffirms the value of clumped isotopes in unravelling the diagenetic history of chalk systems.


| INTRODUCTION
The oxygen isotopic composition of chalk oozes and limestones have long been studied as archives of palaeoceanographic and palaeoclimate trends, and as indicators of diagenetic alteration with progressive burial. The δ 18 O values of relatively shallow buried deep-sea chalks can preserve, albeit with some diagenetic overprinting, the secular patterns of palaeoclimate (temperature and ice volume) and palaeoceanographic water masses (Fantle et al., 2020;Swart, 2015). In contrast, the δ 18 O values of more deeply buried chalks below continental margins and in continental basins show a much stronger diagenetic overprint, expressed as increasingly more negative δ 18 O values with depth ( Figure 1). The progressive decrease in δ 18 O values is attributed to calcite recrystallisation and re-precipitation of pressure-solved CaCO 3 with increased burial stress and temperatures in closed to semi-closed systems (Börre & Fabricius, 1998;Fabricius, 2007;Fabricius & Börre, 2007;Killingley, 1983;Maliva et al., 1991;Schlanger & Douglas, 1974;Scholle, 1977). Deviations in the depth trend in a chalk's δ 18 O values can be informative with respect to a variety of basin phenomena, including fluid flow, spatial and temporal differences in heat flux, cementation by extra-formational CaCO 3 , and uplift or denudation.
The Niobrara Formation in the Denver Basin, USA, is a unique example of a chalk δ 18 O record that deviates from the global trend for deeply buried chalks in petroliferous basins ( Figure 1). Multiple workers have documented δ 18 O values of Niobrara calcite that are significantly lower than other chalk formations worldwide (Hefton, 2015;Pollastro & Scholle, 1986;Scholle, 1977). It remains unclear if the very negative δ 18 O oxygen isotope values of Niobrara calcite is solely the result of elevated burial temperatures (Crysdale & Barker, 1990;Pollastro & Martinez, 1985;Rice, 1984a), 18 O-depleted initial sea waters (Arthur et al., 1985;Pagani & Arthur, 1998;Pollastro & Scholle, 1986;Pratt & Barlow, 1985), an extra-formational fluid with relatively low δ 18 O, or a combination of any of these possibilities.
This study evaluates these historical hypotheses in order to reconcile the disparate interpretations of Niobrara oxygen isotope geochemistry relative to burial and thermal history. Clumped isotope measurements (Δ 47 ) provide new constraints on pore-fluid oxygen isotopic evolution and palaeotemperatures during calcite diagenesis; insights not available from conventional stable isotope data. Geothermal gradients are constrained by vitrinite reflectance values derived from existing pyrolysis data and reconstructed maximum burial depths. Strontium and neodymium isotopes test the possibility of extra-formational fluid fluxes into the Niobrara's hydrologic system. An equilibrium calcite diagenesis model that simulates the evolution of calcite and K E Y W O R D S carbonate diagenesis, chalk, clumped isotopes, Denver Basin, Niobrara Formation, palaeothermometry F I G U R E 1 Measured chalk δ 18 O values versus current burial depth for Cretaceous chalks. Values plotted as averages for individual wells or fields. Data sources: North Sea- Egeberg and Saigal (1991), pore fluid δ 18 O values during recrystallisation and increasing burial temperatures is then used to evaluate potential diagenetic scenarios. Temperatures derived from Δ 47 and calcite δ 18 O values constrain the validity of the model outcomes. The results show that neither elevated temperatures, laterally varied heat fluxes, nor 18 O-depleted sea water in the seaway by themselves are sufficient to explain calcite δ 18 O values. Instead, a diagenetic system that included the influx of 18 O-depleted meteoric fluids into the basin after the Laramide orogeny is the critical factor. All of these new insights showcase the need to cross-validate interpretations of diagenesis in any basin, as an understanding of thermal history, fluid flux and porosity loss are contingent upon an accurate diagenetic framework.

| Sedimentologic and stratigraphic framework
The Niobrara Formation in the Denver Basin, USA, is one of many Laramide-age foreland basins that segmented the greater Cretaceous Western Interior Basin (Figure 2).
The Niobrara Formation contains late Turonian to earliest Campanian (Cobban et al., 2006;Sageman et al., 2014) interbedded chalks and marls deposited in shallow (≤300 m) waters in the eastern reaches of the Western Interior Seaway (Figure 2A; Kauffman, 1977;Kauffman & Caldwell, 1993;Longman et al., 1998). The basal Fort Hays Member (Figure 3) ranges from 5 to 40 m thick in the Denver Basin, and consists of medium to thick bedded, bioturbated, Inoceramid-rich, nearly pure chalk that is interbedded with centimetre-thick laminated marl. The overlying Smoky Hill Member is between 100 and 150 m thick in the study area and is composed of chalk-marl packages at multiple scales. The multi-scale packaging at the centimetre scale consists of alternating marl-chalk laminations that in turn stack into metre-scale, chalkening upwards and marling upwards cycles (Deacon & McDonough, 2018). At the largest scale, decametre-thick marl-chalk sequences are present and informally referred to as Niobrara A, B, C and D intervals ( Figure 3). In each of these thicker sequences, the marl component consists of marl-dominated, metre-scale packages and the chalk component contains chalk-dominated, metrescale packages. Marls at all scales formed during periods of increased terrigenous input with higher preservation of F I G U R E 2 Location maps of the Denver Basin and sample locations. (A) Location of the Denver Basin (red circle) and study area (black dashed box) in the Cretaceous Western Interior Basin (map ©2014 Ron Blakey, Colorado Plateau Geosystems, Inc.). (B) Current position of the study area and the Denver Basin. (C) Detailed map of the study area. Yellow stars and circle denote cored wells and outcrop respectively. Alphabetical well identifiers correspond to sample labels in Tables 1 and 2. Purple shaded area is the 400 m subsurface contour of the Denver Basin. Pink shaded area denotes the greater Wattenberg field area. White square is location of the burial history plot shown in Figure 5; line of section A-A' is shown in Figure 11.

| Burial and thermal history
The Niobrara Formation in the Denver Basin has a different burial history along the basin's axis and eastward F I G U R E 3 Stratigraphic nomenclature, gamma ray log and δ 18 O values for an example of the Niobrara Formation. The gamma ray track is coloured to reflect carbonate content, where blue is more calcareous (chalky), and brown is more siliceous (marly). The centre depth track is thickness in metres. Boundaries of the informal stratigraphic units are tilted in the nomenclature track to represent lateral thickness variations discussed by Longman (2020). Oxygen isotope values demonstrate only minor differences at the unit scale (right isotope track, Berthoud state 3 well retrieved from https://my.usgs.gov/crcwc/ core/repor t/9202). At the scale of centimetre-spaced lamination (far left, magnified isotope insets from three different wells analysed in this study and projected to their approximate stratigraphic position within the Berthoud state well), occasional shifts of more than 1‰ can be observed over a few millimetres regardless of the well's thermal maturity.
relative to the basin's far western margin that abuts the Colorado Front Range, which was subjected to Laramide structural deformation (Higley & Cox, 2007;MacMillan, 1980;Tainter, 1984). Only the former is of interest to this study and that burial history can be broken down into three major segments ( Figure 5). The first segment is rapid late Cretaceous burial during deposition of the overlying Pierre Shale, Fox Hills sandstones and Laramie Formation, which positioned the Niobrara at a depth of 2,250 m or more in the axis of the Denver Basin by the end of the Cretaceous. Sites farther east would have attained lesser depths due to thinner Pierre Shale. In the second burial segment, Niobrara rocks in the basin's axis experienced 600-800 m of additional burial due F I G U R E 4 Niobrara rock fabrics. (A) Thin section photomicrograph of a chalk sample. White material is compacted foramanifera; black are stylolites and/or elongate organic macerals. Light brown is calcite and the brighter, elongate honey brown ovoids are squashed peloids. (B) Scanning electron micrograph showing a calcitic peloid (outlined in red) now composed of euhedral calcite crystals whose intercrystalline pores are infilled with bitumen. White circular objects are pyrite (p) framboids; also present are organic macerals (om). (C, D) Back-scattered electron image (left) and elemental map (right) of the same field of view in a chalk sample that is ca 94% calcite. Blue is Ca (calcite), red is Si (quartz), black is pyrite, organic macerals or bitumen, and yellow to yellowish greens are Al-Si minerals (clays, feldspars). Ovoid peloids of coccolith debris are hard to distinguish from coccolith matrix, although some peloids are more clearly defined by surrounding silicates and examples are outlined in white. (E, F) Back-scattered electron image (left) and elemental map (right) of the same field of view in a marl sample. Colours as noted for image (D). The marl contains calcitic peloids (blue) in a dominantly clay (yellows) matrix.
to deposition of Palaeocene alluvial deposits (Higley & Cox, 2007), and thereafter remained at a relatively uniform depth to the end of the Eocene. The last burial segment involved deposition of less than a hundred metres of late Cenozoic sediment, followed by erosional uplift over the last 7 million years of 300-1,200 m depending on location within the basin (Higley & Cox, 2007).
Later Niobrara calcite diagenesis was dominated by chemical compaction, which is evidenced by throughgoing microstylolites ( Figure 6C) and volume and shape changes to Niobrara foraminifera ( Figure 6D). Pressure solution drove calcite dissolution due to locally elevated stress at grain-to-grain contacts and/or at pressure solution seams (Pollastro & Scholle, 1986;Scholle, 1977). Calcite liberated by pressure solution was subsequently F I G U R E 5 Niobrara Formation burial history plot (dark grey) for a well located in Silo Field (white box on Figure 2C) along the basin axis (modified from Landon et al., 2001). The burial history has three main segments, indicated by numbered arrows.

F I G U R E 6
Evidence of calcite diagenesis within the Niobrara chalks. (A, B) Transmission electron microscope (TEM) images. To the left is a high-angle annular darkfield (HAADF) image that depicts atomic mass, and to the right is a darkfield 4 image that indicates crystal orientations. The light grey colour in the HAADF image indicates all crystalline material is calcite; pores are dark grey to black depending on their depth. Crystals outlined in yellow consist of radial prisms, each at a slightly different orientation, and with many of the crystals exhibiting epitaxial rhombic terminations. These are interpreted as recrystallised fragments of coccolith plates. (C) Elemental map (colours as described for Figure 4D) showing microstylolites (white arrows) that cross-cut blue calcitic matrix, peloids and foraminifera (f). The stylolites are marked by accumulations of quartz (red) and clay minerals (yellow to yellowish greens). (D) Cathodoluminescence image of a foraminifera shell (black arrows) that has been truncated by pressure solution with the surrounding aphanocrystalline calcite matrix (darker red). Dashed white lines mark the probable original extent of the shell walls, assuming circular chambers. Interparticle pore space within the shell is filled with yellow luminescent calcite cement. (E, F) TEM high-angle annular darkfield (HAADF) image (left) and a darkfield 4 image (right). The light grey colour in the HAADF image indicates all crystalline material is calcite; pores are black. Crystals outlined in thin yellow lines are heterogeneous and composed of many subcrystals with varied orientations in the darkfield 4 image. The heterogeneity is interpreted to reflect original coccolith debris at varied orientations that recrystallised, were overgrown with authigenic calcite and amalgamated into larger calcite crystals. Arrows denote overgrowths on other crystals or original debris. Yet other crystals are homogeneous (hc) with no internal subcrystals of varied orientation. precipitated as calcite overgrowths and cement that also yielded progressively larger amalgamated crystals ( Figure 6E,F). Recrystallisation and pressure-solution related cementation converted the shallow-buried, high porosity Niobrara chalks with excellent grain preservation ( Figure 7A) into low porosity (10% or less) 'welded' fabrics  that consist of euhedral to subhedral masses of interlocking rhombic crystals ( Figure 7B). All prior workers have assumed that all authigenic Niobrara calcite was internally sourced (Pollastro & Scholle, 1986;Scholle, 1977). Unlike some other chalk units, no diagenetic geobodies or cementation fronts (Fabricius, 2007;Smit et al., 2018Smit et al., , 2022 that are characterised by sharp and significant decreases in porosity and more negative δ 18 O values have been observed in the Niobrara Formation. Previously reported δ 18 O values of calcite in the Niobrara (Figures 1 and 3) are commonly between −4 and −11‰ VPDB (Vienna-Pee Dee belemnite), with relatively small intrawell variability but large, depth-related interwell variations (Hefton, 2015;Pollastro & Martinez, 1985;Pollastro & Scholle, 1986). The increasingly negative δ 18 O values with depth is attributed to calcite cementation related to recrystallisation and pressure solution at progressively higher temperatures. One explanation for the most negative δ 18 O values is precipitation at higher temperatures than modern day burial depths suggest (Crysdale & Barker, 1990;Hefton, 2015;Pollastro & Martinez, 1985;Rice, 1984a). Crysdale and Barker (1990) in particular argued that portions of the Niobrara were 30-50°C hotter at maximum burial than they are today. Arthur et al. (1985) also suggested that some of the calcites' 18 O-depletion may have been influenced by clay-related diagenetic fluids.
Another explanation for the negative δ 18 O values is that the Cretaceous sea water in the seaway was lower than the time-equivalent global ocean, which would have meant initially more negative δ 18 O values for the biogenic sediment relative to other Cretaceous chalks (Arthur et al., 1985;Pagani & Arthur, 1998;Pratt & Barlow, 1985). Arthur et al. (1985) calculated that the 18 O-depleted calcite in the Cenomanian-Turonian Bridge Creek Limestone would be in isotopic equilibrium at 25°C with a −3‰ (VSMOW; Vienna standard mean ocean water) Cretaceous sea water, a similar result was also obtained by Slingerland et al. (1996). More recently, however, Petersen et al. (2016) and Jones et al. (2022) used Δ 47 thermometry of molluscs from the Campanian and Cenomanian to argue that the seaway contained normal marine sea water of −1‰ (VSMOW) except in nearshore settings influenced by freshwater.

| ANALYTICAL METHODS
Thirteen wells from Colorado, Nebraska and Kansas and an outcrop in western Kansas ( Figure 2C) were sampled for this study. Six are deep wells (Niobrara currently at >2,100 m) in the greater Wattenberg area, and five of those are proprietary wells (referenced as wells L, M, N, O, P). Four others are shallow (Niobrara currently at >550 m) wells from Nebraska and Kansas, and two are wells in eastern Colorado where the Niobrara is at intermediate depths ( Table 1). Cores of the non-proprietary wells were sampled at the U.S. Geological Survey's Core Research F I G U R E 7 Scanning electron micrographs of Ar-milled surfaces depicting the range of calcitic fabric in Niobrara chalk. (A) Porous chalk from a shallow buried B chalk sample (Schock Errington 1 core, 540.7 m). Mechanical compaction reduced porosity to ca 35% and resulted in a slight horizontal orientation of grains, but there is minimal, if any, calcite cement between the skeletal debris (image accessed from http://my.usgs.gov/crcdata/prod/core/9383/core_analy sis_139751_B518S EM.pptx). (B) Deeply buried (well M, 2,177 m) B-chalk consisting of euhedral to subhedral rhombic calcite crystals in a 'welded' fabric, but still with micro intercrystalline pores (black). Those pores are now filled with porous bitumen.
Center, Denver, CO. A total of 246 samples spanning current burial depths of 0 to 2,376 m were collected and represent the basal Fort Hays, B and C chalks.

| Stable oxygen and carbon isotope measurements
Powdered whole-rock subsamples were collected using a fine-tipped dental drill. Powders may have aggregated some foraminifera and their intraparticle cement with the coccolith-dominated matrix but Inoceramid fragments were avoided. Approximately 110 μg of each sample was weighed into 15 ml Exetainer vials, capped, purged with helium, heated to 70°C, and then acidified via anhydrous phosphoric acid for a minimum of 45 min to liberate CO 2 gas. The evolved CO 2 gas was analysed at the CU Boulder Earth Sciences Stable Isotope Laboratory (CUBES-SIL) using a Thermo Delta V magnetic sector isotope ratio mass spectrometer attached to a GasBench headspace-sampling peripheral. Measured isotope ratios were corrected to the VPDB reference scale using in-house and international carbonate reference materials, and are reported using delta notation (δ 13 C and δ 18 O) in permille. Precision on repeated measurements of standards is better than 0.1‰ for both δ 13 C and δ 18 O. The total dataset is in File S1.

| Clumped isotopes of CO 2
Calcite in 16 Niobrara B, C and Fort Hays chalk samples from the shallow Kansas well and four of the five proprietary wells were analysed for clumped isotope thermometry. This measurement determines the abundance of CO 2 molecules with a mass of 47 amu that contain bonded 13 C and 18 O (Eiler, 2007). The abundance of the mass 47 isotopologue relative to a stochastic distribution of 13 C and 18 O during calcite precipitation (reported as Δ 47 values) reflects temperature of crystallisation and is independent of the δ 18 O and δ 13 C of the fluid from which the carbonate precipitated (Ghosh et al., 2006). This independence makes it possible to calculate δ 18 O values of the parent fluid for each sample using the Δ 47 -derived temperature and the calcite's δ 18 O value measured simultaneously with Δ 47 .
Samples' Δ 47 were measured using a Thermo 253+ mass spectrometer coupled to an automated CO 2 gas preparation device. Six to eight milligrams of micro-drilled powder were acidified in a common acid bath at 90°C. The resulting CO 2 gas was purified to remove water, sulphur species, non-condensables and potential isobaric contaminants. Three replicates of each individual sample were analysed, with each replicate composed of nine separate and sequential measurements ('acquisitions') of sample gas.
Each acquisition, in turn, was composed of seven cycles of alternating purified sample and reference ('working') gasses with 26 s integration times. Each analytical session began with heated (1,000°C) and equilibrated (25°C) reference frame gasses, then ETH1, ETH2 and ETH3 standards in random order, and then blocks of Niobrara samples and randomly interspersed international standards (ETH1-3, IAEA-C1, IAEA-C2 and MERCK). Standard error was better than 0.04‰ for δ 13 C, 0.13‰ for δ 18 O, 0.18‰ for δ 47 and 0.02‰ for clumped isotope ratios in Δ 47 notation.
The raw data were screened for an excess of Δ 48 relative to the gas reference frame as an indicator of potential isobaric interferences on masses 44-47 due to impure CO 2 (Huntington et al., 2009). Analyses with Δ 48 values greater than the Δ 48 values of the heated gasses plus the theoretical spread between 1,000 and 25°C for Δ 48 of CO 2 (0.345‰; Wang et al., 2004) were categorised as having Δ 48 -excess. Samples with Δ 48 -excess in the second analytical session were not re-run. Those from the first analytical session were re-analysed after treating new aliquots of powdered calcite in a buffered solution of hydrogen peroxide, which eliminated excess Δ 48 . Zhang et al. (2020) report no effect on Δ 47 values within analytical precision for the peroxide treatment.
The statistical software packages R and RStudio (RStudio, 2019; R Core Team, 2020), and the IsoReader R package (Kopf et al., 2021), were used to calculate all ratios and perform all corrections directly from the raw data files. Isotope ratio calculations use the IUPAC 17 O parameters (Brand et al., 2010) and apply the 90°C acid fractionation factor of Kim et al. (2007). Raw Δ 47 values of individual cycles that were major outliers (>3 standard deviations) with respect to the mean of all cycles of that analysis were removed if the rest of the cycles demonstrated a normal distribution. The Δ 47 values were converted to ICDES (Intercarb Carbon Dioxide Equilibration Scale; Bernasconi et al., 2021) values using the ETH1 and ETH2 standards to correct for non-linearity effects ('compositional correction') and then corrected for compression of the Δ 47 scale ('scale correction') using ETH1-3 and IAEA-C2 standards in a two-step approach equivalent to the pooled regression approach of Daëron (2021). Values for all standards were from Bernasconi et al. (2021), and were checked via comparison of the resulting Δ 47 values for IAEA-C1 and ETH4. Analytical uncertainty and error associated with the reference frame corrections were fully propagated through each session and these propagated errors are the uncertainties reported for each individual analysis. Final reported Δ 47 values are the weighted average of the Δ 47 values of each replicate of the sample, using the propagated uncertainties for the replicate analyses as the weighting factor. Session 2 replicates with Δ 48 -excess (four from one sample and five from four others, File S2) were excluded from final averages.
Temperatures were calculated from the sample averaged Δ 47 values using the composite calibration of Anderson et al. (2021). Estimates of pore fluid δ 18 O values (‰, VSMOW) in isotopic equilibrium with calcite precipitated at those temperatures were calculated using those temperatures, the bulk calcite δ 18 O values, and the δ 18 O fractionation equation of Kim and O'Neil (1997). All results related to clumped isotopes are in File S2.

| Vitrinite reflectance
Constraints on peak burial temperatures (Tb peak ) were generated from pyrolysis-derived temperatures of maximum hydrocarbon generation (T max ) in the B and C chalks of the shallow Schock Errington #1 well (data accessed at https://my.usgs.gov/crcwc/ core/repor t/9383) and four of the proprietary wells using data generated by private service companies and provided by the wells' owners. The T max values were converted to vitrinite reflectance (%VR 0 ) values using the equation of Jarvie et al. (2001), and then %VR 0 to Tb peak with the equation of Pawlewicz (1986, 1994).

| 87 Sr/ 86 Sr and Nd isotope ratios
Strontium and neodymium isotope ratios were collected from B Chalk, C Chalk and Fort Hays calcite in the shallow Kansas well and in three of the proprietary wells in order to assess whether any fluids other than Late Turonian to earliest Campanian sea water might have been involved in Niobrara diagenesis. Strontium isotopes are a good geochronometer (McArthur et al., 2020) and a possible fluid indicator, as deviations from age-equivalent sea water 87 Sr/ 86 Sr can be a signal of extra-formational diagenetic fluids (Banner, 1995). Neodymium isotopes, which are not a good geochronometer due to the short residence time of Nd in sea water, can be a collaborating signal of extra-formational diagenetic fluids (Grousset et al., 1988).
Approximately 50 mg of carbonate from 12 samples were analysed on VG 54 and VG 354 sector thermal ionisation mass spectrometers at the University of Arizona. Samples were dissolved in weak acid to prevent dissolution of clays or silicates, and spiked after dissolution. Strontium data are presented as the 87 Sr/ 86 Sr ratio, with an average standard error of ±0.0010. Neodymium data are presented in standard εNd notation normalised to a 90 Myr depleted T A B L E 1 Stable oxygen isotope data averaged by well and stratigraphic unit.  mantle value, with standard errors averaging ±0.0017. The total dataset can be found in File S1, as well as more detail regarding the analytical methods.

| Oxygen, carbon and clumped isotopes
Stable oxygen isotope values measured independently of the clumped isotopes range from −4.5‰ to −12.4‰ VPDB ( Figure 8, Table 1). This range manifests as more pronounced 18 O-depletion in wells at greater depth and thermal maturity, as previously observed by Hefton (2015), Pollastro and Martinez (1985); Pollastro and Scholle (1986) and Scholle (1977). Except in the B Chalk of the shallowest wells (i.e. wells B-E; Carbon isotope values are comparatively uniform across all wells and stratigraphic units, with a majority of the total dataset forming a tight cluster between 1‰ and 2‰ VPDB (Figure 8). Similar values have been reported from the Niobrara Formation by Hefton (2015), Pollastro and Martinez (1985), Pollastro and Scholle (1986) and Scholle (1977). The Niobrara's positive δ 13 C values are also typical of Coniacian through lower Campanian deep-sea chalks (Clarke & Jenkyns, 1999;Huber et al., 2018;Jenkyns et al., 1995).
Stable oxygen and carbon isotope values measured concurrently with the clumped isotope data overlap with those of the bulk measurements ( Figure 8). Average Δ 47 values of individual samples ranged from 0.543 to 0.378‰, which correspond to temperatures of 44-145°C and parent fluid δ 18 O values of −1 to 8.3‰ VSMOW ( Figure 9; Table 2). These ranges dominantly manifest as increasing crystallisation temperature with increasing thermal maturity (Table 2), which also covaries with depth. No consistent stratigraphic patterns in temperatures occur within wells between the B, C and Fort Hays samples.

| Borehole estimates of palaeotemperatures
The thermal maturity spectrum for the Niobrara samples ranges from completely immature (0.38%VR o ) to wet gas mature (1.41%VR o ). That range translates to peak burial temperatures (Tb peak ) of 57°C in the shallow well from Kansas to 167°C in well N (Table 2). Variations within wells in Tb peak are minor, ranging from 0 to 4°C. Peak burial temperatures are greater than (i) clumped isotope temperatures by 10-69°C, (ii) current bottom hole temperatures (BHT) by 23-64°C, (iii) BHTs corrected by 18°C for the thermal damping effect of drilling fluids and (iv) projected maximum BHTs by 12-29°C (Table 2). The projected maximum BHTs assume additional burial (File S3) to account for erosion since the late Miocene (Figure 5), and a geothermal gradient calculated from current bottom hole depth and the BHT corrected for thermal damping effects. The BHT correction used the 'last resort' method of Corrigan (2003) because the wells' borehole log data were insufficient to perform more precise Horner or multi-point circulation corrections. the range of early Campanian sea water ( Figure 10A). That is, all ratios are greater than those expected for the Late Turonian through Santonian depositional ages of the Fort Hays, C and B chalks. Fort Hays samples do not exhibit lower 87 Sr/ 86 Sr ratios relative to C or B chalk samples, meaning there is also no stratigraphic pattern in 87 Sr/ 86 Sr ratios. The εNd values fall between −8.5 and −11.3 and covary with 87 Sr/ 86 Sr ratios ( Figure 10B).

THE ANALYTICAL RESULTS
The analytical results provide constraints on the thermal gradients and fluids associated with calcite diagenesis. Those constraints are considered in this section. Interpretation of what the analytical data indicates with respect to the origin of the unusually low δ 18 O values of the Niobrara chalks requires detailed diagenetic modelling based on the thermal and fluid constraints, which is presented in Section 6.

| Thermal history and gradients
Deciphering the Niobrara's thermal history is a requirement for understanding the chalks' δ 18 O values because the latter are dependent on temperatures at the time of the calcites' formation. Each well's thermal history is best reflected in the irreversible reactions of its organic material as those reactions will reflect maximum temperature. Thus, peak temperatures (Tb peak ) calculated from vitrinite reflectance measurements are the best estimates of maximum burial temperatures in the Niobrara Formation.
Geothermal gradients based on the Tb peak values and projected maximum burial depths range from 28 to 48.4°C/ km (Table 2), and are similar to modern geothermal gradients in the study area (30.8-49.9°C/km; Berkman & Watterson, 2010). The lowest gradient is for the Schock Errington well to the east (Figure 1). But there is still a variation of nearly 20°C/km in the five wells within the Wattenberg Field, an observation also noted by Thul and Sonnenberg (2018). The variability in those five wells is primarily due to the ca 35°C variations in Tb peak rather than the 30-200 m differences in maximum projected burial depths. Wells L, O and P also have Tb peak values that are 12-29°C higher than the projected maximum burial temperatures (Table 2). These observations combine to indicate local variability in heat flow and the occurrence of maximum heat flow at a time other than maximum burial.
Prior workers have attributed elevated palaeo heat flows in the Wattenberg Field to heat pulses along the Colorado Lineament that extends eastward into the Denver Basin from the CMB (Higley et al., 2003;Inks et al., 2019;Thul & Sonnenberg, 2018). Magmatism and mineralisation at the north-west end of the CMB trend is evidenced by ore deposits within 160 km or less of the study area. Those ores indicate a Laramide heat pulse at about 69-60 Ma and a mid-Cenozoic event at 37-24 Ma (Chapin, 2012).
The insulating effects of Late Cretaceous and Palaeocene coals over Wattenberg Field, which were eroded in the mid to late Cenozoic, is an alternative explanation for high palaeo heat fluxes (Coskey & Cumella, 2015;Landon et al., 2001). An estimate of those shallow insulators effect on peak temperatures and geothermal gradients illustrates their importance. An estimate of 800 m of eroded Cenozoic sediment is assumed, including ca 30 m of thermally insulating coal distributed in multiple beds over an 80 m interval (File S3). Also assumed is a surface temperature of 20°C (Pucéat et al., 2003); lithology-specific thermal conductivity values (Blackwell & Steele, 1989;Eppelbaum et al., 2014); a baseline thermal gradient of 60 mW/m 2 , which is typical east of the Front Range (Decker et al., 1988), and the burial history curve shown in Figure 5. With these assumptions, peak temperatures at the base of the Niobrara would have been ca 171°C (ca 46°C/km) at maximum burial. That value is compatible with the greatest Tb peak value of 161-167°C observed in the thermally most mature well N (Table 2). With shale in place of the coal, the peak temperature would have been only 153°C (40°C/km). The two calculations suggest that if the coal were laterally variable in thickness and magnitude, then it could explain part, but not all of the 20°C/km variability in geothermal gradient across the studied wells in Wattenberg Field. A transient heat pulse related to the CMB would be implicated to explain the balance. Thus, a range of geothermal gradients and thermal events are considered in the diagenetic modelling presented in Section 6.  Kim and O'Neil (1997

| Interpretation of clumped isotope palaeotemperatures
Given the nano and microcrystalline composition of the Niobrara chalks (Figure 7), the calcite δ 18 O values and Δ 47 temperatures are the weighted averages of all calcite formed during a sample's history. The weighted average could include unaltered biogenic calcite formed in the Cretaceous seaway's water column, any calcite cement precipitated concurrently with mechanical compaction over the first hundreds of metres of burial, all calcite associated with recrystallisation, and calcite cement produced in association with pressure solution. All but the latter should dominate in the shallow buried Kansas samples that show little evidence for pressure-solution cementation ( Figure 7A). In contrast, weighted average δ 18 O values and Δ 47 temperatures in the more deeply buried and highly altered samples from Wattenberg Field ( Figure 7B consideration that can complicate the interpretation of the Δ 47 temperatures. Heating during burial can more randomly distribute the heavy isotopes, which produces artificially 'hot' Δ 47 temperatures, and re-'clumping' the heavy isotopes during exhumation and cooling can produce artificially 'cold' Δ 47 temperatures. Reordering is a kinetic phenomenon that is initiated during burial at a threshold temperate of approximately 130°C (Henkes et al., 2014;Stolper & Eiler, 2015) with an upper temperature boundary at which the primary ordering state is completely erased (Henkes et al., 2014;Passey & Henkes, 2012). The upper boundary depends on the rate of reordering, which increases with both temperature and the amount of time that the calcite is exposed to the elevated temperature. For example, Henkes et al.'s (2014) model for the reordering of brachiopod calcite (microcrystalline calcite) reveals the onset of reordering at 133°C with a rate of 1% reordering per 10 6 years, and the rate accelerating at 150°C to 1% reordering per 10 5 years. Hemingway and Henkes' (2021) disordered kinetic model for reordering predicts even less time and/or lower temperatures for reordering.
Reordering in the Niobrara samples from the Schock Errington, N and L wells is rejected because those wells have peak temperatures less than or no greater than the threshold temperature for reordering. However, partial to complete reordering may have occurred in some well O and P samples as those wells have weighted average peak temperatures of 144-163°C respectively ( Table 2). The burial history ( Figure 5) and observed geothermal gradients of these two wells also suggest that the elevated temperatures would have lasted for millions of years in well O and probably tens of millions of years in well P. Eight of 10 replicates from well P and two of 16 from well O have average Δ 47 temperatures greater than 125°C (Figure 9). One well P replicate (Δ 47 temperature of 169°C) and the two well O replicates (Δ 47 temperatures of 127 and 129°C) are also 20-25°C outliers relative to the rest of the weighted average Δ 47 temperatures in those wells. These three replicates, and possibly all of the other well P replicates, are interpreted to be the most probable samples influenced by reordering. This interpretation is further evaluated by the diagenetic modelling (Section 6).

| Rock-buffered carbon system
The positive δ 13 C values of the Niobrara chalks are interpreted to reflect the carbon isotopic compositions of the original carbonate sediment, with little if any diagenetic overprint. That in turn means a closed, rock-buffered system with respect to carbon throughout the Niobrara's diagenetic history. Authigenic calcite formed during recrystallisation inherited the δ 13 C value of the dissolving calcite with no increasingly lower δ 13 C values with increasing depth and temperature because carbon isotope fractionation between aqueous HCO 3 − and calcite is temperature independent (El-Shenawy et al., 2020; Levitt et al., 2018;Romanek et al., 1992).

| Constraining diagenetic fluid sources
The Sr isotope ratios of Niobrara calcite ( Figure 10A) indicate that Sr was added to the chalks during diagenesis. The additional Sr had an isotopic ratio greater than the original Niobrara calcite, but how much greater is unknown. The Sr isotopic ratio of the diagenetic fluid might fingerprint the origin of that fluid (Banner, 1995), but that ratio cannot be constrained with the data available. The covariance in εNd and 87 Sr/ 86 Sr do suggest a common source of both. Two fluid sources are plausible given the stratigraphy and burial history of the Niobrara. The first possibility is that Sr and Nd were derived from the smectitic clays within the Niobrara marls and/ or the shales and marls in the underlying late Cretaceous sedimentary column (i.e. Graneros and Greenhorn formations; Figure 10A). Smectite diagenesis is documented in the Cretaceous section (Elliott et al., 1991;Pollastro, 1990), and smectite's dehydration from three to two to one layers of structural water, then conversion to illite, progresses with increasing burial temperature and pressure (Tremosa et al., 2021;Vidal & Dubacq, 2009). Smectite water, with a δ 18 O value of +20‰ SMOW (Savin & Lee, 1988;Sheppard & Gilg, 1996;Tremosa et al., 2021), would have altered the ambient connate sea waters' δ 18 O value, with the mixed water then fluxed upward as the basin compacted and dewatered, particularly during the first rapid burial segment prior to Laramide deformation ( Figure 5).
The second possibility is that Palaeocene and Eocene meteoric groundwaters infiltrated into the basin beginning with the Laramide orogeny about 67.5 Ma (burial segment 2 and 3; Figure 5; Jorgensen, 1989). The meteoric recharge primarily entered Cretaceous (i.e. Dakota Group sandstones below the Niobrara Formation) and older Palaeozoic sandstone aquifers and flowed into the basin (Figure 11). The topographically driven flow F I G U R E 1 1 Cross-section through the Denver Basin depicting generalised distribution of stratigraphic units at the end of the Cretaceous, with the uplifted mountain front created by the Laramide orogeny shown schematically to the west. Golden fault marks the boundary of the faulted and uplifted rocks with the Denver Basin. Datum placed on top of the Pierre shale because actual palaeo elevation of land surface is unknown. Meteoric groundwater (blue arrows) was topographically driven from the mountains into the basin (Jorgensen, 1989). Faults are shown schematically and based on Weimer (1996) and Sonnenberg et al. (2016). Cross-section modified after Umari et al. (2018). See Figure 2C for the location of the line of section. system was probably artesian (over pressured) and that in turn drove the meteoric groundwater upwards into the Niobrara chalks along conduits, such as fractures, the deep-seated wrench faults of Weimer (1996), and the shallower polygonal fault system of Sonnenberg et al. (2016), Radiogenic Sr and Nd would have been derived from rock-water interaction with siliciclastic minerals along the flow path. This flow system would have existed until either gas generation in the mid to late Cenozoic (Higley & Cox, 2007;Landon et al., 2001) created over pressure in the Niobrara chalks and underlying mid Cretaceous rocks (Spencer, 1987;Swarbrick & Osborne, 1998;Swarbrick et al., 2002;Weimer, 1996), or tilting and erosion in the last 27 Myr created the current under pressure in the aquifer rocks (Belitz & Bredehoeft, 1988;Umari et al., 2018). This meteoric flux is known to have driven post-oil generation burial diagenesis in Lyons (Permian) oil wells in the basin axis (Lee & Bethke, 1994).

| The isotope model
Equilibrium isotope models simulate the δ 18 O values of calcite and pore fluid during burial and diagenesis based on isotopic mass balances and temperature-dependent isotopic fractionation between fluid and authigenic calcite. The analytical results discussed above help to condition these models in terms of fluids and geothermal gradients. Because all δ 13 C values are rock-buffered, carbon isotopes are not included in the model.
The isotopic evolution of the chalk is modelled to occur by the simultaneous dissolution and reprecipitation (i.e. recrystallisation) of small increments of the total mass of calcite in the chalk. Only primary depositional calcite is assumed to be reactive, and once reacted that mass of calcite does not react again. Bulk δ 18 O values and temperature derived at each iteration are mass-weighted averages of the remaining primary calcite and all the increments of secondary calcite formed at each shallower depth increment.
The biggest challenge of this approach is constraining the incremental rate of reaction (Maliva et al., 1991). Pore-water chemistry provides insight when modelling shallowly buried and young deep-sea chalks (Schrag et al., 1995;Stolper et al., 2018), but in petroliferous basins the pore waters associated with much of the chalks' diagenetic history are no longer present. In such settings, linear incremental reaction rates have been assumed (Egeberg & Saigal, 1991;Maliva et al., 1991), and that is the approach taken here. All simulations assume that porosity loss over the first 600 m of burial was due to mechanical compaction with no recrystallisation or significant calcite cementation. This value was chosen based on Pollastro and Martinez's (1985) and Pollastro and Scholle (1986) statements that mechanical compaction dominated Niobrara porosity loss to about 40%, which occurs at 600 m in the porosity-with-depth curve explained below. Calcite diagenesis is simulated at 10 m intervals from 610 m to maximum burial at 3,250 m, so the linear incremental rate is the mass of original calcite in the system divided by 265 increments. This equates to a recrystallisation rate of about 1.4% per myr.
The total mass of calcite is conserved, meaning a closed system with respect to calcite, but pore volume decreases with porosity in accordance with a Niobrara-specific porosity curve that is based on measured core-plug porosity ( Figure 12). The curve is a non-linear parametric transformation of Schmoker and Halley's (1982) limestone porosity-with-depth equation: where ϕ i is the initial porosity at zero burial depth, −a is an exponential curvature factor, depth' the depth over which porosity decreases by approximately 1/e, and θ the empirically derived rotation factor. The fitted curve yields an initial porosity of 57% at zero depth, a reasonable value for a chalk ooze within many tens of metres of the sea Porosity = ϕ i * ê( − a * depth)∕depth � * cos( ) − depth * sin( ) F I G U R E 1 2 Non-linear porosity model (solid line) fitted to measured Niobrara core plug porosities. The Bayesian Information Criteria (Schwarz, 1978) indicates that the non-linear model, which is a parametric transformation of the depth equation given by Schmoker and Halley (1982), is a more robust representation than a linear model. Depths for shallow wells towards the eastern margin of the basin are plotted at sample depth plus 400 m as a minimum correction for sediment eroded in the last 7 Ma (Higley & Cox, 2007); depths for deeper wells in the basin axis are plotted at sample depth plus 800 m to account for even greater eroded material (File S3).
The isotopic mass balance and fractionation equations, in the sequence in which they are computed, are: where n refers to the current iteration step, n−1 the prior step, MVO f and RMVO c the molar volumes of oxygen in the fluid and reacting increment of calcite, respectively; α the temperature-dependent fractionation factor for calcite-water from Kim and O'Neil (1997); and δ 18 O f , δ 18 O c and δ 18 O sys the oxygen isotopic values of the fluid, calcite and total diagenetic system respectively. The local mass transfer of calcite via pressure solution is treated as an instantaneous process that generates no change in the total molar volume of calcite. All simulations were done with the open-source computing program R (R Core Team, 2020).
Temperature (which affects the fractionation factor, α) is determined at each depth increment by the geothermal gradients established by Tb peak values. Constant gradients of 30 and 40°C per km, which span all wells but one (Table 2), were utilised. They produce maximum temperatures of 115.5 and 148°C, respectively, at maximum burial. A variable heat flux was simulated with an initial geothermal gradient of 35°C/km from 0 to 1,600 m of burial, then a gradient of 65°C/km from 1,600 to 2,300 m that accounts for a Laramide (69-62 Ma; Figure 5) heat pulse (Chapin, 2012), followed by a gradient of 44°C/km from 2,300 m to maximum burial to account for the insulating effect of Late Cretaceous coals. This varied gradient produces a maximum temperature of 167°C at the end of the thermal pulse. Sensitivity tests showed that a mid-Cenozoic heat pulse near maximum burial resulted in no appreciable effect on model outcomes due to low porosity and fluid: rock ratios by that point in the burial history. The sequence of calculations assumes porosity loss before isotopic re-equilibration due to temperature increase. The alternative, calculating temperature increase and isotopic re-equilibration followed by porosity loss, yields a difference in δ 18 O values of just 0.1‰, which is insignificant.
Other inputs to the model are δ 18 O of sea water and sea water temperature, with initial calcite δ 18 O derived from those variables. Simulations were done with sea water isotope ratios of −1, −3 and − 4‰ VSMOW. The −1 value represents global Cretaceous sea water (Lhomme & Clarke, 2005;Shackleton & Kennett, 1975), which appears to have dominated the Western Interior Seaway (Jones et al., 2022). The −3 and − 4‰ values represent a freshening of seaway waters due to meteoric influx as suggested by Arthur et al. (1985) and Pagani and Arthur (1998). All simulations used a sea floor temperature of 25°C, which is the mean sea floor palaeotemperature (relative to the Anderson et al., 2021, Δ 47 -temperature calibration) extrapolated from the data of Jones et al. (2022, their figure 4) during Niobrara C and B chalk deposition (i.e. about 85-86 Ma; Sageman et al., 2014). For simplicity, it is assumed sea surface temperature equals sea floor temperature and the initial δ 18 O value of the pelagic biogenic (depositional) calcite is then calculated from the initial sea water's δ 18 O value at 25°C using the fractionation equation of Kim and O'Neil (1997). A vital effect during calcification by nanoplankton (Tagliavento et al., 2019) was ignored in the calculation of the initial calcite δ 18 O value because that effect is minimised under high pCO 2 conditions and considered minimal prior to the Eocene (Bolton et al., 2012;Hermoso et al., 2016).

| Modelling results
6.2.1 | Closed hydrologic system Figure 13 shows simulation results for a diagenetic system closed to external fluid throughout the Niobrara's entire burial history. The modelled outcomes for weighted average calcite δ 18 O values and Δ 47 temperatures only bracket the measured data for the shallow-buried Schock Errington B and some well L, N and O samples, but only for initial sea water δ 18 O values of −3 to −4‰ (SMOW). None of the simulations yield calcite δ 18 O values less than −9.5‰. The trajectory of all simulations is sharply concave upward and to the right with calcite δ 18 O values nearly unchanging at temperatures greater than 80°C. In contrast, a hypothetical line through the measured data would be broadly concave downward and to the left with calcite δ 18 O values continuing to evolve as temperature went up.
These simulations indicate that a hydrologically closed fluid system that just contains the original Cretaceous sea water cannot explain the δ 18 O evolution of Niobrara chalks. There simply is not enough 16 O in the original carbonate sediments and marine pore fluid to yield the observed δ 18 O values in the more diagenetically altered and deeply buried Niobrara chalks.

| Semi-open, leaky hydrologic system
A semi-open, 'leaky' hydrologic system is required to successfully model the measured data. In the closed system, continued recrystallisation results in the pore fluids' becoming progressively enriched in 18 O (e.g. +1, +3, +5, +7 VSMOW, etc.). In a semi-open system, the incoming fluid must mitigate that large increase in the pore water's δ 18 O value by supplying new 16 O, which happens as long as the new water has a δ 18 O value less than the existing water. Two possibilities raised in Section 5.3 are compatible with the history of the basin. First is compaction-driven marine pore fluid associated with burial segment 1 ( Figure 5). This fluid would have been isotopically evolved from the initial Cretaceous sea water δ 18 O value by the dehydration of smectite. Second is a meteoric fluid introduced during burial segments 2 and 3. In both cases, a small but continuous leakage of the extra-formational fluid is considered. The leakage is simulated over a range of depths by replacing a small percentage of the fluid volume present at the end of each depth iteration with the extra-formational fluid. The temperature of the two fluids is assumed to be equal at each depth iteration.
The results ( Figure 14A) show that this extraformational fluid does result in lower calcite δ 18 O values and higher temperatures at a given depth step relative to the closed system models, yet the models bracket very little of the measured data aside from the shallowest buried well to the east (Schock Errington). The trajectory of each simulation is again concave upward rather than downward because the system is hydrologically closed after 2,250 m of burial. If the models remain open at greater depths to the 0.8‰ compactional fluid, then the simulated trends flatten on Figure 14A (not shown) and pass below all of the measured data. In either case, there is still not enough new 16 O fluxing into the system at higher temperatures and burial depths to yield both the observed weighted average bulk calcite δ 18 O values and temperatures in the more deeply buried wells.

Meteoric water influx
The results of the prior two sets of simulations indicate that a very significant source of 16 O that persists into deep burial is needed to explain the lowest weighted average bulk calcite δ 18 O values and generate the concave downward trajectory of the measured data in calcite δ 18 Otemperature space. As noted in Section 5.4, a source of such a fluid is meteoric groundwater that moved downward from orthographic recharge zones on the western margin of the basin, moved laterally through the Dakota aquifer, δ and then moved upwards along conduits into the overlying Niobrara chalks ( Figure 11). The meteoric groundwaters are assumed to have had initial δ 18 O values of about −9.5‰ SMOW based on unpublished bioapatite analyses from the basal Palaeocene in the Denver Basin (see Acknowledgements). Rock-water alterations would have caused those waters to evolve to greater δ 18 O and rockbuffered δ 13 C values during their transit to the Niobrara chalks. Accordingly, a meteoric fluid δ 18 O value of −6.5‰ SMOW was used. Sensitivity tests indicate a leakage of 2.5% of that pore fluid yields optimum results. The influx was initiated at a burial of 2,250 m, which corresponds to when topographically driven post-Laramide meteoric groundwater began to enter the basin (Jorgensen, 1989). The influx was continued through maximum burial. Prior to 2,250 m, the hydrologic system was assumed closed.
The results of these meteoric influx simulations ( Figure 14B) show reasonably good matches with the measured well data. In particular, they have a concave downward trajectory. Best fits to the weighted average measured data from the deeply buried basin-axis wells occur with original sea waters of −1 ± 1‰ (VSMOW) and varied thermal gradient (dashed lines in Figure 14B). The simulation results of those scenarios pass through most of the data cloud related to wells L, N, O and P and produce weighted average calcite δ 18 O values as low as −11.5‰ at a weighted average temperature of 117°C. The Niobrara calcites with the least negative weighted average δ 18 O values (−5 to −6‰ VPDB) and lowest weighted average temperatures (i.e. the Schock Errington data) are still successfully modelled as having formed from −2 to −3‰ (VSMOW) Cretaceous sea water in the closed system before the meteoric influx. The fact that the Schock Errington witnessed a more negative sea water value than the other wells could reflect its location closer to the eastern shoreline and thus freshwater dilution of the sea water.

Composite extra-formation fluid influxes
A final set of simulations evaluate a semi-open system that includes both the early (1,000 to 2,250 m burial depths) compaction-driven and isotopically evolved sea water (δ 18 O of 0.8‰ VSMOW), the varied geothermal F I G U R E 1 4 Evolution of weighted average Niobrara calcite δ 18 O and temperature values (solid and dashed lines) for semi-open hydrologic systems. Solid and dashed lines, initial sea water δ 18 O values and geothermal gradients as in Figure 13. Measured data and their weighted averages from Figure 9. (A) Original sea water is mixed with smectite dehydration water to generate an isotopically evolved sea water with δ 18 O value of 0.8‰ (VSMOW) that leaks (10% by volume per depth iteration) into the chalks from 1,000 to 2,250 m as the section is compacted and dewatered. Thereafter, the hydrologic system is closed. (B) an initially hydrologically closed system consisting of the initial sea water is succeeded by a semi-open system in which a meteoric fluid (−6.5‰ VSMOW) leaks (2.5% by volume per depth iteration) into the chalks. Leakage is simulated to begin after the Laramide orogeny (burial depth of 2,250 m; Figure 5) and continue until diagenesis ended at maximum burial.
δ B A gradients due to heat pulses and the insulating effects of now eroded coals, and late meteoric leakage (2,250 m to maximum burial). The first of these scenarios utilises a 2% leakage rate for the compactional fluid and a 2.5% leakage rate for the −6.5‰ (VSMOW) meteoric groundwater. The results produce weighted average calcite δ 18 O values as low as −13.5‰, but at cooler weighted average temperatures than observed in the measured data ( Figure 15A). The models' trends are also linear, not concave downward, and pass below the cloud of measured data, indicating that the amount of admixed smectite dehydration water is too great.
A second set of composite models use the same varied geothermal gradients and the same meteoric influx parameters, but only a 0.8% leakage rate for the compactional fluid. Depending on specific temperature/depth and initial sea water δ 18 O value, the results of this scenario ( Figure 15B) produces weighted average calcite δ 18 O values that are just 0.3 to 0.9‰ lower than the simulation with late meteoric influx but no early compaction-driven flow ( Figure 14B). The modelled curves with initial sea water δ 18 O values of 0 to −2‰ VSMOW bracket most of the measured data and are concave downward in calcite δ 18 O-temperature space. A small flux of compaction fluids (1% or less by pore volume) are thus tenable when a meteoric flux is also included, yet it is the dominance of the meteoric flux that is needed to reproduce the concave downward trajectory of the measured data.

| Origin of anomalously negative δ 18 O values in Niobrara chalks
The clumped isotopic data and isotopic modelling of that data resolves the question of why Niobrara chalks have significantly more negative δ 18 O values than other late Cretaceous chalks (Figure 1). The uniqueness of the Niobrara's calcite δ 18 O values is primarily the result of two factors unique to the Denver Basin. Foremost is the influence of 18 O-depleted meteoric groundwaters introduced F I G U R E 1 5 Composite models of the evolution of weighted average Niobrara calcite δ 18 O and temperature values. Solid and dashed lines, initial sea water δ 18 O values and geothermal gradients as in Figure 13. Measured data and their weighted averages from Figure 9. (A) an initially hydrologically closed system consisting of the initial sea water is succeeded by a semi-open system in which the original sea water mixes with smectite dehydration water to generate an isotopically evolved sea water with δ 18 O value of 0.8‰ (VSMOW) that leaks (2% by volume per depth iteration) into the Niobrara chalks from 1,000 to 2,250 m. thereafter an isotopically evolved meteoric fluid (−6.5‰ VSMOW) leaks into the Niobrara (2.5% by volume per depth iteration) until maximum burial. (B) Same fluid types but with only a 0.8% leakage of the original sea water mixed with smectite dehydration water. δ A B into the basin after the Laramide orogeny. Of lesser importance is the pulse of high heat flow associated with magnetism along the CMB during the Laramide orogeny. The absence of very negative δ 18 O values in other Cretaceous chalks presumably reflects the absence of leaky meteoric recharge deep into those rocks while pressure solution and calcite recrystallisation were still ongoing. Late calcite cements may impact calcite δ 18 O values in those other settings (Maliva et al., 1995) but not the overall bulk calcite δ 18 O values of the chalks.
With the exception of the shallow eastern margin of the Denver Basin, the diagenetic models do not indicate that sea water with a δ 18 O value of −3‰ SMOW, as suggested by Arthur et al. (1985), Pollastro and Scholle (1986), and Pagani and Arthur (1998), was the cause of the Niobrara's very negative calcite δ 18 O values. Rather, the modelling indicates that the open seaway waters had the same isotopic value (−1‰, VSMOW) as the global Cretaceous ocean; a conclusion in agreement with recent clumped isotope studies of unaltered calcitic molluscs (Jones et al., 2022;Petersen et al., 2016). For the samples from the basin's axis, the ±1‰ variability between the measured data and the final composite model invoking −1‰ sea water ( Figure 15B) is probably not the result of variation in the seaway's oxygen isotope value. Rather the variability is probably due to differences between wells in terms of actual geothermal gradients (Table 2; Section 5.2) and the specific volumes of early compaction-driven and later meteoric groundwater that interacted with each sample during burial segments 2 and 3. Petersen et al. (2016) did find evidence of freshwater dilution and sea water δ 18 O values of −3 to −4‰ VSMOW along the eastern margin of the seaway, and the isotopic modelling indicates that was the case for the Schock Errington chalks on the eastern margin of the Denver Basin. Closed system alteration in sea waters through shallow burial in sea water with an initial δ 18 O value of −2 to −3‰ VSMOW reproduces the measured Δ 47 temperatures and calcite δ 18 O values in those chalks. Their alteration ceased prior to, or was unaffected by, the Laramide heat pulse and the meteoric influx. Aside from the lower initial sea water δ 18 O values, the modelling indicates that all shallow-buried Niobrara chalks on the eastern margin of the basin (Table 1) had an alteration history similar to other Cretaceous chalk deposits.
The Laramide heat pulse is also modelled to have generated a peak temperature of 167°C for some millions of years. That is sufficient heat and time above 150°C to partially reorder 13 C-18 O bonds (Hemingway & Henkes, 2021;Henkes et al., 2014). This supports the inference in Section 5.2 that the elevated clumped isotope temperatures in some well P and O samples reflect partial reordering. It would also account for the failure of all the simulations (Figures 13-15) from matching weighted average Δ 47 temperatures in the P well even when the most negative calcite δ 18 O values were successfully simulated.
The diagenetic history revealed by the data and modelling presented in this study, and especially the fluid history associated with diagenesis, is far more complex than that envisioned by earlier studies of Niobrara chalks that only possessed conventional stable isotope data (Hefton, 2015;Pollastro & Martinez, 1985;Pollastro & Scholle, 1986;Scholle, 1977). Like other recent efforts on chalk diagenesis that have utilised clumped isotopes (Smit et al., 2022;Tagliavento et al., 2021), the results of this study highlight the need to improve the 'chalk burial diagenesis model' with the use of clumped isotope constraints on diagenetic temperatures.

| Model sensitivity and limitations
There are seven major input parameters to the simulations that control the temperatures and calcite δ 18 O values at maximum burial. These are: the initial fluid δ 18 O value and temperature; the thermal gradient; the fluid systems behaviour (closed or semi-open and timing of changes from one to the other); the depth at which diagenetic alteration begins; fluid: rock ratio (controlled by pore volume at a constant reaction rate); and the influxing fluid's δ 18 O value and volume. The interaction of these variables is complex (Figure 16), in part because the δ 18 O value and F I G U R E 1 6 Relationship between the seven model input parameters and the general trend of model outcomes. Key to the right includes the sensitivity of key model outcomes to those input parameters. Specific outcomes are initial temperature (T i ), weighted average temperature at maximum burial (T f ), indices or curvature (k 1 , k 2 ), initial calcite (δ 18 O i ), and the final weighted average calcite δ 18 O. See text for discussion of line segments a-d. crystallisation temperature of the calcite are weighted average values.
The weighted average temperature at any burial step is primarily controlled by the thermal gradient, which is why a range of thermal gradients were modelled. A secondary control is initial temperature-any change in that value results in an equal change in final weighted average temperature at maximum burial. Initial temperature on the Cretaceous sea floor, however, is reasonably well constrained by Jones et al. (2022) with just a few degrees of uncertainty.
The depth at which calcite diagenesis begins is also a secondary control on the weighted average temperature. If calcite recrystallisation begins at a depth less than the 600 m modelled, the weighted averages temperatures at any one depth are lower whereas weighted average δ 18 O calcite values are greater ( Figure 16, T f ). Conversely, if calcite recrystallisation begins below 600 m, the weighted average temperature is higher and weighted average δ 18 O calcite values are lower. The effect, however, is small when surface and bottom waters have the same temperature (Stolper et al., 2018), as assumed herein and in other chalk recrystallisation models. For the geothermal gradients used herein, the effect at maximum burial is to only change δ 18 O values 0.2-0.3‰ and temperature by 2-3°C per 150 m difference in the onset of diagenesis.
The fluid system is the dominant control on the weighted average calcite δ 18 O values. Because each iteration of the model reacts the same volume of initial (depositional) calcite, it is the changes in the fluid volume (i.e. pore volume from Figure 12) and fluid δ 18 O values that have the most critical influence on the total oxygen isotope system at any given iteration step. As a result, in simulations with a significant difference between the δ 18 O values of calcite and fluid, particularly when pore volume is still large, the fluid effects on fractionation are more impactful than the temperature effect. This is why the modelled trend lines in δ 18 O-temperature space were concave downward only with the meteoric fluid. For fluids with positive δ 18 O, the fluid impact on fractionation can even result in increasingly higher δ 18 O calcite values with burial, which is a direct reversal of the expected isotopic evolution.
Fluid system effects are visualised in Figure 16 relative to the curvature of the calcite δ 18 O-temperature trend. A closed system simulation creates a concave-up curve (Figure 16, dashed line a), the intensity of which (κ 1 ) is controlled by the depth at which calcite diagenesis begins. In contrast, a semi-open fluid system results in a flat (low κ) to concave-down (high κ 2 ) curve, with the orientation relative to the δ 18 O axis (dashed lines b and c, Figure 16) controlled by the δ 18 O values of the influxing fluid. The thermal gradient and fluid: rock ratio (pore volume) control the height of the curvature on the modelled trend (Figure 16, box d), but because fluid system effects are more significant, it can be difficult to visually detect temperature-controlled inflection points in more complex simulations ( Figure 15B vs. Figure 14B).
In terms of specific input parameters, the greatest potential uncertainties in the isotopic models occur with respect to the parameters that affect the δ 18 O value of the initial calcite, plus the uncertainty in the influx rate and δ 18 O values of the various fluids modelled as 'leaking' into the Niobrara during burial diagenesis. Each is considered below.
The initial pore fluid (sea water) δ 18 O value has a significant impact on model outcomes, which is why a range of values were considered. It had greater impact than thermal gradients, as evidenced by Figures 13-15, where different thermal gradients for a specific initial sea water δ 18 O value are more tightly spaced with respect to each other than simulations of different initial sea water δ 18 O values at any one thermal gradient.
The initial sediment δ 18 O value, derived from the initial sea water's δ 18 O value and temperature, also impacts the model results because the calcite contains a high percentage of the total molar volume of oxygen in the initial diagenetic system. Again, the range in initial sea water δ 18 O values accounts for this uncertainty. No variation in initial sea water temperature was considered due to the strong evidence for a 25°C sea water temperature (Jones et al., 2022). However, another factor that could affect uncertainty in the initial calcite δ 18 O value is the potential role of a vital effect on fractionation during coccolith precipitation by nanoplankton. No vital effect was considered because they are poorly constrained and considered minimal during the Cretaceous (Bolton et al., 2012;Hermoso et al., 2016). Nonetheless, Tagliavento et al. (2019) argued that the vital effect with respect to bulk Latest Cretaceous chalks could be as low as −2‰, whereas the observed bulk δ 18 O value of shallowly buried and pristine appearing Maastrichtian and Campanian deep-sea coccolith debris is in equilibrium with non-biogenic calcite if the bulk vital effect is between +1 and +1.5‰ (Berrera & Savin, 1999;Huber et al., 2002). Adding a negative vital effect reduces but does not eliminate the need for late 18 O-depleted fluids. A positive vital effect requires either significantly more 18 O-depleted fluids or limits the successful simulations to an initial sea water value of −3‰ rather than −1‰ (VSMOW).
Lastly, the δ 18 O values and influx rates of extraformational fluids in the hypothesised semi-open hydrologic system are poorly constrained, but they do identify the fluid δ 18 O values and volumes that are required to match the observed clumped isotope data. The bigger uncertainty is whether 18 O-depleted fluids such as meteoric groundwater circulated into the basin and reached the Niobrara chalks after the Laramide orogeny. Addressing that uncertainty is beyond the scope of this work.

| CONCLUSIONS
Multiple hypotheses have been proposed to explain the anomalously low δ 18 O values of Niobrara chalks relative to other Cretaceous chalks. New clumped isotope data coupled with numerical models simulating calcite recrystallisation indicate that the freshening of the Cretaceous sea water in the Western Interior Seaway is not needed to explain the very low calcite δ 18 O values. Rather those values are primarily the result of meteoric recharge deep into the Denver Basin after the Laramide orogeny. This created a semi-open hydrologic system that steadily 'leaked' 18 O-depleted (meteoric) groundwater into the Niobrara chalks through maximum burial. The influx of these 18 O-depleted fluids kept pore fluid δ 18 O values from significantly increasing as calcite recrystallisation continued with burial, which in turn allowed calcite δ 18 O values as low as −11.5‰ to be achieved. A basement heat pulse, the thermal insulating effects of now eroded coals, and a small flux of compaction-driven Cretaceous sea water that was isotopically evolved by smectite dehydration all contributed to, but are not the main drivers of the anomalously low calcite δ 18 O values. This work confirms that the broad diagenetic processes reflected by δ 18 O values of the Niobrara chalks are comparable to other chalks worldwide. However, the detailed insights derived from clumped isotope analyses highlight the attributes of the diagenetic fluids and fluxes that made the δ 18 O record of Niobrara chalk unique with respect to other Cretaceous chalks.