Cretaceous oceanic anoxic event 2 (OAE2) is thought to have been contemporary with extensive volcanism and the release of large quantities of volcanic CO2 capable of triggering marine anoxia through a series of biogeochemical feedbacks. High-resolution reconstructions of atmospheric CO2 concentrations across the initiation of OAE2 suggest that there were also two distinct pulses of CO2 drawdown coeval with increased organic carbon burial. These fluctuations in CO2 likely led to significant climatic changes, including fluctuations in temperatures and the hydrological cycle. Paleofire proxy records suggest that wildfire was a common occurrence throughout the Cretaceous Period, likely fueled by the estimated high atmospheric O2 concentrations at this time. However, over geological time scales, the likelihood and behavior of fire are also controlled by other factors such as climate, implying that CO2-driven climate changes should also be observable in the fossil charcoal record.

We tested this hypothesis and present a high-resolution study of fire history through the use of fossil charcoal abundances across the OAE2 onset, and we compared our records to the estimated CO2 fluctuations published from the same study sites. Our study illustrates that inferred wildfire activity appears to relate to changes in CO2 occurring across the onset of OAE2, where periods of CO2 drawdown may have enabled an increase in fire activity through suppression of the hydrological cycle. Our study provides further insight into the relationships between rapid changes in the carbon cycle, climate, and wildfire activity, illustrating that CO2 and climate changes related to inferred wildfire activity can be detected despite the estimated high Cretaceous atmospheric O2 concentrations.

Periods of oceanic anoxia have occurred many times in Earth’s history, and they were particularly common throughout the Cretaceous Period (Jenkyns, 2010). These anoxic events are evidenced in the geological record by common occurrences of widespread black shales and excursions in the organic and carbonate δ13C record, indicating significant disruption to the global carbon (C) cycle (Jarvis et al., 2011).

Cretaceous oceanic anoxic event 2 (OAE2; Schlanger and Jenkyns, 1976) occurred ca. 94.45 Ma, spanning the Cenomanian-Turonian boundary (ca. 93.90 Ma; Meyers et al., 2012). OAE2 is hypothesized to have been triggered by an influx of volcanically derived CO2, driving increases in global temperatures, enhancing the hydrological cycle, and increasing continental weathering rates (Du Vivier et al., 2015; Jenkyns et al., 2017). These changes are hypothesized to have preconditioned the oceans in the runup to the event, creating low-oxygen conditions. At the onset of OAE2, these changes, coupled with rising sea level, are hypothesized to have enabled a flux of otherwise limiting nutrients to the marine system, enhancing primary production and Corg burial, and initiating the event (Mort et al., 2007; Owens et al., 2016). During the initiation, enhanced preservation of organic matter is hypothesized to have taken place in euxinic hotspots (Owens et al., 2016), which constituted a small proportion of the global seafloor area (<10%; Dickson et al., 2016).

Support for this scenario comes from a high-resolution stomatal pCO2 proxy analysis, covering the runup to and/or onset of OAE2 (Barclay et al., 2010). Using the stomatal proxy, Barclay et al. (2010) suggested that CO2 levels began to gradually rise prior to the start of OAE2, reaching a peak of 500 +400/–180 ppm. However, this overall rise in pCO2 was interrupted by two inferred rapid phases of CO2 drawdown, taking just ∼10–100 k.y. (Barclay et al., 2010), which also appear to be expressed in the global δ13C record as small positive excursions. These rapid changes in pCO2 are hypothesized to have led to significant climate changes, causing large shifts in Northern Hemisphere temperatures and the hydrological cycle (Jenkyns et al., 2017). By contrast, in the Southern Hemisphere, recent TEX86 data (where TEX86 is an index of tetraethers consisting of 86 carbon atoms) suggest a period of sustained warmth throughout the OAE2 interval, with little climate fluctuation (Robinson et al., 2019).

Within the Earth system, combustion is an important process (Li et al., 2014) that influences global ecosystems, climate, and the C cycle by generating feedbacks that influence Earth’s biogeochemical cycles (e.g., Kump, 1988; Lenton and Watson, 2000; Handoh and Lenton, 2003; Bowman et al., 2009). A key example of this can be seen in fire’s role in preventing atmospheric O2 levels from rising too high, and thus maintaining concentrations within habitable bounds. Variations in the occurrence and abundance of fossil charcoal, which is a by-product of paleowildfire that is used as a proxy for changes in fire activity, have been shown to fluctuate in accordance with changes in climate (Daniau et al., 2007, 2013), as well as with broad changes in atmospheric O2 (Belcher and McElwain, 2008; Glasspool and Scott, 2010; Baker et al., 2017).

For the Cretaceous Period, biogeochemical models estimate high atmospheric O2 concentrations between 23 vol% and 27 vol% (Lenton, 2013; Mills et al., 2016; Lenton et al., 2017). Laboratory experiments show that given a source of ignition, the probabilities of ignition and fire spread beyond 23 vol% are high, and they become much less sensitive to changes in O2 than at lower O2 concentrations (Belcher et al., 2010). This observation suggests that at the high O2 concentrations modeled for the Cretaceous Period, any small O2-driven fire changes are unlikely to be observable in the fossil fire record, as burn probabilities remain at 99% between 23% and 27% O2. Instead, variations in fire and the associated abundance of charcoal in the fossil record are more likely to be driven by other factors, such as climatic variations, superimposed upon the general high fire activity inferred for much of the Cretaceous (Glasspool and Scott, 2010; Belcher and Hudspith, 2016).

The time period running up to OAE2 is thus a good test case for assessing the probability of shorter time-scale, CO2-driven climate responses to Earth system perturbations. Here, we analyzed the abundance of fossil charcoal across these previously identified intervals of CO2 fluctuation, and we show that CO2-climate driven variations are indicated in the fossil record of fire.

Study Location

