Precession-driven climate cycles and time scale prior to the Hirnantian glacial maximum

The theory of astronomical climate forcing has revolutionized our understanding of Cenozoic climate systems and is the basis for unprecedented continuous time scales (astrochronologies) with precision down to ~10⁴ yr (Zachos et al., 2001). Pre-Cenozoic astrochronologies face several challenges relating to (1) uncertainties in the deep-time astronomical solutions and parameters (Berger and Loutre, 1994; Waltham, 2015); (2) less-complete and less-well-preserved strata; and (3) the sparsity of geochronologic anchor points. Consequently, Paleozoic astrochronologies are typically based on identification of the stable 405 k.y. eccentricity cycle instead of shorter astronomical cycles, which have the potential to provide an order-of-magnitude increase in temporal resolution. However, the prevalence of eccentricity-based astrochronologies is mechanistically difficult to explain because most of the insolation power lies in the obliquity and precession bands—not in the eccentricity.

Several researchers have interpreted the record of eccentricity (~100 and ~405 k.y.) and long obliquity (~1.2 m.y.) cycles in the Upper Ordovician reference outcrop sections of Anticosti Island, Québec, Canada (Fig. 1; Long, 2007; Elrick et al., 2013; Ghienne et al., 2014; Mauviel et al., 2020). However, extrapolating these interpreted accumulation rates for the relatively homogeneous upper Katian subsurface lithology results in total time spans of tens of millions of years, which is inconsistent with integrated stratigraphic constraints indicating an estimated duration of only 4–5 m.y. (Cooper and Sadler, 2012; McLaughlin et al., 2016). While acknowledging challenges to pre-Cenozoic astrochronologies, our study makes use of new subsurface records from recent drilling on Anticosti Island and explores several astrochronological scenarios within the available stratigraphic constraints. Pre-Hirnantian glacial buildup is well established (Vandenbroucke et al., 2010; Pohl et al., 2016), and correctly documenting Late Ordovician astronomical cycles is crucial for constructing high-resolution time scales and studying dynamic changes in climate and biodiversity (Saupe et al., 2020). In general, reliable identification of high-frequency Paleozoic astronomical cycles can pave the way toward a greatly enhanced understanding of the interplay between biotic evolution and environmental change.

The lower Paleozoic mixed siliciclastic-carbonate succession of Anticosti Island contains some of the thickest, most fossiliferous, well-exposed, and little-altered Upper Ordovician sections preserved on Earth (Fig. 1; Desrochers et al., 2010; Finnegan et al., 2011). The succession was deposited within a structural embayment along the eastern margin of Laurentia on the distal portion of a Taconic-Acadian foreland located in the paleo(sub-)tropics (10°S–25°S; Fig. 1B; Blakey, 2013). Our study focuses on the subsurface part of the upper Katian Vauréal Formation, consisting of predominantly gray, interbedded micrite, calcarenite, and marl deposited in a mid- to outer-shelf environment (Fig. S1 in the Supplemental Material¹; Long, 2007). The lithological variations of interest in this study are multimeter bed bundles, and not the centimeter- to decimeter-thick limestone-marl couplets that are potentially early diagenetic in origin (Fig. 1C; Nohl et al., 2019). Previous outcrop studies suggested that eccentricity is the dominant astronomical signal in the Vauréal Formation (Long, 2007; Elrick et al., 2013; Mauviel et al., 2020). Updated biostratigraphy together with new chemostratigraphy indicate that the Vauréal Formation belongs to Ka4, a stage slice estimated at a total duration of 4–5 m.y. (Fig. 1D; see the Supplemental Material; McLaughlin et al., 2016).

The La Loutre #1 core (LL1; 49°35′18″N, 63°38′14″W) was drilled by Consortium Hydrocarbures Anticosti ~10 km southwest of the well-studied New Associated Consolidated Paper (NACP) core (49°37′20″N, 63°26′18″W; Fig. 1A; McLaughlin et al., 2016). The two cores were correlated using the informal lithological units V2–V5 defined by McLaughlin et al. (2016) (see also Figs. 1A and 2; Fig. S2; see the Supplemental Material). The greater thickness of lithological units in the LLI core is attributed to separation by ~10 km along the south-directed depositional dip of the Anticosti Basin. The NACP core and the western outcrop sections were correlated based on their bulk δ¹³C_carb (carb—carbonate) records (Figs. 1A and 2; Mauviel and Desrochers, 2016; McLaughlin et al., 2016). Additionally, 472 new NACP samples were measured for bulk δ¹³C_carb and δ¹⁸O_carb to increase the resolution of the record (see the Supplemental Material). The LL1 potassium (⁴⁰K perent) record as measured by natural gamma-ray (NGR) logging was used as proxy for time-series analyses. This proxy reflects the multimeter cycles of carbonate versus clay lithology of interest and is continuously recorded with the highest available resolution (10 cm). Evolutive harmonic analysis (EHA) (Thomson, 1982) and TimeOpt, eTimeOpt, and timeOptTemplate (Meyers, 2015, 2019) analyses were done with “Astrochron” (https://cran.r-project.org/web/packages/astrochron/index.html; Meyers, 2014) in R (R Core Team, 2017). TimeOpt is a statistical optimization method that can simultaneously consider power spectra distributions and amplitude modulation patterns (Meyers, 2015, 2019). Accumulation rate reconstructions were also done in MATLAB using ACEv.1 and spectral moments (Sinnesael et al., 2016, 2018). Code and data are available in the Supplemental Material.

