Geologic settings

The Kowala section (50.79° N, 20.57° E) was deposited in the intrashelf Chêciny-Zbrza Basin. The F–F interval at Kowala consists of thin-bedded marls with intercalated limestone beds and cherty limestone strata with marly intercalations. The LKW and UKW horizons are marked by two increases in total organic carbon (TOC) content, but the two black shale horizons that typify most F–F sections in Central Europe are absent here19. The Kowala section holds a previously published astrochronology, which is based on the stacking pattern of the marl-limestone alternations16. Both Kellwasser-equivalent TOC maxima, as well as four younger black shale horizons in the Famennian part of the section, were each deposited during just a fraction of one 405-kyr cycle according to that chronology. This finding was later confirmed for the Hangenberg shale, through high-precision U–Pb dating54. For the F–F interval, high-resolution δ13C carb (N = 77) and low-resolution δ13C carb data (N = 34) are available from De Vleeschouwer et al.16 and low-resolution δ13C carb (N = 38) data are available from Joachimski et al.19. To our knowledge, no magnetic susceptibility series exists for the F–F interval of the Kowala section.

Section C (53.17° N, 118.18° W) is located along the southeast margin of the Ancient Wall platform in the western part of the Western Canada Sedimentary Basin55. The section records upper-slope and prograding outer-ramp facies. The base of Section C is late Frasnian in age and consists of clay- and locally silt-rich shales of the Mount Hawk Formation, which is characterized by high magnetic susceptibility values throughout the basin18,55. A sea-level deepening event in the very late Frasnian resulted in the progradation of ramp facies and low magnetic susceptibility values at Section C, reflecting the increased carbonate input. The early Famennian Sassenach Formation records an influx of siliciclastics into the basin and related higher magnetic susceptibility values. From Section C, 50-cm spaced magnetic susceptibility data (N = 365)18 and δ13C org (N = 52)17 are available.

During the Late Devonian, the Fuhe section (24.99° N, 110.41° E) was located in a tectonically active basin, surrounded by carbonate platforms. Pelagic nodular limestones intercalated with calciturbidites and debris flow deposits dominate the F–F sequence. Chen and Tucker56 interpreted lithological cycles in the Fuhe section as the expression of precession cycles, stacked into meter-scale cycle-sets interpreted as 100-kyr eccentricity cycles. The UKW event is clearly marked by a +3‰ excursion in the high-resolution 10-cm spaced δ13C carb (N = 54) and δ13C org (N = 52) records3. The LKW event occurs in the lower part of the section, in the upper rhenana conodont zone, and is characterized by a positive excursion in the low-resolution δ13C carb (N = 73) and δ13C org (N = 54)20 (Fig. 1).

The Sinsin roadcut section (50.28° N, 5.24° E) consists of shales, marls, and (nodular) limestones, deposited on the continental shelf below fair-weather wave base57. In the framework of this study, we sampled the section at 5-cm resolution for magnetic susceptibility (N = 223) and at 10-cm resolution for δ13C carb (N = 116). The carbon isotope chemostratigraphy exhibits a gradual shift towards more positive δ13C carb values (Fig. 1). This +3‰ excursion concurs with the upper linguiformis conodont Zone58 and is therefore interpreted as the expression of the UKW. In the same interval as the positive δ13C carb excursion, TOC values peak up to 0.2 wt% in a 30-cm thick interval directly underlying the F–F boundary59. A previously published cyclostratigraphy does not exist for the Sinsin section. Kaiho et al.59 demonstrated the imprint of astronomical forcing in the Sinsin section through wavelet analysis, but did not delineate individual cycles.

The CG-1 core (42.93° N, 93.08° W) was deposited in a shelf margin environment and consists of fossiliferous shales and argillaceous limestones of the Lime Creek Formation. In the framework of this study, we sampled CG-1 at an average resolution of 11 cm for magnetic susceptibility (N = 329) and 26 cm for δ13C carb (N = 137). The H-32 core (40.47° N, 91.47° W) is the basinal, deep-water time-equivalent section of CG-1 and consists of the dark-colored Sweetland Creek shale. The overlying Grassy Creek shale contains seven layers of volcanic ash, none of which yielded datable zircons. We sampled H-32 at 3 or 4-cm intervals for magnetic susceptibility (N = 294) and at irregular intervals for δ13C carb (N = 165) and δ13C org (N = 70). The carbon isotope chemostratigraphy for both sections exhibits distinct positive δ13C excursions (Fig. 1). The +3‰ excursion in the CG1 core begins within the upper part of Frasnian Zone 12. We thus interpret this excursion as the expression of the LKW. At H-32, a +2‰ δ13C org excursion spans Zones 13b and 13c. We interpret this excursion as the UKW as the F–F boundary occurs at the base of the first ash bed just below 1000 cm. Between 800 and 900 cm, a smaller positive excursion is observed in both δ13C carb and δ13C org , which we interpret as the LKW. The H-32 and Cerro Gordo Project Hole 1 (CG-1) drill cores have never been studied in terms of astronomical forcing before.

Magnetic susceptibility

Magnetic susceptiblity measurements for the Fuhe and Sinsin sections were taken with a KLY-3S kappabridge at the University of Liège (Belgium). Samples from the CG-1 and H-32 cores were measured at the University of Wisconsin at Milwaukee (USA) with a MFK1 kappabridge. Both devices have a noise level <2 × 10−8 SI. Each magnetic susceptibility data point is the average of three measurements and has been corrected for sample mass. Samples typically have a mass larger than 10 g and were weighed with a precision of 0.01 g.

Carbon isotope stratigraphy