Deposits representing the runup to and the beginning of OAE2 can be found well exposed within the middle and upper members of the newly renamed Naturita Formation, previously referred to as the Dakota Formation (Carpenter, 2014), in southern Utah, located ∼500 km west of Pueblo, Colorado, the global stratotype section and point (GSSP) for the Cenomanian-Turonian boundary (Fig. 1; Kennedy et al., 2005; Barclay et al., 2010). Here, high rates of tectonic subsidence and/or high rates of sediment input provide a stratigraphically expanded record of the time period (Laurin and Sageman, 2007), along with high-resolution biostratigraphy and tephrochronology (e.g., Elder, 1985, 1991; Elder and Kirkland, 1985; Elder et al., 1994). During the Cretaceous, this site was located on the western margin of the Western Interior Seaway (Laurin and Sageman, 2007).

The sediments deposited within the studied section display a range of depositional environments indicative of sea-level fluctuations influencing a coastal setting. Facies include back barriers composed of sandstone, claystone, and coal beds, as well as streambeds and floodplains, lagoons and peat swamps, shorefaces, tide-influenced estuaries, and semi-enclosed bays (Laurin and Sageman, 2007; Barclay et al., 2010).

The middle member of the Naturita Formation is largely composed of nonmarine sandstone, mudstone, organic-rich claystone, coal horizons, and carbonaceous mudstone, with occasional marine strata representative of “short-lived marine incursions” (Barclay et al., 2015, p. 215). The upper member of the Naturita Formation is characterized by “marginal marine sandstones, mudstones and shale” with an occasional organic-rich claystone (Barclay et al., 2015, p. 216). The high sedimentation rate of these nearshore depositional sites enabled us to capture potential small-scale variations in wildfire activity, and it allows for a detailed comparison with the associated pulses in CO2 previously recorded by Barclay et al. (2010).

Sample Collection and Processing

The materials used in this study were previously collected and used in Barclay et al. (2010), with the majority of the samples taken from the same set used in Barclay et al. (2010) to produce their CO2 curve. This approach allows for good comparison between the charcoal and stomatal index data.

Samples covering the onset of OAE2 were collected from the Big Hill and Kanarra Mountains localities, which are located at the western margin of the Markagunt Plateau, in a nonmarine to paralic setting of the middle and upper members of the Naturita Formation (formerly Dakota Formation; Figs. 1, 2, and 3). Some of the stomatal index and charcoal abundance data from the pre-OAE2 interval come from localities on the Kaiparowits Plateau (Cottonwood Canyon, Wahweap Wash, and Bald Knoll Mine). The superposition of samples from different localities was inferred from sequence stratigraphic relationships aided by bentonite and biostratigraphic correlation (Laurin and Sageman, 2007; Uličný, 1999). Because sedimentation rates vary along the proximal-distal transect, the depth-domain position of samples relative to a datum is not constant between sections and cannot be used to construct a composite depth domain section without a risk of introducing errors to the temporal sequence of samples. Hence, we avoided combining samples from different localities in the depth domain (i.e., assigning stratigraphic height to a composite section). Instead, the composite logs of δ13C, stomatal index, and charcoal abundance were plotted against time. Construction of the age model is described below.

Samples for this study were chosen from two distinct lithologies: organic-rich claystone and mudstone (with the exception of one siltstone; Fig. 3; Table DR21). Sampling from predominantly just two lithologies was intended to minimize any influence on charcoal abundances from the varying depositional environment, although the changes in inferred depositional environment were also taken into account.

Fossil charcoal was extracted from rock samples using standard palynological acid maceration techniques. The remaining organics were then sieved using a 125 µm mesh, and the >125 µm fraction was collected. All particles within the ∼20 g rock sample were quantified using a binocular microscope, to capture a local- to regional-scale representation of paleowildfire activity (Tinner et al., 2006). Charcoal fragments were identified using key features commonly accepted throughout the literature, which include opacity, cuboid or angular shape, and black/grayish black/dark gray with a graphite-like sheen, in contrast to brownish hues of unpyrolyzed vegetal matter (Patterson et al., 1987; Jones and Chaloner, 1991; Rhodes, 1998; Scott, 2010; Fig. DR1 [see footnote 1]).

Due to the variations in depositional setting (Fig. 3), and to ensure that charcoal abundances were not biased due to overall changes in the flux of organic material to the sites, the total organic matter (TOM) content left after acid maceration was weighed. Charcoal counts were then expressed as per gram of organic matter (g/TOM), as well as per 10 g of bulk rock (Fig. DR2, footnote 1). The similar pattern of variations in charcoal abundances observed in the two graphs (charcoal count vs. charcoal abundance per g/TOM) corroborates the calculated abundance per g/TOM data, indicating that the abundance per g/TOM data provide a good reflection of variations in charcoal during this period.

Age Model

Paralic strata of Big Hill and Kanarra Mountains have been correlated to the astronomical time scale of the Bridge Creek Limestone (Sageman et al., 2006; Ma et al., 2014) using sequence stratigraphic, biostratigraphic, and bentonite and limestone marker-bed correlation (Laurin and Sageman, 2007). In total, eight control points were applied here to the upper member of the Naturita Formation at the reference section of Big Hill (Table DR1, footnote 1). The underlying fluvial strata were further constrained by a bentonite U-Pb age (BH2015; Laurin et al., 2019; Table DR1, footnote 1). The age of the individual samples was estimated using a linear interpolation between age-control points listed in Table DR1 (see footnote 1).

We note that numerical ages presented here can be subject to updates in the future following publication of new radioisotopic and/or astrochronological results. However, none of the potential updates should affect the superposition of samples from Big Hill and Kanarra Mountains and the structure of variations in charcoal abundance discussed here.

A detailed correlation of pre-OAE strata between the Big Hill reference section and localities in the Kaiparowits Plateau (Cottonwood Canyon, Wahweap Wash, and Bald Knoll Mine) is hampered by the lack of suitable correlation markers in the nonmarine succession. Age data from the Kaiparowits Plateau were therefore estimated using a different set of age-control points, including U-Pb ages of two bentonites from the Henrieville locality (Barclay et al., 2015) and a prominent transgressive surface at the top of unit 6A (Uličný, 1999), which is interpreted as coeval with the onset of paralic sedimentation (upper member of the Naturita Formation) at Big Hill (Laurin and Sageman, 2007; Laurin et al., 2019). Correlation between Henrieville and other localities in the Kaiparowits Plateau was based on Uličný (1999). Due to large lateral differences in sedimentation rates reflecting local erosion and differential compaction in the fluvial and estuarine strata of the Kaiparowits Plateau, an additional uncertainty (∼±100 k.y.) should be added to both floating and anchored ages of Henrieville, Cottonwood Canyon, Wahweap Wash, and Bald Knoll Mine samples. These uncertainties, however, do not affect the timing and structure of charcoal maxima discussed in this paper.