The initial target for analysis was the V4–V5 stratigraphic interval of the Vauréal Formation, given that deeper intervals contained multiple hardgrounds and sharp contacts (e.g., the V2–V3 limit), and given that the overlying Ellis Bay Formation shows a pronounced facies shift toward more calcareous shales and was not completely recorded by the LL1 logging (Fig. 2; McLaughlin et al., 2016). Most lithological variation in the V4–V5 interval of LL1 (1030–335 m) is characterized by 5–10-m-thick cycles of clay-rich and clay-poor beds (Figs. 3B and 3C; Fig. S3), with a gradual decrease in thickness up section. These cycles are less clear in the lowermost (1030–880 m) and uppermost (403–335 m) intervals of the record (Fig. S3). Therefore, our analysis focused on the 880–403 m interval, which shows clear uninterrupted regular cycles (Figs. 2 and 3). Parallel to the decrease in cycle thickness up section, the amplitude changes attenuate, suggesting a gradual decrease in accumulation rate. The NGR signatures also appear to be bundled (5–6) into larger ~25–50-m-thick packages, with more clearly developed 5–10 m cycles occurring in the middle of the bundles (Fig. 3E).

Integrating the biostratigraphic (conodont, graptolite, and chitinozoan) and chemostratigraphic (⁸⁷Sr/⁸⁶Sr and δ¹³C_carb isotope) constraints with the numerical Ordovician 2012 Geological Time Scale (Cooper and Sadler, 2012; see the Supplemental Material), McLaughlin et al. (2016) inferred that the entire ~1100-m-thick Vauréal Formation represents a few million years, implying that the 5–10 m cycles each represents a few tens of thousands of years. This interpretation is also consistent with the 2020 Ordovician Geological Time Scale (Goldman et al., 2020). These independent time constraints exclude eccentricity time scales for the 5–10 m cycles. However, they do not allow us to discriminate between precession and obliquity. Nonetheless, both precession and obliquity should be amplitude modulated by eccentricity and long-obliquity, respectively. A bundling of 5–6 cycles into larger cycles could be an indication of a Late Ordovician precession (~16–21 k.y.; Berger and Loutre, 1994; Waltham, 2015) and short eccentricity (~100 k.y.) ratio. However, a similar ratio exists between Late Ordovician obliquity (~30–33 k.y.) and the 173-k.y.-long obliquity cycle (Boulila et al., 2018). An independent astronomical characteristic that can distinguish between both bundling scenarios is the identification of the longer 100 and 405 k.y. eccentricity cycles in a precession-eccentricity–dominated record. This hypothesis can be tested statistically using TimeOpt.

The eTimeOpt analysis supports a precession-eccentricity–dominated signal with accumulation rates decreasing from ~60 to 30 cm k.y.^–1 (Fig. 3A). This trend agrees with visual inspection (manual tuning) and additional independent numerical accumulation rate reconstructions (Fig. S4). We used this trend to additionally constrain the timeOptTemplate analysis, statistically testing a precession-eccentricity signal and reconstructing precession and eccentricity in the time domain. The resulting precession band-pass filters and reconstructed eccentricity (combined fitting of the precession-band amplitude and a user-defined eccentricity model containing the 405, 125, and 95 k.y. periods) show clear precession amplitude modulations by the ~100 and 405 k.y. eccentricity periods (Fig. 3D). Changes in accumulation rates follow the 25–50 m (~100 k.y. eccentricity) cycles, with times of high eccentricity corresponding with higher accumulation rate (Fig. 3; Fig. S4). Such an interpretation suggests a positive phase relationship between the occurrence of thicker clay-rich layers and eccentricity maxima. Located at the southeast margin of the large Laurentian landmass at paleo(sub-) tropical latitudes (Fig. 1B), the Anticosti Basin could have been prone to monsoons. We suggest a monsoon-driven climatic control on the input of siliciclastic detrital material to explain the observed variations in clay content on astronomical time scales, as commonly documented throughout the Phanerozoic (cf. De Vleeschouwer et al., 2012; Wang et al., 2014).

We estimated the total duration of the 880–403 m interval in several ways. The timeOptTemplate approach, optimizing the precession-eccentricity relationship, resulted in a total duration of 1138 k.y. (Fig. S5), whereas the short-long eccentricity optimization gave a total duration of 1245 k.y. (Fig. S6). Manually counting 66 lithological cycles multiplied by an average precession duration of 18.725 k.y. for 445 m.y. ago (Waltham, 2015) led to a total duration of 1236 k.y. A possible explanation for this slight offset might be that some exceptionally thick cycles are not single precession cycles; instead, they are obliquity cycles that become more prominent during eccentricity minima when the amplitude of the climatic precession is relatively low (Fig. 3E; Fig. S3).

