Solar System chaos and the Paleocene–Eocene boundary age constrained by geology and astronomy

See allHide authors and affiliations

Science  30 Aug 2019:
Vol. 365, Issue 6456, pp. 926-929
DOI: 10.1126/science.aax0612

Filling a dating hole

The periodic nature of Earth's orbit around the Sun produces cycles of insolation reflected in climate records. Conversely, these climate records can be used to infer changes in the dynamics of the Solar System, which is inherently chaotic and not always similarly periodic. A particular obstacle is the lack of well-defined planetary orbital constraints between 50 and 60 million years ago. Zeebe and Lourens found an astronomical solution for that interval showing that the Solar System experienced a specific resonance transition pattern. These data provide a measure of the duration of the Paleocene-Eocene Thermal Maximum.

Science, this issue p. 926


Astronomical calculations reveal the Solar System’s dynamical evolution, including its chaoticity, and represent the backbone of cyclostratigraphy and astrochronology. An absolute, fully calibrated astronomical time scale has hitherto been hampered beyond ~50 million years before the present (Ma) because orbital calculations disagree before that age. Here, we present geologic data and a new astronomical solution (ZB18a) showing exceptional agreement from ~58 to 53 Ma. We provide a new absolute astrochronology up to 58 Ma and a new Paleocene–Eocene boundary age (56.01 ± 0.05 Ma). We show that the Paleocene–Eocene Thermal Maximum (PETM) onset occurred near a 405-thousand-year (kyr) eccentricity maximum, suggesting an orbital trigger. We also provide an independent PETM duration (170 ± 30 kyr) from onset to recovery inflection. Our astronomical solution requires a chaotic resonance transition at ~50 Ma in the Solar System’s fundamental frequencies.

Numerical solutions for the Solar System’s orbital motion provide Earth’s orbital parameters in the past, which are widely used to date geologic records and investigate Earth’s paleoclimate (111). The Solar System’s chaoticity imposes an apparently firm limit of ~50 million years before the present (Ma) on identifying a unique orbital solution, as small differences in initial conditions and parameters cause astronomical solutions to diverge around that age [Lyapunov time ~5 million years (Myr); supplementary materials] (4, 6, 12, 13). Recent evidence for a chaotic resonance transition (change in resonance pattern; see below) in the Cretaceous (Libsack record) (9) confirms the Solar System’s chaoticity but does not provide constraints to identify a unique astronomical solution. The unconstrained interval between the Libsack record (90 to 83 Ma) and 50 Ma is too large a gap, allowing chaos to drive the solutions apart (supplementary materials). Thus, proper geologic data around 60 to 50 Ma are essential to selecting a specific astronomical solution and, conversely, the astronomical solution is essential to extending the astronomically calibrated time scale beyond 50 Ma.

We analyzed color reflectance data (a*, red-to-green spectrum) (7, 8) at Ocean Drilling Program (ODP) Site 1262 (supplementary materials), a*-1262 hereafter, a proxy for changes in lithology (7). The related Fe-intensity proxy (8) gives nearly identical results (fig. S4). We focus on the section at ~170 to 110 m (~58 to 53 Ma), which exhibits an exceptional expression of eccentricity cycles at Site 1262 (7, 8, 10, 14, 15), less so in the preceding (older) section. Our focus interval includes the PETM and Eocene Thermal Maximum 2 (ETM2), extreme global warming events (hyperthermals), considered the best paleo-analogs for the climate response to anthropogenic carbon release (1618). The PETM’s trigger mechanism and duration remain highly debated (1921). Thus, in addition to geological and astronomical implications, unraveling the chronology of events in our studied interval is critical for understanding Earth’s past and future climate.

We developed a simple floating chronology, attempting to use a minimum number of assumptions (supplementary materials). We initially used a uniform sedimentation rate throughout the section (except for the PETM) and a root mean square deviation (RMSD) optimization routine to derive ages (for final age model and difference from previous work, see supplementary materials). No additional tuning, wiggle-matching, or preexisting age model was applied to the data. Using our floating chronology, the best fit between the filtered and normalized data target a** (Fig. 1) and a given astronomical solution was obtained through minimizing the RMSD between record and solution by shifting a** along the time axis (offset τ) over a time interval of ±200 thousand years (kyr), with ETM2 centered around 54 Ma (supplementary materials). Before applying the minimization, both a** and the solution were demeaned, linearly detrended, and normalized to their respective standard deviation (Fig. 1).

Fig. 1 Data analysis and comparison of color reflectance a* to our astronomical solution ZB18a.

(A) a* at ODP Site 1262 (blue-green), interpolated, demeaned, detrended record (Norm. a*) including PETM stretch (light-blue); scaled long/short eccentricity cycle filter (blue/gray), PETM onset (up-triangle), PETM recovery inflection (down-triangle), Elmo (square). mcd, meters composite depth. As primary CaCO3 variations within the PETM interval are not preserved due to dissolution and erosion, the interval was cropped. (B) Sum of long- and short-cycle filter outputs in the time domain (data target a**, light blue) and normalized eccentricity of Earth’s orbit from our astronomical solution ZB18a (purple). a** and ZB18a were demeaned, detrended, and normalized to their respective standard deviation before optimization (RMSD minimization between a** and solution by stretch-shift operation, see text). Across the cropped PETM interval, a** provides no actual information and is omitted. Up-triangle and error bar indicate our new age for the PEB (PETM onset) of 56.01 ± 0.05 Ma. Arrow indicates the prolonged high-eccentricity period before the PETM (see text).