We found that the abundance of charcoal per g/TOM at the study site is generally low up to an estimated age of 94.436 Ma (Fig. 4), averaging 189 particles per g/TOM. Charcoal abundances then increase to 5965 particles per g/TOM at 94.436 Ma. Above this rise, charcoal abundances then decline although remaining above “background” abundances at 94.420 Ma to 1446 particles, and increasing again at 94.417 Ma to 6947 particles per g/TOM. At 94.388 Ma, charcoal abundances decline to 186 particles per g/TOM and remain below 35 particles throughout the rest of the section.

It can be seen that the period of increased fire activity appears during a phase characterized by an overall trend of declining pCO2 and following the culmination of a negative excursion in δ13Corganic, at the beginning of a positive δ13Corganic shift. This period is marked as phase 2 (P2) on Figure 4, and it represents the onset of OAE2, which is commonly recognized by a positive shift in carbonate and organic δ13C values (e.g., Jarvis et al., 2011, and references therein; Du Vivier et al., 2014, 2015; Jenkyns et al., 2017). During this phase, the peaks in charcoal abundance appear to occur after the two large pulses of CO2, when CO2 concentrations are estimated to have started to drop. The first pulse in CO2 peaked at ca. 94.474 Ma and corresponds to a low charcoal abundance of 24 particles from the same bulk sample; CO2 concentrations are then estimated to have fallen at 94.436 Ma, correlating to a rise in charcoal abundance to 5965 particles at the same time (also from the same sample). The second CO2 pulse peaking at 94.420 Ma corresponds to a lower charcoal count of 1446 particles, before a decline in estimated CO2 occurred at 94.382 Ma. During this interval of declining CO2, charcoal abundances at 94.417 Ma exhibit a high abundance of 6947 particles, yet they decline to 186 particles at 94.385 Ma.

CO2-Driven Climatic Forcing on Fire

Within the Naturita Formation, Barclay et al. (2010) identified two phases in estimated pCO2 concentrations; phase 1 (P1) corresponds to an overall rising trend in pCO2 between ca. 95.194 Ma and 94.474 Ma that is suggested to have caused enhanced greenhouse climatic conditions. This phase is characterized by generally low charcoal abundances (Fig. 4), except for a brief enhancement at 94.182 Ma, where abundances increase to 916 particles from 143. This appears to correspond closely with a slight decline in stomatal index CO2 estimates at 94.484 Ma. Phase 2 (P2) sees an overall declining trend in pCO2 from 94.474 Ma (Fig. 4; Barclay et al., 2010) that is interrupted by two pulses of CO2 increase. Charcoal abundances peak during P2, occurring shortly after peaks in pCO2 at points where CO2 declines and is sequestered from the atmosphere. Hence, there appears to have been large variations in fire activity at the point when the Earth system was at its most highly disturbed and in a state of flux, where major increases in fire activity seem to correspond to declines in CO2, with the exception of one sample at 94.385 Ma. This suggests a potentially complex relationship between CO2 fluctuations and variations in fire activity.

Are Taphonomic Factors Leading to a Loss or Gain in Charcoal Abundance?

The changing depositional environment inferred for the studied section (Fig. 3; Uličný, 1999; Laurin and Sageman, 2007) could be considered to have influenced the charcoal signal preserved in these sediments. Throughout the studied section, the depositional environment varies between fluvial, tidally influenced estuarine, fluvial-floodplain, semi-enclosed bay, and lagoon/peat swamps (Uličný, 1999; Laurin and Sageman, 2007).

Peat swamps, floodplain, and lagoon settings are ideal catchment areas for charcoal produced by local fires (e.g., Ribeiro et al., 2016; Mackenzie and Moss, 2017). Fluvial environments may act to transport charcoal and other terrigenous material away to nearshore, shallow-marine settings (Nichols et al., 2000; Forbes et al., 2006), such that shallow-marine settings and major river mouths receive the highest concentration of sediment transported by riverine and runoff processes (Nichols et al., 2000). Hence, a higher charcoal abundance may be expected within samples collected from the inferred semi-enclosed bay environment, compared with those from interpreted fluvial deposits. This could therefore explain the higher charcoal abundances observed between 94.417 Ma and 94.436 Ma collected from the Big Hill “Lower” section, and the lower abundances found deposited in the inferred fluvial environment between 95.194 Ma and 94.955 Ma collected from Wahweap Creek, and at 94.520 Ma and 94.482 Ma collected from the Kanarra Mountains and Cottonwood Canyon sections. However, this cannot explain the low charcoal abundance found within the inferred lagoon/peat swamp environments, nor the sample collected at 94.466 Ma in the Kanarra Mountains section, interpreted by Laurin and Sageman (2007) to have been deposited within an intertidal to shallow subtidal restricted estuarine and/or back-barrier tidal flat setting. Furthermore, charcoal abundances between 94.436 Ma and 94.417 Ma vary, yet the depositional environment remains constant in the same Big Hill “Lower” section. Hence, it is important to note that although the increase in charcoal abundance is based on only two sample horizons, variations in charcoal abundances neither correlate with any variations in total organic matter during this period (Fig. DR3A, footnote 1), nor to any changes in the sample collection site, nor to variations in sample lithology (Fig. DR3B, footnote 1). This suggests that the signal is a genuine signal of variations in fire activity that is unrelated to changes in the depositional environment alone, and thus there are likely other factors influencing the wildfire record during this time.

Environmental Shifts across OAE2 and their Potential to Drive Variations in Fire Activity

Phase 1: Rising CO2 Concentrations and Pre-OAE Environmental Changes

