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

Paleozoic astrochronologies are limited by uncertainties in past astronomical configurations and the availability of complete stratigraphic sections with precise, independent age control. We show it is possible to reconstruct a robust Paleozoic ∼ 10 4 -yr-resolution astrochronology in the well-preserved and thick Upper Ordovician reference record of Anticosti Island (Canada). The clear imprint of astronomical cycles, including ∼ 18 k.y. precession, potential obliquity, and short and long eccentricity, constrains the entire Vauréal Formation ( ∼ 1 km thick) to only ∼ 3 m.y. in total, representing ∼ 10 times higher accumulation rates than previously suggested. This ∼ 10 4 yr resolution represents an order of magnitude increase in the current standard temporal resolution for the Katian and even allows for the detection of sub-Milankovitch climate-scale variability. The loss of a clear precession signal in the uppermost Vauréal Formation might be related to contemporaneous global cooling prior to the Hirnantian glacial maximum as indicated by the δ 18 O record. Complementary to the study of cyclostratigraphy of longer and often simplified records, it is important to recognize stratigraphic hiatuses and complexities on the ∼ 10 4 yr scale to achieve robust sub-eccentricity-scale Paleozoic astrochronologies.


INTRODUCTION
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 4 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 orderof-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 Anti-costi 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.

GEOLOGICAL SETTING AND STRATIGRAPHY
The lower Paleozoic mixed siliciclastic-carbonate succession of Anticosti Island contains some of the thickest, most fossiliferous, wellexposed, 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 Material1 ; 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).

CYCLOSTRATIGRAPHIC ANALYSIS
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 ( 87 Sr/ 86 Sr and δ 13 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-eccentricitydominated 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. - (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 precessionband 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-) 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 timeOpt-Template approach, optimizing the precessioneccentricity 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).

IMPLICATIONS
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 δ 18 O records, the trends in the Anticosti Island bulk δ 18 O record are consistent with results of studies using clumped oxygen isotopes and well-preserved brachiopod and conodont δ 18 O data (Finnegan et al., 2011;Goldberg et al., 2021), as well as with the averaged δ 18 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 δ 18 O conodont apatite variations interpreted by Elrick et al. (2013) could reflect glacio-eustatic fluctuations.
Alternatively, the δ 18 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 δ 18 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).

CONCLUSIONS
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 4 -yrresolution 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 "Cenozoicstyle" scenario reconstructions.
Figure 1.(A) Map of the bedrock geology of Anticosti Island, Québec, Canada (after Mauviel et al., 2020).(B) Map of the Late Ordovician paleogeographic reconstruction for southern Laurentia (modified from Blakey, 2013).(C) Photo of differential weathering profile between more carbonate versus shale intervals of the upper Vauréal Formation in Vauréal Canyon representing precession (P) cycles bundled in short eccentricity (SE) cycles (person next to red bar for scale).(D) Stratigraphic column of carbonate-dominated Anticosti Island subsurface (modified from McLaughlin et al., 2016).

Figure 2 .Figure 3 .
Figure 2. Compiled stratigraphy and proxy data for the Upper Ordovician section of Anticosti Island, Québec, Canada.Dashed lines represent boundaries between informal lithological units as defined by McLaughlin et al. (2016).The V4-V5 boundary is mainly characterized by a decrease in clay content with more pure limestones in unit V5 (thinner dashed line).Red rectangle (880-403 m) indicates the interval selected for detailed time-series analysis (cf.Fig. 3).SMOW-standard mean ocean water; EB-Ellis Bay; V-Vauréal; K pXRF-potassium measured by portable X-ray fluorescence.