Centennial to millennial variability of greenhouse climate across the mid-Cenomanian event

Climate change models with elevated atmospheric CO₂ indicate that natural forcing contributions from total solar irradiance (TSI) variations are likely to be small (Myhre et al., 2013). However, modeling suggests that TSI variations have the potential to influence ocean temperature (Seidenglanz et al., 2012). TSI variations occur across centennial to millennial (C-M) time scales, although knowledge is limited by the availability of data. Holocene solar activity proxies reveal variations that include the Schwabe (11 yr), Gleissberg (80–100 yr), de Vries-Suess (200–250 yr), Eddy (∼1000 yr), and Halstatt (or Bray) (∼2400 yr) cycles (Usoskin, 2017). Importantly, these cycles intersect frequency bands of unforced internal climate oscillations generated by global climate models (Dijkstra and von der Heyt, 2017; von der Heydt et al., 2021). A fundamental question is whether C-M paleoclimate variability reflects responses to unforced internal climate oscillations and/or forcing from TSI variations.

There is a rapidly growing database of Holocene C-M paleoclimate records, from tree rings, speleothems, ice, and varved lake and marine sediment sequences (http://pastglobalchanges.org/data/databases). Many of these records exhibit cycles with periodicities reminiscent of those of the solar activity proxies indicated above. In contrast, pre-Holocene C-M paleoclimate records from sediment are exceedingly rare, due to low sediment accumulation rates, postdepositional mixing, incomplete preservation, and the lack of precise chronologies.

We investigated C-M paleoclimate variability immediately preceding and during the lower part of the mid-Cenomanian event (MCE; ca. 96.6–96.2 Ma; Eldrett et al., 2015a), a “greenhouse” period characterized by high atmospheric CO₂ concentrations (500–1500 ppmv) and warm global sea-surface temperatures (>35 °C; Wang et al., 2014; O’Brien et al., 2017). Our goal was to characterize C-M variability in the climate system from the marine perspective, compare its spectral features with the Holocene TSI record, and gauge the relative contributions of external TSI forcing versus internal climate variability.

The C-M paleoclimate record examined in this study is from the Shell Iona-1 research core (west Texas), located in the southern gateway of the Cretaceous Western Interior Seaway (KWIS, North America; Fig. S1 in the Supplemental Material¹). The studied interval is characterized by suboxic-anoxic depositional conditions under a restricted and stratified water column (Eldrett et al., 2017), which experienced temporary photic zone euxinia (Sun et al., 2016). Meter-scale alternations in redox state are also recorded between limestone (L) and marlstone (M) beds, with L-beds reflecting periods of greater water-mass ventilation and current activity in response to precession forcing (Eldrett et al., 2015a, 2015b). We analyzed nine M-beds (M1–M3, M6–M11) spanning just prior to and across the lower part of the MCE interval. (Fig. S2). Within each M-bed, millimeter-scale light-dark laminae have been preserved. Sedimentation rates for these nine M-beds were obtained by Ma et al. (2020); the time interval represented by M11, selected for more detailed analyses, was estimated as 18.933 k.y., with a sedimentation rate of 1.02 cm/k.y. (Ma et al., 2020).

We conducted detailed petrographic and micropaleontological analyses on M11 (Fig. 1; Figs. S4 and S5). Major and minor elements were measured using an X-ray fluorescence (XRF) device on an Itrax core scanner at 200 μm resolution (Δt = ∼20 yr; Fig. S6) at the National Oceanography Centre, University of Southampton, Southampton, UK. Grayscale data were measured at 230 μm resolution (Δt = ∼23 yr) along the image shown in Figure S2.

Principal component analysis (PCA; prcomp function in R, https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/prcomp) was conducted on six major elements (Ca, Si, Fe, K, Ti, and S). Multitaper power spectral analysis with classical red noise modeling and harmonic F-ratio testing was applied with three 2π prolate tapers to evaluate significant cycles in the time series shown in Figure 1. (Scripts are available at: http://mason.gmu.edu/∼lhinnov/cyclostratigraphytools.html.)

We then compared the M11 time series with the Holocene TSI time series for the past 9400 yr (Steinhilber et al., 2013) to determine similarities and/or differences among the records.

The M11 succession represents mainly hemipelagic deposition. M11 preserves millimeter-scale alternations between white to light-gray carbonate-rich laminae (termed “light”) and dark-gray to black organic carbon–rich laminae with 4–16 wt% total organic carbon (termed “dark”). The light laminae are composed primarily of foraminifera, whereas dark laminae are composed of abundant fecal pellets, some foraminifera, and a lack of bioturbation features. Very few light-dark laminae sets exhibit internal truncations, downlapping structures, and/or low-angle ripples that would indicate advective sediment transport, and so they are interpreted as single bed-load depositional events (Minisini et al., 2018). Erosional surfaces at the bases of light laminae were not observed to cut through multiple laminae, indicating limited impact of hiatuses on the time-series analysis.

The light laminae contain diverse assemblages of foraminifera including larger adult-sized and smaller juvenile-sized Hedbergellids and Heterohelicids from the upper mixed layer, Whiteinellids (lower mixed layer), and Rotaliporids (deep dwellers from within or below the thermocline; Gebhardt et al., 2010) combined with minor occurrences of calcispheres, inoceramids, crinoids, and collophane (Fig. S5). The dark laminae contain rare benthic foraminifera, and low-diversity planktonic foraminifera assemblages dominated by juveniles of the r-strategist (capable of large population increases in a short period of time) Heterohelicid. The mixed grain-size population of planktonic foraminifera in the light to dark laminae sets is not consistent with hydrodynamic sorting associated with event deposition and bed-load advective transport.

The light laminae are relatively enriched in Ca, are depleted in the other five major elements (Fig. S6), and have low concentrations of redox-sensitive trace metals and slightly increased Mn compared to the dark laminae (Fig. S7). The dark laminae are relatively depleted in Ca, are enriched in the other five major elements (Fig. S6), and contain relatively high abundances of redox-sensitive trace metals V and Ni (Fig. S7).

The first two principal components account for 68.6% (PC1) and 17.2% (PC2); i.e., a total of 85.8% of the total variance (Fig. S8b). For PC1, the six elements have comparable loadings in terms of absolute values, but with a negative Ca loading (high Ca with low PC1, where PC1 is interpreted as a signal of carbonate productivity; Figs. S8a and S8c). In the PC2 axis (Figs. S8a and S8c), Ca and Ti are almost zero; Fe and S have negative loading values; and Si and K have positive loading values. These results suggest that PC2 is associated with redox elements (Fe and S) versus terrigenous elements (Si and K).

The M11 proxy power spectra are arranged in order of increasing high-frequency power contributions (Fig. 2; Figs. S9–S11). PC2 has the simplest spectrum, dominated by a 10,784 yr (110 mm) cycle and a low-power 412 yr (4.2 mm) cycle. PC1 expresses a 398 yr (4.06 mm) cycle, with lower power from 316 yr (3.22 mm) to 198 yr (2.02 mm), and still lower power from 167 yr (1.70 mm) to 109 yr (1.11 mm). High-power peaks at 1647 yr (16.8 mm) and 4843 yr (49.4 mm) do not exceed the 99% significance level for red noise. The grayscale spectrum has sharp peaks at 393 yr (4.01 mm), 211 yr (2.15 mm), 168 yr (1.71 mm), and 127 yr (1.30 mm), a high-power peak at 1686 yr (17.2 mm), and another peak at lower power at 906 yr (9.24 mm), neither exceeding the 99% significance level. The grayscale spectrum resolves variability out to 78 yr (0.80 mm) cycles and is not as “red” as the PC1 spectrum. Finally, the TSI spectrum has a large peak at 203 yr and lower-power peaks in the band from 151 yr to 87 yr. Peaks at periods longer than 203 yr do not exceed the 99% significance level. Spectral analyses on eight M-beds above M11 also showed ∼2 mm cycles that represent an ∼200 yr cycle (Fig. S3; Tables S1 and S2). Furthermore, a simple analysis of number of laminae pairs (average thickness 1.96–2.34 mm) in each M-bed suggests that the duration of the individual laminae pairs is ∼200 yr (Table S1).

The light laminae contain foraminifera assemblages that are interpreted to occupy habitats throughout the water column, suggesting relatively oxygenated, biologically favorable environmental conditions within the water column. Low concentrations of redox-sensitive trace metals and slightly increased manganese also point to oxygenated environmental conditions. The dark laminae have low-diversity planktonic foraminifera assemblages dominated by juveniles of the r-strategist Heterohelicid, rare benthic foraminifera, and a lack bioturbation of features. These features and the high abundance of redox-sensitive trace metals are consistent with an environmentally stressed eutrophic environment in the water column, with suboxic-anoxic conditions at the seafloor (Brumsack, 2006; Eldrett et al., 2014, 2015b, 2017; Minisini et al., 2018).

PC1 is interpreted as a carbonate productivity signal, correlating well to the grayscale time series, with high Ca (low PC1) correlating to high grayscale values. PC1 displays dramatic swings between high and low values, suggesting rapid changes from anoxic to oxygenated conditions. Individual layers of foraminifera are responsible for the high magnitudes of these changes, many of which are interpreted to occur at a quasi-annual time scale (Fig. 1). PC2 is interpreted as a seafloor redox signal, with Fe and S indicating the presence of pyrite, and with higher values usually (but not always) coinciding with dark laminae (Fig. 1).

The 80–100 yr and 200–230 yr cycles in the PC1 and grayscale spectra are reminiscent of similar cycles in the TSI spectrum attributed to Gleissberg and de Vries–Suess cycles, respectively (Fig. 2). The M11 spectra also share a strong peak with a 400 yr cycle that has strong side lobes, suggesting some frequency instability; the TSI spectrum has two peaks at 477 yr and 351 yr, possibly a 400 yr cycle analogue. A strong 400 yr cycle appears in Holocene paleoclimate records in North America (Yu and Ito, 2002) and Asia (Wu et al., 2009), matched to atmospheric radiocarbon activity (a solar activity proxy) time series. The TSI, grayscale, and PC2 spectra (Fig. 2) also share a broad and low-power ∼900–1000 yr cycle.

The M11 centennial-scale cycles suggest variations in water-column stratification, redox conditions, and marine community structure in response to changes in storm activity and intensity. Storms periodically reached sufficient intensity to trigger temporary destratification of the water column before returning to the stratified and anoxic environment that normally characterized this location (Fig. 3). The rapid, short-lived, ∼0.5 W/m² decreases in the TSI time series mimic spikes in grayscale intensity (to white) and PC1 minima (Fig. 1). If Cenomanian and Holocene TSI variations were tantamount, lower TSI would be linked with light laminae and times of increased storminess and water-column ventilation. More intense storms could have been generated by stronger winds produced by an increased meridional thermal gradient and ocean-land thermal contrast when TSI was low (Weng, 2012).

The M11 spectra share multiple spectral features with the Holocene TSI spectrum. This raises an interesting (two-part) question: If the M11 spectra reflect mid-Cenomanian internal climate variability, with greenhouse conditions and a very different ocean-continent configuration from the Holocene, should the same internal climate oscillations be expected as those in the Holocene? Alternatively, if the Holocene TSI time series has captured permanent features of solar variability, should they be expected to occur embedded in the deep-time paleoclimate record? A solar activity proxy time series for the mid-Cenomanian is needed to answer these questions. However, the affinity of the M11 spectra with the Holocene TSI spectrum hints at a solar forcing origin. Both the MCE and oceanic anoxic event 2 (OAE2) are proposed to have occurred during 2.4 m.y. eccentricity minima (Mitchell et al., 2008). One possibility is that during a 2.4 m.y. eccentricity minimum, seasonality was not as pronounced, and solar cycles may have a relatively stronger expression in the depositional environment. Testing this hypothesis requires exploring similar variability in other records that have a wide temporal and spatial coverage.

Semi-precession-scale climate change is predicted for intertropical locations from 23.5°S to 23.5°N that experience two insolation maxima per year (Berger and Loutre, 1997; Berger et al., 2006). Semi-precession cycling was discovered in the equatorial Pacific thermocline (0–300 m) over the last glacial cycle (Jian et al., 2020) and in Cretaceous ocean deposits at a paleolatitude of 30°S (Park et al., 1993). These discoveries raise expectations for finding semi-precession signals in near-equatorial oceans throughout deep time. The strong 10,784 yr cycle in the PC2 redox proxy is one such candidate signal. The Iona core had a paleolatitude of 30°N, i.e., not intertropical; however, the southern KWIS gateway was connected to the intertropical Gulf of Mexico and beyond, so it was not immune to circulation originating from elsewhere. In fact, analysis of water-mass proxies in the Iona core indicate influence primarily from Tethyan oxygen minimum zone water masses arriving from the south (Arthur and Sageman, 2005; Eldrett et al., 2017). It is possible that these water masses carried a semi-precession signal with them from intertropical insolation forcing that affected dissolved oxygen.

Time-series analysis of multivariate paleoclimate proxy data obtained from a marine marlstone with interbedded dark/light laminae formed during the greenhouse climate just prior to and during the lower part of the mid-Cenomanian event (ca. 96.5 Ma) reveal significant (>99%) cycling with periods of 80–100 yr, 200–230 yr, and 400 yr reminiscent of Holocene solar activity cycles. Laminae deposition is interpreted to reflect variations in water-column stratification, redox conditions, and marine community structure influenced by solar activity–forced changes in storm activity and intensity. Such cycles may be unique to shallow seas, where bottom sediments are sensitive to storms. Principal component analysis of major elements measured by XRF core scanning further revealed a strong 10,784 yr cycle in Fe and S, taken to be a redox proxy that records the incursion of Tethyan oxygen minimum zone water sources.

This study was supported by the National Natural Science Foundation of China (grants 42172137, 41888101, 42050104 and 42050102), the Everest Scientific Research Program of Chengdu University of Technology (grant 2020ZF11402), and the U.S. National Science Foundation through grants EAR-1151438 (to S.R. Meyers) and OAC-1835717. We thank Ian Croudace for Itrax data. Finally, we thank David De Vleeschouwer, Lucas Lourens, Hong Ao, Nicolas Thibault, and four anonymous reviewers for their constructive comments.