OAE2 is hypothesized to have been triggered by a large magmatic episode influencing the Earth system up to ∼200 k.y. prior to the onset of OAE2 (Du Vivier et al., 2014; Ostrander et al., 2017). The initiation of the magmatic episode is hypothesized to have driven a rise in global temperatures, driven by volcanic CO2 outgassing (Frijia and Parente, 2008). As global temperatures increased, the hydrological cycle is hypothesized to have intensified, enabling an increased delivery of nutrients to the oceans, which in turn fueled primary productivity and Corg burial and primed the oceans for deoxygenation (Ostrander et al., 2017).

Records of several isotopic systems (Sr—Frijia and Parente, 2008; Os—Turgeon and Creaser, 2008; Du Vivier et al., 2014; Ca—Du Vivier et al., 2015; Li—Pogge von Strandmann et al., 2013; Tl—Ostrander et al., 2017) have all been used to support a major magmatic episode and release of CO2 into the atmosphere and/or anoxia a few tens of thousands of years prior to the start of the positive carbon isotope excursion that marks the OAE. Small, local perturbations of sedimentary Hg have also been identified during the onset of OAE2, further supporting the hypothesis of a magmatic episode (Scaife et al., 2017; Percival et al., 2018).

The influence of this magmatic episode appears to link with the timing of rising CO2 concentrations in the studied section. During the runup to the OAE, stomatal index proxy estimates from Barclay et al. (2010) suggest a long-term rise in CO2 beginning up to ∼500 k.y. prior to the OAE and peaking at an estimated 500 +400/–180 ppm ∼26 k.y. prior to the start of the OAE (Fig. 4). Throughout this period of generally rising pCO2 concentrations, charcoal abundances at the study site are consistently low, with the exception of a slight increase from 129 particles at 94.487 Ma to 916 particles at ca. 94.482 Ma, which then return to 24 particles at 94.474 Ma. This exception appears to occur close to a slight decline in estimated CO2 occurring at 94.484 Ma.

Increasing CO2 concentrations have been associated with rising global temperatures and also acceleration of the hydrological cycle (Bazzaz, 1990; McElwain et al., 2005; Steinthorsdottir et al., 2012), suggesting that an increasingly moist environment during this time may have suppressed fire activity, producing low charcoal counts.

In the runup to OAE2, temperature proxy data for sea-surface temperatures (SSTs) from the TEX86 method and δ18O data support the idea that global temperatures rose in the prelude to OAE2 (Forster et al., 2007; Sinninghe Damsté et al., 2010; Jarvis et al., 2011). Here, both temperature proxies, taken from European sites and TEX86 from Ocean Drilling Program (ODP) sites at (equatorial) Demerara Rise and around the Newfoundland Basin (proto–North Atlantic, ODP 2176, ∼30°N paleolatitude), indicate a warming trend in the lead up to OAE2. The SST record at Bass River, New Jersey (van Helmond et al., 2014), which is situated ∼2800 km east of our study site, also shows a rising temperature trend, with palynological (terrestrial/marine [T/M] ratios and Paleohystrichophora infusorioides cyst abundances per dry gram of sediment) and inorganic geochemical (Ti/Al) records suggesting enhanced runoff closer to OAE2 (van Helmond et al., 2014). At Bass River, humid and potentially fire-suppressing “greenhouse conditions” are thought to have prevailed in the prelude to the OAE (van Helmond et al., 2014), supporting the hypothesis that increased moisture may have suppressed wildfire activity at this time.

The Cretaceous Period is considered by many to have been a period of high pO2 concentrations (e.g., Berner et al., 2003; Bergman et al., 2004; Mills et al., 2016). At high atmospheric pO2 concentrations (>23%), changes in fuel moisture become the more important driver of variations in fire activity (Watson, 1978; Watson and Lovelock, 2013) because changes in pO2 alone above 23 vol% lead to little variation in burn probability (see fig. 4 in Belcher et al., 2010). However, as pO2 rises, increasingly wet fuel is able to carry a fire because the moisture of extinction (the moisture level at which a fire can no longer spread) is raised.

Belcher and Hudspith (2016) modeled the differences in surface fire behavior for different fuel types in both current ambient and in “super ambient” (26%) O2. They showed that increasing pO2 enables fires to continue to spread above 39% (modern-day) moisture of extinction, and the spread rate is also increased (see fig. 3a in Belcher and Hudspith). However, there are limits to fire spread even at high pO2, with Belcher and Hudspith (2016, their table S2) estimating the moisture of extinction to be ∼80% at 26% pO2.

Belcher and Hudspith (2016) used the equation Mex = 8O2 – 128 (Watson and Lovelock, 2013) to estimate the moisture of extinction for different atmospheric O2 concentrations in their models. By assuming that pO2 was between 25% and 27% during OAE2 (Mills et al., 2016), the moisture of extinction would likely have been between 77% and 88%. If the Belcher and Hudspith (2016) surface fire spread data for different fuel types are extrapolated across this range, a rise of just 10% fuel moisture would likely cease fire spread according to their low-fuel-load models (i.e., conifer litter [CL2] and weedy understory [WUn2]) and would halve the spread rate in higher-fuel-load models (e.g., dense shrub and fern understories [SUn2 and FUn2] models). Furthermore, the intensity of surface fires (see fig. 3C in Belcher and Hudspith, 2016), would also fall significantly by approximately two thirds, from ∼900 kWm−1 to 300 kWm−1, as fuel moisture increased.

Most surface fires are started in dead and drier surface fuels. These surface fires then have to be able to dry the overstory canopy fuels below the required moisture of extinction before significant crown fires can be ignited and carried. This is significant because fire-line intensity relates to flame heights and the ability of heat energy from surface fires to scorch the forest canopy and transition to crown fires (see Belcher and Hudspith, 2016). Lowering of fire-line intensity leads to lower flame heights, thus limiting the ability of surface fires to dry, scorch, and ignite the canopy. Therefore, if the enhancement of the hydrological cycle led to a rise in fuel moisture of just 10%, in the phases leading up to onset of the OAE, it would be anticipated to have a strong negative influence on large-fire potential, despite the estimated high pO2, and may therefore explain the relatively low abundances of charcoal in phase 1.