Bulk carbonate δ13C from the H-32, CG-1 and Sinsin sections and δ13C org from H-32 were measured at the Vrije Universiteit Brussel (Belgium). δ13C carb was measured using a NuPerspective IRMS interfaced with a NuCarb automated carbonate device. We use the internal NCM standard calibrated against NBS19 (+2.09‰). The long-term standard deviation of the internal standard is 0.05‰ for δ13C carb . Powdered samples for δ13C org were weighed (~10 mg) in silver capsules. The samples were decarbonated by adding 5% HCl in steps of 4 h and reaction at 50 °C, till no further reaction occurred. The δ13C org was measured using a Thermo1112 flash elemental analyzer coupled via a Con-Flo III interface to a Thermo DeltaV isotope ratio mass spectrometer (IRMS). The international reference standard IAEA-CH-6 (−10.449‰) was used for calibration (1σ ≤ 0.2‰).

Conodont faunas

The H-32 and CG-1 zonal nomenclature, presented in this study, follows the newest Frasnian60 and Famennian61 standard zonations.

We took 118 samples for biostratigraphy throughout the CG-1 core. All samples are typical of the Frasnian Polygnathus biofacies. Samples from the lower 23 m of the Juniper Hill and lower Cerro Gordo members include Palmatolepis semichatovae and Pa. kireevae correlated with Frasnian Zone 1162,63. The first appearance of Polygnathus planarius is interpreted as marking the base of Zone 12 based on its known range in western Canada. The occurrence of Polygnathus samueli high in the Cerro Gordo Member, just below the onset of the LKW δ13C excursion, is high within Frasnian Zone 12. The joint occurrence of Ancyrognathus asymmetricus and Icriodus alternatus in the Owen Member is interpreted as Frasnian Zone 13a, as the latter first occurs very high in Zone 12 and ranges into Zones 13a and 13b. Stratigraphically, this level occurs above peak values of the LKW δ13C excursion.

We took 54 samples for biostratigraphy throughout Core H-32. Most of the Sweetland Creek Shale at H-32 correlates with Frasnian Subzones 13a to 13b. The overlying Grassy Creek shale, just below the F–F boundary, corresponds to Frasnian Subzones 13b–13c (i.e., the linguiformis interval). The conodont yields of H-32 samples are generally very small or barren, but nevertheless, clear zonal indicators allow for the construction of the biostratigraphic scheme presented in Fig. 1. The Zone 11 and Zone 12 interval in the Sweetland Creek Shale is condensed in less than 3 meters. Most of the Sweetland Creek Shale at H-32 represents an expanded record of the LKW through the very Early Famennian.

Age modeling and spectral analyses

Tie-points between H-32, CG-1, Sinsin, Fuhe, Kowala, and Section C are obtained by visually correlating distinct features in magnetic susceptibility (red ties in Fig. 5) and carbon isotope geochemistry (blue ties in Fig. 6). In a first step, we assign an age relative to the F–F boundary to each of these tie-points according to the 405-kyr astrochronologic framework of Section C (Figure 7 in De Vleeschouwer et al.15). These ages are listed in black in Figs. 5 and 6. In a second step, we apply a Monte Carlo procedure to distort time-differences between consecutive tie-points. The goal of this procedure is to slightly modify the tie-point ages, so to get a better expression of the astronomical frequencies in the power spectra of the F–F proxy-series. We run the Monte Carlo procedure 5000 times in our search for the best astronomical fit. The tie-point ages for the best-fitting Monte Carlo simulation are listed in green in Figs. 5 and 6. The fit of each Monte Carlo simulation is evaluated as follows: we compute the power spectra for all susceptibility time series and consider the spectral power in the 100-kyr eccentricity (0.008–0.013 cycles/kyr), obliquity (0.025–0.035 cycles/kyr) and precession (0.045–0.065 cycles/kyr) bands. A good astronomical fit exists when high spectral power occurs in each of these bands. We measure this fit by subtracting the maximum AR1 confidence level occurring within each band from 100%. Low percentages thus represent high spectral power within the considered frequency band. Additionally, we measure the misfit between the frequency of maximum AR1 confidence level, and the expected orbital frequency according to Berger et al.26 (at 375 Ma: 0.0105 per kyr for 100-kyr eccentricity, 0.031 per kyr for obliquity and 0.055 per kyr for precession). This frequency misfit is also expressed as a percentage. Finally, we measure the total astronomical fit of the Monte Carlo simulation in question by averaging the spectral power and misfit measures among all frequency bands and among all sections. Our age modeling strategy is illustrated in detail using three distorted eccentricity-tilt-precession (ETP) series in Supplementary Figs. 3–5.

All spectral analyses in this study were carried out using the MTM with three 2π -tapers64, as implemented in the R-package “astrochron”65. The confidence levels were calculated applying the conventional AR1 method for estimating the red noise spectrum66. The R-code used for age modeling and spectral analyses is available through GitHub (“Code availability”).

Code availability

Software S1. R-script for the illustration of the age modeling strategy adopted in this study, using distorted eccentricity-tilt-precession (ETP) series is available at https://github.com/dadevlee/Late_Devonian. Software S2. R-script for age modeling of Late Devonian data is available at https://github.com/dadevlee/Late_Devonian. Software S3. R-script for spectral analyses and the generation of the data in Figs. 2, 3 and 4 is available at https://github.com/dadevlee/Late_Devonian

Data availability

Magnetic susceptibility and carbon isotope series are available from https://github.com/dadevlee/Late_Devonian/tree/master/Data_depth_domain and https://doi.pangaea.de/10.1594/PANGAEA.882366