It turned out that one additional step was necessary for a meaningful comparison between a** and astronomical solutions. Relative to all solutions tested here, a** was consistently offset (shifted toward the PETM after optimizing τ) by about one short eccentricity cycle for ages either younger (some solutions) or older than the PETM (other solutions). The consistent offset relative to the PETM suggests that the condensed PETM interval in the data record is the culprit, for which we applied a correction, also obtained through optimization. At Site 1262, the PETM is marked by a ~16-cm clay layer (<1 weight % CaCO3), largely due to dissolution and some erosion across the interval (16, 22), although erosion of Paleocene (pre-PETM) sediment alone cannot account for the offset of about one short eccentricity cycle (supplementary materials). Sedimentation rates were hence nonuniform across the PETM interval (8, 10, 16), and primary lithologic cycles from variations in CaCO3 content are not preserved within the clay layer. Thus, we corrected the condensed interval by stretching a total of k grid points across the PETM by Δz for a total length of ΔL = kΔz and included k as a second parameter in our optimization routine (Fig. 1). Essentially, the correction for the reduction (gap) in carbonate sedimentation across the PETM is determined by the entire record except the PETM itself (supplementary materials). In summary, we minimized the RMSD between data target and solution by a stretch-shift operation, i.e., we simultaneously optimized the number of stretched PETM grid points (k) and the overall time shift (τ) between floating chronology and solution.

Our new astronomical solution, ZB18a [computations build on our earlier work (6, 23, 24), supplementary materials], agrees exceptionally well with the final a** record (Fig. 1B) and has the lowest RMSD of all 18 solutions published to date that cover the interval (Table 1). The 18 solutions were computed by multiple investigators, representing different solution classes due to initial conditions, parameters, etc. (supplementary text S6 and S7). Based on ZB18a, we provide a new astronomically calibrated age model to 58 Ma (Fig. 1B and supplementary materials) and a revised age for the Paleocene–Eocene boundary (PEB) of 56.01 ± 0.05 Ma (see supplementary materials for errors). Our PEB age differs from previous ages (8, 2527) but is close to approximate estimates from 405-kyr cycle counting across the Paleocene (28) (supplementary materials).

Table 1

RMSD between a** record and selected astronomical solutions.†

View this table:

Contrary to current thinking (8, 14, 20, 27, 29), the PETM onset therefore occurred temporally near, not distant, to a 405-kyr maximum in Earth’s orbital eccentricity [Fig. 1, compare also (10)]. As for ETM2 and successive early Eocene hyperthermals (7. 29, 30), this suggests an orbital trigger for the PETM, given theoretical grounding and extensive, robust observational evidence for eccentricity controls on Earth’s climate (2, 710, 14, 20, 2632). Note, however, that the onset does not necessarily coincide with a 100-kyr eccentricity maximum (see below). Our analysis also provides an independent PETM main phase duration of 170 ± 30 kyr from onset to recovery inflection (for tie points, see fig. S6 and supplementary materials). This duration might be an underestimate given that sedimentation rates increased during the PETM recovery (compacting the recovery would require additional stretching of the main phase). Our duration is significantly longer than 94 kyr (20) but agrees with the 3He age model at Site 1266 (167±2434 kyr) (21) and is consistent with >8 cycles in Si/Fe ratios at Zumaia (31), which, we suggest, are full (not half) precession cycles.

If high orbital eccentricity (e) also contributed to the long PETM duration (e ≈ 0.025 to 0.044 during PETM), then the potential for prolonged future warming from eccentricity is reduced due to its currently low values (e ≈ 0.0024 to 0.0167 during next 100 kyr). A similar argument may hold for eccentricity-related PETM trigger mechanisms. The PETM occurred superimposed on a long-term, multimillion year warming trend (7, 30). Our solution ZB18a shows a 405-kyr eccentricity maximum around the PETM but reduced 100-kyr variability (Fig. 1B). Eccentricity in ZB18a remained high before the PETM for one short eccentricity cycle (Fig. 1B, arrow), suggesting that the combination of orbital configuration and background warming (30, 32) forced the Earth system across a threshold, resulting in the PETM. Although similar thresholds may exist in the modern Earth system, the current orbital configuration (lower e) and background climate (Quaternary/Holocene) are different from 56 Ma. None of the above, however, will directly mitigate future warming and is therefore no reason to downplay anthropogenic carbon emissions and climate change.

Our astronomical solution ZB18a shows a chaotic resonance transition (change in resonance pattern) (33) at ~50 Ma, visualized by wavelet analysis (34) of the classical variables:

h=e sinϖ;p=sin(I/2)sinΩ ,(1) where e, I, ϖ, and Ω are eccentricity, inclination, longitude of perihelion, and longitude of ascending node of Earth’s orbit, respectively (Fig. 2). The wavelet spectrum highlights several fundamental frequencies (g’s and s’s) of the Solar System, corresponding to eigenmodes. For example, g3 and g4 are loosely related to the perihelion precession of Earth’s and Mars’ orbits (s3 and s4 correspondingly to the nodes). The g’s and s’s are constant in quasiperiodic systems but vary over time in chaotic systems (supplementary materials). The period P43 associated with g4g3 switches from ~1.5 to ~2.4 Myr in ZB18a ~50 Ma, characteristic of a resonance transition (Fig. 2, arrow) (33). An independent analysis of the a*-1262 record recently also reconstructed P43 ≈ 1.5 Myr (35) within the interval ~56 to 54 Ma. However, our individual g-values from ZB18a are different from the reconstructed mean values, although within their 2σ error bounds (supplementary materials).

Fig. 2 Wavelet analysis of astronomical solution.

Wavelet analysis (34) of (A) h=e sinϖ and (B) p = sin(I/2) sin Ω from our astronomical solution ZB18a (see text). g’s and s’s indicate fundamental frequencies of the Solar System’s eigenmodes (multiple frequencies are expressed in the spectrum of a single planet). For example, g3 and g4 are loosely related to the perihelion precession of Earth’s and Mars’ orbits. The wavelet amplitude (red, peaks; blue, valleys) in, e.g., the g3 and g4 frequency band is modulated by the period 1/(g4g3) ≈ 2.4 Myr for ages younger ~45 Ma, where g3 ≈ 1/74.61 kyr–1 and g4 ≈ 1/72.33 kyr–1 (6). Correspondingly, 1/(s4s3) ≈ 1.2 Myr. However, in our solution, the period associated with g4g3 switches from ~1.5 Myr to ~2.4 Myr across the resonance transition around 50 Ma (arrow). The ratio (g4g3):(s4s3) ≈ 1:2 after ~45 Ma and closer to 1:1 before but appears irregular.

Notably, parameters required for long-term integrations compatible with geologic observations (e.g., ZB18a versus a**, Fig. 1) appear somewhat incompatible with our best knowledge of the current Solar System. For instance, ZB18a is part of a solution class featuring specific combinations of number of asteroids and solar quadrupole moment (J2), with J2 values lower than recent evidence suggests (supplementary materials). In addition, the La10c solution (33) with a small RMSD (Table 1) used the INPOP08 ephemeris, considered less accurate than the more recent INPOP10 used for La11 (13). However, La10c fits the geologic data better than La11 does [Table 1 and (27)].

The resonance transition in ZB18a is an unmistakable manifestation of chaos and is also key to distinguishing between different solutions before ~50 Ma, e.g., using the g4g3 term. This term modulates the amplitude of eccentricity and, e.g., the interval between consecutive minima in a 2-Myr filter of eccentricity (Fig. 3). Other solutions such as La10c (33) also show a resonance transition around 50 Ma. However, the pattern for ZB18a is different before 55 Ma, which is critical for its better fit with the data record from 58 to 53 Ma (smaller RMSD; Table 1 and Fig. 1). For example, P43 ≈ 2 and ~1.6 Myr at ~59 and ~56 Ma in La10c but is rather stable at ~1.5 to 1.6 Myr across this interval in ZB18a. Briefly, to explain the geologic record, our astronomical solution requires that the Solar System is (i) chaotic and (ii) underwent a specific resonance transition pattern between ~60 and 50 Myr BP.

Fig. 3 Resonance transition in selected astronomical solutions.

Interval between consecutive minima (Δtmin) in 2-Myr Gaussian filter (±60%) of Earth’s orbital eccentricity for selected solutions (6, 33). The rise of Δtmin around 50 Ma in ZB18a and La10c indicate resonance transitions. However, note distinct pattern of ZB18a before 55 Ma. Hence our solution ZB18a (closest agreement with the data record, Fig. 1) requires that the Solar System underwent a specific chaotic resonance transition pattern between ~60 and 50 Myr BP.

Supplementary Materials

Materials and Methods

Supplementary Text

Figs. S1 to S8

Tables S1 to S3

References (3659)

References and Notes

Acknowledgments: Funding: R.E.Z. acknowledges support from the U.S. National Science Foundation (OCE16-58023). This research was supported by the Netherlands Organisation for Scientific Research (NWO-ALW 865.10.001) and by the Netherlands Earth System Science Centre (NESSC 024.002.001) to L.J.L. We are grateful to the International Ocean Discovery Program for drilling and exploring ODP Site 1262. Author contributions: R.E.Z. and L.J.L. devised the study and wrote the manuscript. L.J.L. was instrumental in guiding the cyclostratigraphic analyses and selecting the a*-1262 record. R.E.Z. led the numerical integrations for the orbital motion of the Solar System. Competing interests: The authors declare no competing interests. Data and materials availability: Solutions for Earth’s orbital eccentricity and inclination are available at and We provide results from 100 to 0 Ma but caution that the interval 100 to 58 Ma is unconstrained due to chaos.
View Abstract

Stay Connected to Science

Navigate This Article