Interestingly, prior to the peak in pCO2 concentrations at ca. 94.474 Ma, a brief decline in pCO2 concentrations occurred at ca. 94.484 Ma (Fig. 4), estimated at a mean decline of ∼52 ppm (Hypodaphnis zenkeri estimated decline from the stomatal index proxy of 59 ppm; Laurus nobilis estimated decline of 45 ppm). This decline in CO2 occurs at a similar time to where we see a slight increase in charcoal abundances to 916 particles at ca. 94.482 Ma, before swiftly returning to low abundances at 94.474 Ma, corresponding with the high inferred CO2 concentrations occurring within the same bulk sample at 94.474 Ma. Although based on only one data point, this is particularly interesting because it may indicate that geologically swift fluctuations in inferred CO2 in the runup to OAE2 may have been capable of leading to observable short-term variations in fire activity, possibly driven by regional fluctuations in climate.

During this time, SST calculations using the TEX86 proxy in the Newfoundland Basin, offshore Canada (Forster et al., 2007), and the Bass River borehole, New Jersey (van Helmond et al., 2014), suggest SST continued to increase throughout the runup to the OAE (van Helmond et al., 2014). However, within the Bass River borehole, New Jersey, variations in SST can be observed, showing brief periods of reversal within the overall warming trend, e.g., at 595.1 m below seafloor (mbsf), 594.6 mbsf, 593.6 mbsf, and even at the inferred onset of the OAE at 593.3 mbsf (see van Helmond et al., 2014, their fig. 2). Furthermore, van Helmond et al. (2014) indicated that temperature and hydrology changed consistently across OAE2, suggesting that this small variation in inferred fire activity may not be an anomaly; however, further study would be required to corroborate whether these small shifts in SST correlate to the small shifts in CO2 and charcoal found within our studied section.

Phase 2: Declining CO2 Concentrations and OAE Environmental Changes

The second phase (P2) begins with the period of peak pCO2 at 94.474 Ma and ends with a CO2 minimum at 94.382 Ma (Fig. 4). Throughout phase P2, pCO2 generally declines, except for a brief enhancement of pCO2 at ca. 94.420 Ma. However, phase P2 suggests evidence of a large increase in fire activity. During P2, the abundance of charcoal per g/TOM suddenly increases to 5965 particles per g/TOM per year at 94.436 Ma. Abundances appear to decline briefly at 94.420 Ma to 1446 particles g/TOM coincident with the brief rise in pCO2, and increase and peak at 6947 particles g/TOM at 94.417 Ma, before falling to just 186 particles g/TOM at 94.385 Ma (Fig. 4). The overall decline in pCO2 concentrations observed at the study site during P2 has been attributed to increased C-sequestration by primary production that was capable of drawing down twice the magnitude of CO2 compared to the first CO2 drawdown within P1 (Barclay et al., 2010). It is during this phase that the shift toward positive δ13C isotopes, correlative between terrestrial and marine records, marking the start of the OAE occurs.

By the initiation of OAE2, anoxic conditions are hypothesized to have expanded by ∼40% since the start of deoxygenation (beginning ∼43 k.y. prior to the event; Ostrander et al., 2017). During the onset, enhanced production and preservation of Corg occurring in euxinic “hotspots” are hypothesized to have driven the positive shift in C-isotope values that marks the onset of OAE2 (Owens et al., 2016).

Interestingly, during the initiation of OAE2, geochemical sediment analysis indicates a noticeable increase in P at the GSSP at Pueblo, where all P species measured (PFe, Pauthigenic, Pdetrital, Porg) reach peak values at the onset of the event (from 0.10 to 0.20 [PFe], 0.50 to 5.00 [Pauthigenic], 0.25 to 0.50 [Pdetrital], and from 0.08 to 0.20 [Porg]; Mort et al., 2007; Flaum, 2008), likely driven by increased delivery and regeneration from anoxic sediments (Monteiro et al., 2012). This is intriguing because increases in P fluxes to the oceans have been associated with wildfires, and fire seasons (Kump, 1988; Mahowald et al., 2008).

A rise in sea level is inferred at the same time that could have aided the delivery of nutrients from the continents to the oceans (Bjerrum et al., 2006). Changes in sea level are also speculated to have caused marine hiatuses in key sections, e.g., those observed in the Eastbourne section (at the base of the Plenus Marl Member; Jarvis et al., 2011) and the Portland core (within the uppermost Hartland Shale). The Portland core has an estimated hiatus of up to ∼40 k.y. (Du Vivier et al., 2014, 2015). When the section in this study is correlated to the Portland core (Fig. 2), it is clear that the period of high charcoal found in the Naturita Formation likely corresponds to this hiatus, spanning ∼40 k.y. (Figs. 2 and 4). Therefore, one hypothesis for the rise in charcoal abundances during this time could be an increased flux of charcoal to the study sites associated with a rise in sea level (Laurin and Sageman, 2007; Barclay et al., 2010). Here, a rise in sea level may have enabled removal of charcoal from the newly flooded land surface, providing a “spike” in charcoal abundances during the initiation of OAE2. However, this seems unlikely because charcoal abundances prior to the onset are consistently low, making leaching from nearby land surfaces an unlikely source of charcoal. The lack of charcoal abundances found throughout the pre-OAE period similarly makes reworking at the OAE onset an unlikely source of charcoal.

Over geological time scales, atmospheric CO2 concentrations are arguably a primary driver of SST fluctuations, with increasing CO2 linked with rising SSTs and drawdown periods associated with cooling (Jarvis et al., 2011, and references therein).

Interestingly, across the onset of OAE2, SST records from the European sites at Pont d’Issole, SE France, and Gröbern, N Germany, which exhibit no apparent hiatus, indicate a brief period of cooling in SSTs (Jarvis et al., 2011). Similarly, the Bass River borehole, New Jersey, shows a small decline in SST across the OAE2 onset, before immediately rising by ∼0.5 °C during the OAE (van Helmond et al., 2014). The section at Pont d’Issole suggests a more significant cooling period across the OAE onset, with an estimated fall in SST of ∼4 °C (Jarvis et al., 2011, fig. 8), while the Bass River borehole section shows smaller declines of just ∼0.24 °C (van Helmond et al., 2014). These variations indicate that temperature and likely climate may have been regional in distribution, and they suggest that rising global temperatures observed in the runup to OAE2 slowed, and possibly even briefly reversed across the initiation of the OAE.