Overall, the studied 477-m-thick subsurface interval is interpreted to represent ~1.2 m.y., recording three 405 k.y. eccentricity cycles, 12 ~100 k.y. eccentricity cycles, and 60–66 precession cycles of ~18 k.y., with a potential imprint of obliquity cycles during eccentricity minima. Extrapolating the resulting accumulation rates across the entire Vauréal Formation suggests a represented total duration of ~3 m.y. This interpretation is in closer agreement with the available integrated stratigraphic constraints (McLaughlin et al., 2016) compared to the previous outcrop-based cyclostratigraphic interpretations assigning eccentricity-scale durations to comparable meter-scale lithological cycles.

Interestingly, the clear precession imprint disappears at ~400 m LL1 core depth toward the top of the Vauréal Formation (Fig. 2; Fig. S3). This interval coincides with a period of cooling (Fig. 2). Although caution is needed while interpreting Paleozoic bulk δ¹⁸O records, the trends in the Anticosti Island bulk δ¹⁸O record are consistent with results of studies using clumped oxygen isotopes and well-preserved brachiopod and conodont δ¹⁸O data (Finnegan et al., 2011; Goldberg et al., 2021), as well as with the averaged δ¹⁸O_conodont apatite data from Elrick et al. (2013) (Fig. 2). Late Katian cooling could have shifted the position of the Intertropical Convergence Zone and monsoon circulation so that the paleolocation of the Anticosti Basin was less influenced by monsoon climate forcing (Armstrong et al., 2009). Additionally, the increasingly large ice buildup on Gondwana could have shifted the dominant global astronomical signature toward stronger obliquity (Herrmann et al., 2003) or eccentricity imprints as suggested for the uppermost Vauréal Formation or Hirnantian Ellis Bay Formation (Long, 2007; Mauviel et al., 2020). Thus, the weakened precession signal might reflect a stronger nonlinear response of the climate system to the precession forcing, thereby enhancing the strength of the modulating eccentricity signal. In such a scenario, the δ¹⁸O_conodont apatite variations interpreted by Elrick et al. (2013) could reflect glacio-eustatic fluctuations.

Alternatively, the δ¹⁸O_conodont apatite variations could include more than glacio-eustatic fluctuations. Next to a control on supply of siliciclastic material, a strong monsoon circulation could have had an important influence on the continental runoff of freshwater and local evaporation rates. Simulations for the late Carboniferous North American Midcontinent epicontinental sea, a setting that could be compared with the Late Ordovician Laurentian epicontinental sea (Fig. 1B), demonstrate how such seawater freshening might considerably impact marine biogenic δ¹⁸O proxy interpretations (Macarewich et al., 2021). An additional possibility is a link with sub-Milankovitch climate variability. Our reconstructed accumulation rates result in some of the highest-time-resolution windows on the Ka4 slice worldwide. Spectral analyses of the precession-calibrated time series as shown in Figure 3E demonstrate a consistent cycle of ~1–1.5 k.y. (see the Supplemental Material), corresponding with the Pleistocene Dansgaard-Oeschger oscillation period (Dansgaard et al., 1989). Even though Paleozoic boundary conditions were fundamentally different compared to the Cenozoic (Elrick and Hinnov, 2007), evidence for a Paleozoic record of millennial-scale climate variability is growing (e.g., Da Silva et al., 2018).

We demonstrated how identification of the unique amplitude modulation characteristic of precession and eccentricity can be used as a powerful tool to develop Paleozoic ~10⁴-yr-resolution astrochronologies. Using such astronomical signal properties in combination with integrated stratigraphic constraints allowed us to move beyond the limits of stable 405 k.y. eccentricity astrochronologies, while recognizing the stratigraphic limitations imposed by the occurrences of unconformities or hiatuses and intervals where cycles might not be clearly expressed. Such an approach is complementary to longer, but less precise, cyclostratigraphic studies. Our unprecedented temporal resolution provides a snapshot into late Katian climate dynamics, paving the way for future “Cenozoic-style” scenario reconstructions.

We thank the Research Foundation–Flanders (M. Sinnesael: Ph.D. fellowship FWOTM782), the King Baudouin Foundation (Brussels) (T. Vandenbroucke and J. De Weirdt: Professor T. Van Autenboer Fund), Bijzonder Onderzoeksfonds—Universiteit Gent (BOF-UGent) (T. Vandenbroucke: BOF17/STA/013), the Natural Sciences and Engineering Research Council of Canada (A. Desrochers: NSERC Discovery Grant), and the Vrije Universiteit Brussel (VUB) Strategic Research Program (P. Claeys) for funding, and M. Elrick and F. Hilgen for constructive reviews. This work contributes to International Geoscience Programme projects IGCP 652 (Reading geologic time in Paleozoic sedimentary rocks) and IGCP 653 (The onset of the Great Ordovician Biodiversification Event).