The Bass River borehole is located nearly ∼3000 km east of our study site, and it is the closest site preserving a SST record across the onset of OAE2. The variations in fire activity, as indicated by changes in the abundance of charcoal, appear to relate directly to changes in pCO2 across the onset of OAE2, with the exception of the lower abundance found at 94.385 Ma. The first “peak” in charcoal of 5965 particles occurs at 94.436 Ma and corresponds with a fall in stomatal index–estimated CO2 from 461 ppm to 389 ppm and from 515 ppm to 464 ppm in Hypodaphnis zenkeri and Laurus nobilis, respectively. A decline in charcoal abundance then occurs at 94.420 Ma, along with a rise in inferred CO2 of 56 ppm (mean from Hypodaphnis zenkeri and Laurus nobilis), before charcoal abundances rise again to 6947 particles at 94.417 Ma as CO2 declines again. It should be indicated that these fluctuations in stomatal index occur in the same samples as we have extracted our charcoal. Therefore, it appears that, at least regionally, changes in CO2 appear to have directly influenced wildfire activity recorded at the same study section. This suggests a CO2-climate-driven influence on wildfire across the onset of OAE2, and that the drawdown in CO2 across the onset of OAE2 may have not only initiated temporary cooling, but also suppressed the hydrological cycle, enabling wildfire activity to briefly increase.

This seems to be a realistic interpretation because drier climate conditions can lead to a rise in wildfire activity by reducing the moisture content of fuels (Viegas et al., 1992; Collins et al., 2007). A decline in the moisture content of vegetation, below the required moisture of extinction, would enable an increase in the rate of fire spread and fire intensity and also increase the likelihood that a surface fire would transition to a crown fire (Belcher and Hudspith, 2016). It could therefore be suggested that, during the phases of CO2 drawdown, climate conditions became more favorable to wildfire activity. The decline in charcoal abundance at 94.420 Ma, coeval with evidence for a small pulse in CO2, may have been driven by a small, swift reversal in the cooling and drying trend and may have provided the ideal conditions for regrowth of significant fuel following the period of higher fire activity.

Geologically swift reversals in climate trends can be observed throughout the OAE2 period, with a short-term cooling identified at the base of the Metoicoceras geslinianum zone in European section Pont d’Issole (see fig. 8 in Jarvis et al., 2011) estimated to have lasted ∼40 k.y. (Jarvis et al., 2011). Similarly, the Plenus cold event, occurring toward the top of our sampled section (Fig. 2), is estimated to have lasted <60 k.y. before subsequent warming (Jarvis et al., 2011), and it has recently been suggested to have been characterized by two cooler intervals, punctuated by an intervening warmer period (Jenkyns et al., 2017). Thus, there is good evidence that shifts in climate are possible within the short time scales of tens of thousands of years observed here in the charcoal record.

However, the return to lower charcoal abundances of just 186 particles at 94.385 Ma appears anomalous when compared with the previous higher charcoal abundances found during this inferred CO2 drawdown phase, and it may therefore suggest a more complex relationship.

Several studies have suggested an important link between Milankovitch cyclicity and astronomically paced insolation and the timing of OAE2 initiation. Milankovitch cycles exert a control on the seasonal and latitudinal distribution of solar radiation, and they can thus influence the amplitude of seasonal variations over varying periodicities of precession (∼20 k.y.), obliquity (∼40 k.y.), and eccentricity (∼95–124 k.y. and ∼405 k.y.; Berger et al., 1993; Laskar et al., 2004). Changes in the precession cycle exert a control over the strength of hemisphere-specific seasonality, in particular, influencing temperature and precipitation patterns (e.g., Kutzbach and Gailimore, 1994; Flögel, 2002). The climatic effect of precession is linked to Earth’s orbital eccentricity, whereby changes between more and less elliptical orbits modulate the amplitude of precession-paced seasonality fluctuations.

Obliquity determines the strength of the temperature gradient between latitudes (Raymo and Nisancioglu, 2003), as well as seasonal insolation in high-latitude regions and integrated annual insolation at all latitudes (see discussion in Huybers, 2006).

Sites situated around the Cretaceous Western Interior Seaway, including the GSSP at Pueblo, show evidence of Milankovitch-scale forcing around the time of OAE2 (e.g., Sageman et al., 1997; Flögel, 2002; Eldrett et al., 2015; Charbonnier et al., 2018). This evidence suggests that despite the inferred warm, equitable climates of the Cretaceous greenhouse world, the regional climate around the Western Interior Seaway appears to have responded to changes in solar insolation driven by Milankovitch cycles (Eldrett et al., 2015).

The initiation of OAE2 is suggested to have coincided with a maximum in the 405 k.y. eccentricity cycle (Laurin et al., 2016). During periods of precession-paced insolation maxima, increased storm frequency and intensity, precipitation, and temperatures are hypothesized in the Western Interior Seaway (Eldrett et al., 2015). By contrast, during periods of precession-paced insolation minima, lower temperatures, precipitation, and runoff are inferred (Eldrett et al., 2015). High eccentricity at the onset of OAE2 (Laurin et al., 2016) therefore implies that our study site may have experienced pronounced extremes (both maxima and minima) in seasonality.

The two charcoal samples at 94.436 Ma and 94.417 Ma exhibiting higher abundances (5965 particles and 6947 particles, respectively) during the OAE initiation occur over an estimated ∼19 k.y., sitting within the precession period band of 17 k.y. to 21.5 k.y. The following low charcoal abundance of 182 particles, also occurring during the period of inferred CO2 drawdown, is found ∼32 k.y. following the high charcoal sample at 94.417 Ma. These time scales suggest that precessional forcing may have been capable of influencing the fire record, possibly through regional climate fluctuations. Here, the low charcoal abundance at 94.385 Ma equates to ∼1.5 precession cycles following the high charcoal abundances and would therefore correspond to an opposing precessional index (e.g., perihelion vs. aphelion seasons), enhanced by the maximum 405 k.y. eccentricity (Laurin et al., 2016). Hence, a switch between extreme seasonality and decreased seasonality may explain a dampening of the inferred wildfire record at this point. Here, the higher charcoal abundances may also be explained by a reduction in precipitation occurring during periods of precession-paced insolation minima, while the lower charcoal abundance may have occurred during precession-paced insolation maxima, a period of hypothesized increased precipitation and runoff (Eldrett et al., 2015). If correct, this may suggest that Milankovitch climate-related forcing on the wildfire record may be apparent during the onset of OAE2, coinciding with a period of inferred declining CO2.

This hypothesis would appear to support a recent prediction by Charbonnier et al. (2018), who inferred that detrital input (using changes in aluminum content, titanium content, and magnetic susceptibility) during OAE2 was likely driven by intensification of continental weathering and acceleration of the hydrological cycle related to increased volcanic activity and the release of greenhouse gases, which may have been modulated by the obliquity Milankovitch cycle.

Following the CO2 drawdown at 94.436 Ma, the Barclay et al. (2010), pCO2 proxy data suggest that concentrations did not remain at minimum values, but instead, they likely increased again within the OAE. This likely indicates a return to warm, humid, fire-suppressing greenhouse conditions over the remainder of the anoxic interval. Support for this also comes from geochemical proxy estimates for SSTs (e.g., Jarvis et al., 2011) and continental weathering (e.g., δ7Li isotope values; Jenkyns et al., 2017), which indicate that much of the anoxic interval was characterized by warm, humid, greenhouse conditions that likely enhanced the hydrological cycle, increasing weathering and runoff, and facilitating the delivery of limiting nutrients to the marine system (Jarvis et al., 2011; van Helmond et al., 2014; Jenkyns et al., 2017). These conditions are anticipated to have prevailed for the remainder of the OAE (Jarvis et al., 2011; van Helmond et al., 2014), with the exception of the geologically brief Plenus cold event.

One explanation could therefore be that the return to low inferred fire activity was due to the return to warm, humid conditions, where high fuel moisture once again suppressed fire ignition and spread, as evidenced by the consistently low charcoal abundance found in the remainder of the OAE section that was sampled, with the exception of one sample at 94.230 Ma, which correlates with the Plenus cold event (Fig. 2).

Another explanation could be that the samples collected happen to coincide with periods of increased precipitation and runoff driven by Milankovitch cyclicity. However, the decline in charcoal abundance found at 94.420 Ma, coinciding exactly with the estimated rise in CO2 from stomatal indexes between the two high charcoal counts, does not occur over the precessional astronomical time scale and instead occurs just 3 k.y. prior to the high charcoal abundance found at 94.417 Ma. Hence, we propose that “volcanically induced” CO2 fluctuations likely were the primary driver for variations in inferred fire activity observed across the OAE2 onset. Across the OAE onset and during the period of inferred declining CO2 concentrations, Milankovitch forcing associated with the OAE onset may have been capable of modifying the wildfire signal at this crucial time.

Implications of Increased Fire Activity across the Onset of OAE2

Fire is an important feedback mechanism in the natural world, and it has been shown to play a crucial role in the long-term regulation of atmospheric O2 concentrations (e.g., biogeochemical model COPSE; Bergman et al., 2004). This regulation is driven by fluctuations in the amount of Corg burial, which in turn is driven by the flux and availability of limiting nutrients such as P. The key role that fire plays as a feedback process via a control over P weathering/distribution to the marine system has been illustrated by Kump (1988) and Lenton (2013), where debate remains as to whether fires aid in the redistribution of P from the land to the oceans (e.g., Kump, 1988), or whether fires act to limit P weathering through the suppression of land vegetation (e.g., Lenton and Watson, 2000; Bergman et al., 2004).

In the modern day, studies have shown that, over short time scales, the concentration of P in nearby streams increases significantly following a fire, suggesting that fire can actually aid in the redistribution of P from the land to nearby streams, which feed into the marine environment. For example, Spencer and Hauer (1991) found that P concentrations in nearby streams increased from 3 μg P/L prefire to 135 μg P/L within 24 h postfire. Similarly, Taylor et al. (1993) found an increase in P from 0.05 mg/L prefire and flood to 0.10 g/L postfire and flood. Fires can also influence the supply of P via aerosol inputs to the atmosphere through smoke, causing atmospheric deposition of P to open-ocean sites. For example, the Amazon Basin appears to be losing P to the atmosphere, where it is being deposited in neighboring regions and adjacent oceans downwind. Of this flux, 23% is attributed to fire (Mahowald et al., 2005).

At the onset of OAE2, bottom waters at the U.S. Geological Survey (USGS) #1 Portland core site are interpreted to have been relatively oxic due to the presence of extensive bioturbation (Savrda, 1998) and a decline in the redox-sensitive trace-element molybdenum (Mo; Meyers et al., 2005). Although seemingly contradictory with the start of an anoxic event, locally the Western Interior Seaway basin is hypothesized to have been well mixed, preventing stratification and anoxia in bottom waters (Flaum, 2008). A switch to increasingly oxic waters could increase P released through the initial oxidation of organic matter, in particular producing an increase in precipitated Pauthigenic (Mort et al., 2007).

Pdetrital values, which Flaum (2008) argued can be utilized as a proxy for riverine and hence terrestrial P input to marine surface waters, are also high across the onset of the OAE at Pueblo (Mort et al., 2007). Although it is likely that the flux in P observed during the initiation of the OAE2 (Mort et al., 2007) may thus have been largely driven by local conditions, including oxidation of organic matter, enhancement of the hydrological cycle, and a rise in sea level (causing a flux in Pdetrital), it is interesting that, coinciding with the flux in P found at nearby study sites (Mort et al., 2007), and the OAE initiation, there is also evidence for an increase in wildfire activity, which could also contribute to the Pdetrital flux.

This phenomenon of increased charcoal abundances at the onset of OAE2 has also been observed at the start of the Jurassic Toarcian OAE (T-OAE; Baker et al., 2017). This is despite the fact that the T-OAE and Cretaceous OAE2 occurred during very different periods of Earth’s history, in particular, with respect to atmospheric O2 concentrations. The T-OAE is hypothesized to have occurred during a time when estimated atmospheric O2 concentrations were ∼20% (Lenton, 2013). By comparison, OAE2 occurred during a time when O2 concentration estimates were much higher (Wilson and Norris, 2001; Mills et al., 2016). Despite this difference, the charcoal record across OAE2 indicates that fire activity responds similarly when compared to that observed by Baker et al. (2017) across the initiation of the T-OAE. For example, fire activity rises across the initiation of the T-OAE, at both the Peniche locality and the Mochras core, as indicated by enhanced charcoal abundances at the culmination of the early positive δ13C trend. During the early stages of OAE 2 in the Naturita Formation, inferred fire activity appears to have responded similarly, with increased abundances occurring during the initiation.

This is particularly interesting because it suggests that CO2-induced climate changes may have been capable of driving variations in paleofire activity during the initiation and early stages of OAEs, irrespective of the different “background” atmospheric O2 concentrations. Furthermore, a rise in fire activity during the onset of OAEs may have aided in the short term, the delivery of limiting nutrients, e.g., including Pdetrital, from the continents to the oceans, tipping the ocean more rapidly toward anoxia.

We therefore suggest that future studies may want to consider the importance of changes in wildfire as being critical toward altering Earth system balances across OAEs, such that variations in fire may provide a link that tips the oceans into anoxia (as described herein), as well as later aiding recovery (Baker et al., 2017).

Could Increased Wildfire Activity Have Aided in the Movement toward Negative C-Isotope Values?

The rise in inferred wildfire activity during the early stages of OAE2 appears to have occurred coeval with a period of generally low C-isotope values (Fig. 4; Barclay et al., 2010). Sustained extensive biomass burning over thousand-year time scales is hypothesized to have been capable of producing pronounced negative excursions in the δ13C record (e.g., Finkelstein et al., 2006). Finkelstein et al. (2006) conducted model runs of biomass burning under Cretaceous conditions of 2× present atmospheric level CO2 and 20 °C ocean temperatures, over varying time scales, hypothesizing that sustained combustion of 0.08% of the global peat reservoir over ∼1000 yr could be capable of producing a negative excursion of up to –2.4‰ within the atmospheric reservoir. Hence, one hypothesis could be that the increase in inferred wildfire activity, and hence biomass burning, may have also aided in the movement toward negative C-isotope values during this time.

However, the movement toward “peak” negative C-isotope values occurred prior to the first peak in charcoal abundances, occurring ca. 94.474 Ma (with a δ13C value of –25.27‰) and reaching –25.99‰ at ca. 94.449 Ma. The first peak in charcoal abundances occurred at ca. 94.436 Ma, remaining high until 94.417 Ma. During this period of high charcoal, δ13C values moved toward more positive values, from –25.99‰ at ca. 94.449 Ma to –25.5‰ at ca. 94.417 Ma.

One explanation could be that the release of isotopically light C from biomass burning may have dampened the return to positive C-isotope values, although not by enough to prevent the δ13C values from rising. Owens et al. (2018) estimated the amount of organic carbon burial from known global distributions of organic-rich facies deposited during OAE2. The model results indicated that the known global deposits of organic matter would only account for ∼40% of the recorded δ13C excursion observed. Hence, it appears that even when expanded to unknown portions of the world’s oceans, significant amounts of organic carbon must have been buried in reservoirs not currently constrained. It is possible that this “unaccounted” reservoir may therefore have been larger than that predicted by Owens et al. (2018) in order to be able to further account for and buffer the influx of negative C-isotope values from biomass burning. However, further study would be required to test the geographical extent of the inferred increase in wildfire, compared with the local-regional scale observed within the Naturita Formation, as such an influence on the C-isotope record would likely require a larger, more global signal of increased wildfire occurrence that should be observable within the fossil fire record.

Throughout the prelude to and early stages of OAE2, fluctuations in volcanically derived CO2 are hypothesized to have led to significant changes in the global climate. The influx of CO2 from large igneous province activity during this time was followed by a drawdown period occurring across the OAE initiation, which was associated with increased C sequestration by primary producers (Barclay et al., 2010) and increased silicate weathering (Jenkyns et al., 2017). Here, we find that the period leading up to the OAE initiation, coeval with evidence of rising CO2, was characterized by consistently low charcoal abundances, while the inferred CO2 drawdown appears to have been coeval with increased abundances of charcoal and therefore enhanced fire activity. A brief pulse in CO2 is hypothesized to have occurred during the overall drawdown phase, which occurred coeval with a brief decline in charcoal abundances, before abundances then rose again as CO2 continued to decline. Published proxy data indicate that during this phase of CO2 drawdown, a brief period of global cooling, and thus drying of the climate, may have occurred (Jarvis et al., 2011; van Helmond et al., 2014). Hence, it appears that CO2- and climate-related variations in inferred fire activity may be observable across the onset of OAE2.

Our study highlights that climate-driven changes in wildfire activity remain possible even under the modeled high pO2 concentrations of the Cretaceous Period, and these results implicate changes in wildfire activity as being drivers of both OAE initiation and recovery (cf. Baker et al., 2017). Our analysis provides strong support for the occurrence of geologically rapid climate shifts in the lead up to and during the early stages of OAE2, where changes on land (e.g., fire activity) likely contributed significantly to the onset of OAE2.

This work was supported by the Natural Environment Research Council, through funding a studentship grant (grant number NE/L 501669/1) to Baker. Belcher acknowledges funding through a European Research Council (ERC) Started Grant (grant ERC-2013-StG-335891-ECOFLAM). Laurin acknowledges funding from the Czech Science Foundation, project 17-10982S. Barclay acknowledges a National Science Foundation (NSF) grant EAR-0643290 to B. Sageman and J. McElwain that supported his dissertation work, during which the samples that were analyzed in this study were collected, and the Field Museum, where the samples are curated. Barclay also acknowledges support from The Peter Buck Foundation and the department at the Smithsonian’s National Museum of Natural History. We thank Mark Grosvenor for his assistance in the laboratory. We also thank one anonymous reviewer and J. Owens for their useful comments on the manuscript.

1GSA Data Repository item 2019243, comprising raw data for charcoal counts and age control points, charcoal image, and cross plots of charcoal abundance vs. organic matter contents, is available at or by request to
Science Editor: Bradley S. Singer
Associate Editor: Ganqing Jiang
Gold Open Access: This paper is published under the terms of the CC-BY license.