Climate Formula

Earth’s climate-relevant orbital forcing arises from the gravitational interplay of all eight planets. Their orbital and rotational cycles synchronise on the Solar System Resonance Cycle of 8H = 2,682,536 years, so every climate-relevant cycle lands at an integer divisor of 8H. The canonical climate formula is a three-layer decomposition: L1 = orbital forcing on the 8H integer-divisor lattice (32 components); L2 = climate-internal periodic response (the 405-kyr carbon-cycle thermostat and its harmonics); L3 = Heaviside step transitions at major Cenozoic boundary-condition changes (PETM, EOT, Mi-1, MMCT, iNHG, MPT). Sequential ridge regression fits the formula independently per climate regime, validated on four independent proxy records.

Climate cluster (companion pages): Climate Summary (synthesis) · L1 Attribution (per-integer dual attribution) · Insolation Null Test (ΔR² ≈ 0 evidence) · Related Work (peer-reviewed lit comparison). Reproducibility scripts in the 3d repository .

Canonical 3-layer Climate Formula on LR04 δ¹⁸O over the post-MPT window. Red = formula (L1 + L2 + L3, sequential ridge fit); black = observed LR04. R² = 0.87 within this regime.

Headline findings

1. R² up to 0.87 post-MPT LR04 — full per-regime table in §Sequential ridge regression.
2. Every significant LR04 peak lands on an integer divisor of 8H. 32 L1 integers carry the signal; 0 of 32 fully agree with Berger on both planet and mechanism — see L1 Attribution.
3. Mars and Jupiter co-dominate the per-planet fingerprint. Jupiter: 4 direct L1 matches. Mars: 3. Detail in §Per-Planet Contributions.
4. The 100-kyr glacial cycle is an inclination-side eigenmode beat, not direct eccentricity forcing. Empirical centroid at the s₁−s₄ nodal eigenmode beat at n = 25 = 107.3 kyr. Detail in §The 100-kyr Glacial Cycle.
5. Pre-MPT and post-MPT differ in climate sensitivity, not orbital forcing. Per-regime R²: pre-iNHG 0.43, iNHG-MPT 0.72, post-MPT 0.87.
6. Frequency-resolved ice-albedo decomposition. Ice-albedo share of total forcing: 57 % at the 100-kyr band, 64 % at obliquity, 69 % at precession, 46 % at long periods (>130 kyr). The +31 pp obliquity ice-share shift across the MPT (95 % pre-MPT → 64 % post-MPT) quantitatively confirms Willeit et al. 2019 . Detail in §Ice-Albedo Decomposition.
7. Forward projection: next natural glaciation peak ~60,500 AD (58,500 yr from now); strongest in the next quarter-million years at ~198,500 AD; warmest interglacial at ~166,000 AD. Detail in §Forward Projection.

The canonical three-layer climate formula

Architecture

$C(t) \;=\; c_0 \;+\; \underbrace{\sum_{n \in N_{L1}} \left[ a_n \cos\!\left(\tfrac{2\pi n t}{8H}\right) + b_n \sin\!\left(\tfrac{2\pi n t}{8H}\right) \right]}_{\text{L1: orbital forcing}} \;+\; \underbrace{\sum_{p \in P_{L2}} \left[ \alpha_p \cos\!\left(\tfrac{2\pi t}{p}\right) + \beta_p \sin\!\left(\tfrac{2\pi t}{p}\right) \right]}_{\text{L2: carbon thermostat}} \;+\; \underbrace{\sum_{i \in T_{L3}} \gamma_i \cdot H(t - t_i)}_{\text{L3: boundary-condition steps}}$

• L1 = 32 integers (25 canonical Berger / Mars–Jupiter beats + 6 precession-band sidebands + 1 Berger-quintet completion at n=141). See L1 Attribution for the per-integer dual attribution.
• L2 = three periods: 404.5 kyr (silicate-weathering thermostat fundamental), 202.25 kyr (2nd harmonic), 134.83 kyr (3rd harmonic) — confirmed carbon-amplified across LR04, CENOGRID, EPICA.
• L3 = up to 6 Heaviside transitions: PETM (56 Ma), EOT (34 Ma), Mi-1 (23 Ma), MMCT (14 Ma), iNHG (2.7 Ma), MPT (1.0 Ma), applied only when inside the fit window.
• t = age in kyr BP (positive = past, negative = future, t = 0 ≈ 2000 AD); 8H = 2,682,536 years; C(t) = normalised δ¹⁸O / δ¹³C / CO₂.

Sequential ridge regression, fit per regime

Coefficients fit by sequential ridge regression: L1 first (ridge λ = 1), then L2 on the residual, then L3 on the L1+L2 residual. Sequential fitting resolves the L1↔L2 collinearity at short windows; ridge resolves the L1↔L1 collinearity within the lattice (post-MPT condition number ~632, max VIF ~7×10⁴ under plain OLS). The formula is fit independently per climate regime — lattice positions are the same across regimes; per-line amplitudes vary.

At Quaternary scales (post-MPT LR04, EPICA CO₂) L1 carries the heavy lifting (~83–87% R² from the lattice alone). At deep-time scales (CENOGRID, CenCO2PIP) L3 step transitions dominate (~60% of variance over 67 Myr) — lattice oscillations average toward zero across many cycles, but discrete step-changes at major transitions carry the deep-time signal.

The same 8H lattice fitted to atmospheric CO₂ (EPICA Dome C, Bereiter 2015) — independent proxy, independent chronology, entirely different recording physics from benthic δ¹⁸O. R² = 0.83 from L1 alone, 0.85 with L2 added: the orbital structure that drives ice volume also drives atmospheric carbon-cycle dynamics.

What survived investigation

The formula commits to 32 + 3 + 6 = 41 structural components. Excluded candidates:

• 13H = 4.36 Myr Boulila libration — failed cross-window stability (amplitude CV 42–50%; phase circular std 97.9° ≈ uniform); not deployed.
• 9-Myr long-period carbon resonance — ΔR² CENOGRID δ¹³C = +0.078 with carbon-amplification ratio 4.2; investigated candidate, cross-window stability test pending.
• 405-kyr 4th and 5th harmonics (101 and 81 kyr) — collinear with L1 lattice integers 8H/27 and 8H/33; absorbed by L1.
• 23 off-lattice Laskar eigenmode beats — all fail the promotion criteria (ΔR² > 0.005 + carbon-amplification > 0.5).
• Piecewise / polynomial / LOESS detrending — superseded by L3 step covariates (5–10× better R² on CENOGRID).
• Linear-response ODE models for L2 — failed (all L2 ↔ L1 amplitude correlations |r| < 0.4); L2 dynamics are nonlinear (threshold / hysteresis / saturating).

Full inventory: doc 92 §9 .

The 32 L1 lattice integers (summary)

Each n corresponds to a celestial-mechanics quantity of one of two kinds: (i) an eigenmode beat (Laskar 2004 secular eigenfrequency combination — composite, multi-planet; planet labels follow Berger’s convention but the model treats them as math objects, see Eigenfrequencies); or (ii) a direct planet cycle from the 8H period table. The dominant peaks and their interpretations:

The full per-integer table with all 32 components, ranked amplitudes, dual attribution (Berger eigenmode vs Holistic top-1 multi-planet beat), and the precession-band sideband completions is the canonical reference at L1 Attribution.

Several integers correspond directly to specific planet cycles from the 8H period table: n=9 (Mercury Cassini lock), n=16 (Mars Axial / Jupiter Obliquity), n=21 (Mars Obliquity / Jupiter Axial — the “Mars–Jupiter Axial-Obliquity Swap”), n=39 (Jupiter ecliptic perihelion), n=65 (Saturn ecliptic = Jupiter ICRF perihelion — Law 6), n=68 (Mars ICRF perihelion), n=110 (Venus ICRF = Venus Obliquity). All 32 land on integer divisors of 8H because 8H is the natural synchronisation period of the solar system.

Closure: no orphan peaks off the 8H lattice

The strongest possible test of the 8H framework is its closure: do all significant spectral peaks in LR04 land on integer divisors of 8H, or are there orphan peaks at non-integer positions that would imply forcing from outside the planetary-eigenmode framework?

Closure test: fit all 200 integer divisors of 8H jointly to LR04 (joint OLS, R² = 0.443), then scan the residual for non-integer power. Result:

• Residual amplitude at integer positions: 0.000 (machine zero, by construction).
• Noise floor at random non-integer positions: 0.029 median, 0.127 (95th percentile).
• 14 residual peaks above noise threshold sitting >0.3 from any integer — every one between two adjacent integer divisors.

The biggest “orphan” (n = 65.45, period 40.99 kyr, amp 0.544) sits between n=65 (k+s₃ obliquity, 41.27 kyr) and n=66 — the arithmetic-mean position of the non-stationary obliquity band (Jensen’s-inequality phenomenon). The other 13 orphans similarly sit between adjacent integer beats — the expected fingerprint of cycle-length dispersion and spectral leakage. Zero orphan peaks land in “empty” regions of the 8H lattice. The 8H integer-divisor framework captures the full frequency structure of LR04’s significant spectral content.

Framework validation and the 405-kyr investigation

Beyond the closure test, the framework has been subjected to a battery of independent falsifiable tests. Most predictions hold; one specific signal — the 405-kyr long-eccentricity cycle — does not fit the orbital lattice and is best explained as a climate-system internal phenomenon.

Super-cycle hypothesis on deep-time geology: NULL. A pre-registered test of whether 20 major Phanerozoic geological/climatic events cluster preferentially at integer multiples of 8H returned p = 0.233 (primary 8H) and p = 0.504 (sharpened H). The Plio-Pleistocene 1×8H + 1×8H pattern is most parsimoniously a climate-response amplification artifact combined with statistical coincidence, not a deep-time geological clock. The framework describes orbital geometry that paces possible climate transitions, not a structural clock that causes geological events.

Fourteen follow-up tests on framework predictions (chronology validation, statistical baselines, bispectral coupling, external-publication convergence, cross-validation, phase prediction, cross-proxy replication, 67-Myr Cenozoic generalization, per-line Thomson MTM F-test significance, time-frequency centroid stability, head-to-head against Laskar 2004 placements) returned 16 clean positives, 2 partials, 5 nulls. Each null traces cleanly to one of four mechanistic causes (Rayleigh resolution at T < 8H, MPT amplitude non-stationarity, proxy-specific spectral concentration, non-linear climate lags); none falsifies the framework. Full per-test reproduction in doc 91 .

The 405-kyr cycle is L2, not L1

Standard Milankovitch attributes the 405-kyr line to the Laskar g₂−g₅ “Venus-Jupiter apsidal beat”. In this framework’s planet motions, Venus and Jupiter have different apsidal periods than Laskar’s eigenfrequencies (Venus: 8H/6 = 447 kyr retrograde; Jupiter: 8H/39 = 68.8 kyr; their actual beat is ~58–79 kyr, not 405 kyr). So the standard attribution does not carry over.

Mathematically: 405 kyr is unreachable on the 8H lattice. A systematic search across all pair and triplet beats among the 46 cycles in the 8H period table finds zero combinations within ±3% of 405 kyr. The closest matches cluster at 8H/7 = 383 kyr (5.4% below) and 8H/6 = 447 kyr (10.4% above).

Carbon-cycle signature. CENOGRID full Cenozoic places the line at 404.5 kyr with FWHM = Rayleigh resolution (clean single-frequency line). The δ¹³C/δ¹⁸O amplitude ratio is 1.53× — vs 0.49 at obliquity and 0.66 at precession. The signal lives preferentially in the carbon record. Phase residuals from a linear orbital model correlate at only r = 0.21 between δ¹³C and δ¹⁸O; about 96% of phase variance evolves independently between the two proxies.

The 405-kyr cycle behaves as an entrained internal oscillator: the silicate-weathering thermostat has a narrow resonance peak near 400 kyr (Walker, Hays & Kasting 1981; Pälike et al. 2006 “heartbeat of the Oligocene climate system”); long-period orbital eccentricity supplies energy that synchronises but does not rigidly drive it. The L2 layer captures it with 405 / 202 / 135 kyr (fundamental + 2nd + 3rd harmonics). The 9-Myr long-period candidate is investigated but not yet deployed; sub-harmonic (“17 × precession”) alternatives are ruled out by three independent tests, and the g₄−g₃ 2.4-Myr beat shows no comparable carbon-amplification (ratio 0.20 vs 1.53), refining the interpretation: the silicate-weathering thermostat has a narrow resonance peak near 400 kyr, not a broad long-period amplification.

Per-Planet Contributions

Cross-referencing each L1 lattice integer against the full 8H period table reveals which planets contribute directly versus via eigenmode beats. The per-integer dual attribution is in L1 Attribution; this section is the aggregate per-planet view.

Direct L1 matches per planet

Mars and Jupiter co-dominate. Mars also participates in three additional climate peaks via eigenmode beats: n=25 (s₁−s₄ nodal, 100-kyr centroid), n=28 (g₄−g₅ eccentricity, Berger 95k), and n=53 (Mars-Uranus s-beat near Mars eccentricity at 8H/52). The Mars-Jupiter resonance lock at 8H/36 (Mars.Peri_ecl = Jupiter.AscNode = 74.5 kyr) is the cleanest planet-pair resonance in the framework, sitting one Rayleigh element off the LR04 n=35 peak at 76.6 kyr.

Neptune and Uranus appear only pre-MPT

The pre-MPT window (1,200–3,000 kyr BP) reveals three additional peaks involving outer planets: n=14 (g₂−g₈ Venus-Neptune), n=30 (g₃−g₇ Earth-Uranus), n=38 (s₈−s₃ Neptune-Earth). Neptune’s eigenmodes are the slowest (g₈ ≈ 0.67″/yr, |s₈| ≈ 0.69″/yr); beats with inner planets produce ~70–200 kyr periods with small amplitudes. Pre-MPT, ice sheets were smaller and less hysteretic — slow outer-planet beats registered above noise. Post-MPT, large ice sheets impose a ~100-kyr hysteretic response that dominates and filters out slower outer-planet signals.

Cross-planet obliquity validation

Mercury’s 0.09% match against the independent Cassini-state forced-obliquity calculation (Bills & Comstock used Cassini-state dynamical theory, corroborated by Yseboodt & Margot 2006, Peale 2005, Bois & Rambaux 2007) is the model’s tightest cross-validation against non-Holistic references.

The 100-kyr Glacial Cycle

The 100-kyr glacial cycle is a multi-planet eigenmode-beat band, not direct eccentricity forcing. LR04’s 100-kyr band is a broad single peak spanning ~80–125 kyr — not the split structure (95k + 125k) that direct eccentricity forcing would produce. Energy is spread across three adjacent lattice integers:

• n = 25 = 107.3 kyr — s₁ − s₄ nodal eigenmode beat (energy-weighted centroid; a nodal / orbital-plane beat between two inner-rocky-planet eigenmodes)
• n = 28 = 95.8 kyr — g₄ − g₅ eccentricity beat (Berger’s classical 95-kyr peak; highest amplitude in the post-MPT spectrum)
• n = 22 ≈ 122 kyr — s₂ − s₄ nodal eigenmode beat (smaller contribution)

The centroid identification is consistent in spirit with Muller & MacDonald 1997 ’s “inclination, not eccentricity” framing, with a specific eigenmode identification not previously stated. Berger 1978’s single-planet attribution names s₁−s₄ “Mercury-Mars”; the Holistic model accepts the eigenmodes as math objects but does not endorse the single-planet attribution — see L1 Attribution.

Earth’s H/3 does not drive climate directly

Earth’s own inclination precession (H/3 = ~111,772 years = 111.77 kyr) is a real integer-divisor position on the 8H lattice (n = 24) but does not play a direct climate-driving role. The L1 fit on LR04 places near-zero amplitude at n = 24. The 100-kyr-band signal is carried by the adjacent multi-planet eigenmode beats above (n = 22, 25, 28). At T = 1.2 Myr the Rayleigh resolution at P = 110 kyr is ΔP ≈ 10 kyr, so H/3 sits within the same broad ~80–125 kyr band as the signal-carrying integers — but the L1 amplitude fit cleanly separates them and Earth’s own H/3 drops out.

More broadly: the model’s Earth-precession periods (Fibonacci divisors H/3, H/5, H/8, H/13, H/16) all sit at near-zero amplitude in LR04 (0.014 to 0.078, none above 0.085). This is structurally consistent: the Fibonacci H/n periods are Earth’s geometric precession in its own orbital frame; climate observes Earth in the heliocentric frame where planetary-induced orbital-plane motion shifts the signal to integers adjacent to the Fibonacci anchors (Berger’s k+g_j and k+s_j sub-peaks structure).

Eccentricity attribution has specific empirical headwinds

Direct-eccentricity attribution faces three discriminating empirical failure modes that the inclination-side framework does not:

Climate couples to ecliptic-frame periods, not ICRF

When another planet’s gravitational coupling registers in the paleoclimate spectrum, it appears at that planet’s ecliptic-frame perihelion period, not its ICRF-frame period. Mars’s ecliptic perihelion (8H/36 = 74.5 kyr) sits close to the LR04 n=35 peak at 76.6 kyr (permutation-null z ≈ +4.4, post-MPT); Saturn’s ICRF perihelion (8H/169 = 15.9 kyr) — a pure inertial-frame period with no ecliptic counterpart — shows no detectable LR04 signal (z ≈ −1). The apparent exceptions where an ICRF period does appear (Earth ICRF = H/3 inclination precession; Jupiter ICRF = 8H/65 obliquity beat) are exactly the cases where the ICRF period coincides with an independently climate-active ecliptic-frame feature.

This is physically expected: the ICRF has no dynamical coupling to the solar system’s internal secular motions, while Earth’s climate is set entirely by ecliptic-frame geometry. It is also why the 8H-lattice secular periods on Earth’s orbital plane (Jupiter 8H/39, Saturn −8H/65; see Fibonacci Laws — Law 6) are the climate-relevant attribution, rather than the Fibonacci anchors.

Ice-Albedo Decomposition

The climate formula’s frequency decomposition makes a quantitative prediction: how much of the temperature signal at each orbital band comes from ice-albedo feedback versus CO₂/greenhouse-gas forcing. Standard paleoclimate references give a single ice-albedo number averaged across the Last Glacial Maximum (Hansen et al. 2013  reports ~3.5 W/m² ≈ 55 % of LGM ΔF). The 8H-lattice fit produces a frequency-resolved decomposition — ice-share separated by orbital band — which has not been published in this form.

Method: Monte Carlo over the four calibration constants — ice fraction f_ice of LR04 δ¹⁸O (default 0.6, range [0.5, 0.7]), sea-level-per-mille ρ_SL (default 120 m/‰, [100, 140]), ice-albedo forcing per metre of sea level κ_F (default 0.04 W/m²/m, [0.030, 0.050]), and the other-GHG multiplier (CH₄ + N₂O + H₂O combined forcing relative to CO₂; default 0.5, [0.4, 0.6]). Per-band amplitude weights come from the L1 lattice fit. Results (post-MPT regime, 0–800 kyr BP):

The decomposition feeds the Charney equilibrium climate sensitivity calculation. Marginalized over the four calibration constants (N = 5,000 MC draws), the ΔT-weighted Charney sensitivity from the 8H-lattice frequency-resolved decomposition is 3.63 K, 90 % CI [3.01, 4.31] K — consistent with all major paleo and Bayesian-synthesis estimates (Hansen 2013 , Sherwood et al. 2020 , PALAEOSENS / Köhler et al. 2017 , IPCC AR6). Full synthesis on Climate Summary; literature comparison on Related Work.

The same decomposition computed pre-MPT reveals a sharp obliquity-band ice-share shift — the topic of §MPT obliquity-band ice-share shift below.

Pre-MPT vs Post-MPT — same forcing, different response

Canonical formula on full LR04 with three independent per-regime fits stitched together. Stitched R² = 0.93. The amplitude growth from left to right is the MPT story: orbital forcing is essentially stationary across 5.3 Myr, but ice-sheet hysteresis amplification grew dramatically after iNHG (~2.7 Ma) and again after MPT (~1 Ma).

Three structural transitions in the full LR04 record: iNHG (~2.7 Ma, Northern Hemisphere ice sheets establish), Early Pleistocene “41-kyr world” (2.7–1.2 Ma, obliquity dominates), MPT (~1.2–0.7 Ma, shift to “100-kyr world” as ice-sheet hysteresis crosses the Willeit threshold).

The orbital forcing — planetary eigenmodes, Earth’s precession, obliquity oscillation — is essentially constant over the past 5 Myr. What changed is the climate system’s sensitivity: long-term CO₂ decline, growing Northern Hemisphere ice sheets, and ice-sheet hysteresis crossing threshold (~1 Ma). The L1 lattice positions are stationary across the full 5 Myr; per-line amplitudes change across regimes, which is why the formula is fit per regime.

Windowed amplitude analysis across the MPT (pre-MPT 1,500–2,500 kyr vs post-MPT 0–1,000 kyr) gives concrete growth ratios:

The 41-kyr peak decreased — consistent with the standard “ice-sheet saturation silences the obliquity pacemaker” framing (Willeit 2019). The 100-kyr and 23-kyr bands both grew.

MPT obliquity-band ice-share shift

The MPT was not just an amplitude redistribution; it was a mechanism shift in how the climate system responded to obliquity forcing. The frequency-resolved ice-albedo decomposition (§Ice-Albedo Decomposition) computed separately for pre-MPT (≈ 1–2 Myr BP) and post-MPT (0–800 kyr BP) windows shows the ice-share at the obliquity band drops by +31 percentage points across the MPT transition:

This quantitatively confirms Willeit et al. 2019 ’s framing: pre-MPT was the obliquity-paced ice world — ice sheets responded almost entirely at obliquity periods, with CO₂ playing a smaller relative role. The Pleistocene MPT pushed ice-sheet hysteresis past the Willeit threshold; the obliquity ice-coupling weakened to today’s 64 %; and the 100-kyr regime took over.

The 100-kyr band stays flat across the MPT. This is the discriminating evidence that the obliquity shift is genuine signal rather than a uniform proxy bias — a systematic artifact would shift all bands equally. Boron-isotope δ¹¹B-derived CO₂ reconstructions (Chalk et al. 2017 , Dyez et al. 2018 , Martínez-Botí et al. 2015 , de la Vega et al. 2020 ) cross-validate the post-MPT decomposition: boron and ice-core CO₂ agree on Charney ECS to within ~10 %. The same boron data independently confirm the CenCO2PIP 100-kyr-smoothing artifact at deeper time (iNHG-MPT obliquity ECS drops from 6.4 K to 1.8 K when boron CO₂ replaces CenCO2PIP); since the smoothing affects all bands roughly equally, the differential ice-share signal (obliquity +31 pp vs 100-kyr flat) survives the artifact correction.

Forward Projection

Forward projection within the post-MPT regime. Past (left of ‘today’) overlays formula on LR04; future (right of ‘today’) is the formula extrapolated forward 250 kyr. Next predicted natural glaciation peak: ~60,500 AD (58,500 yr from now). Warmest interglacial: ~166,000 AD. Strongest glaciation in window: ~198,500 AD. Forward projection stays in scope only within the post-MPT regime.

The formula is fit on the post-MPT regime alone (0–1 Ma, R² = 0.87) and extrapolated forward 250 kyr. Ridge regression bounds the peak-to-peak forward range to ~4.7‰ normalised (plain OLS would extrapolate to ~24‰, unphysical).

Predicted glacial maxima:

Predicted interglacial peaks: 92,500 yr (~94,500 AD, C = −0.62); 120,000 yr (~122,000 AD, −0.54); 145,000 yr (~147,000 AD, −1.44); 164,000 yr (~166,000 AD, −2.31 — warmest in window).

The 9-kyr LGM offset in past validation (predicted ~29 kyr BP vs observed ~20 kyr BP) suggests surface-temperature response can lag the orbital clock by several kyr; the actual surface peak-warmth and glacial-onset dates may differ from the orbital values above.

Comparison with other published forecasts

All five identify that we sit in an unusually weak eccentricity minimum and predict a longer-than-typical interglacial. Without anthropogenic forcing, the consensus next-onset range is ~50–80 kyr. Mechanism differs: Berger-Loutre / Loutre-Berger / Tzedakis attribute the 100-kyr glacial cycle to direct eccentricity forcing; Ganopolski adds nonlinear CO₂-ice-sheet hysteresis; the Holistic formula attributes the 100-kyr-band centroid to the s₁−s₄ nodal eigenmode beat (§100-kyr cycle) with eccentricity-attribution failing three discriminating empirical tests.

Distinctive commitment: the Holistic formula commits to a fixed integer-divisor lattice (positions established from observed orbital motion, not fitted to climate). Climate peaks at non-lattice frequencies would refute the L1 layer; the closure test confirms zero such peaks in LR04.

Validation against past glacials

The LGM offset is expected: pure orbital forcing predicts the insolation drive, not the ice-sheet response; ice sheets carry memory and lag by several kyr. The MIS 6 match within 2 kyr is excellent.

Pacing shift and the late-Pliocene analogue

The forward projection departs from the regular ~100-kyr pacing of the past 700 kyr. Strong-to-strong interval is 138 kyr (58.5 → 196.5 kyr from now), with intermediate wiggles alternating between obliquity-band (43–48 kyr) and precession-band (22–25 kyr) intervals. Mean of all 5 local maxima: 34.5 kyr.

Why. The L1 layer is exactly 8H-periodic by construction — every L1 component is an integer divisor of 8H, so C_L1(t) = C_L1(t + 8H). This means the L1 orbital backbone for the next 250 kyr is byte-identical to the L1 backbone of 2.43–2.68 Ma BC — the orbital forcing is mathematically the same window. The late-Pliocene LR04 record is therefore the direct empirical anchor for what an orbital-coupled response to this exact L1 backbone looks like:

The two responses differ dramatically because climate sensitivity changed at the MPT (ice-sheet hysteresis, CO₂ feedbacks). The late-Pliocene LR04 record shows what an orbital-coupled response to this exact L1 backbone looks like — moderate-amplitude, 41-kyr-paced, without the extreme glacial swings of the post-MPT regime.

Climate response is a separate question. Whether Earth’s surface climate actually tracks this mixed orbital pacing depends on ice-sheet hysteresis. Three scenarios: (A) ice sheets shrink below the Willeit 2019 threshold and the response re-couples to the orbital clock (a “return to 41-kyr world” at the L1 backbone level); (B) post-MPT hysteresis persists and the climate keeps producing ~100-kyr cycles even with weaker orbital triggers (skipped or subdued glacials, matching the very-mild intermediate peaks in the projection); (C) anthropogenic CO₂ overrides orbital pacing entirely for the near-term (Ganopolski 2016 ). The 8H formula speaks to the orbital half only.

Limitations

Within-regime descriptive power is high; cross-regime predictive power is zero — formulas trained on pre-MPT data give catastrophically negative R² when applied to post-MPT, and vice versa. The L1 lattice positions are regime-independent, but per-line response amplitudes evolve with boundary conditions (CO₂, ice-sheet area). Forward projection therefore stays in scope only within the post-MPT regime (no scheduled boundary-condition shift in the next 250 kyr).

Anthropogenic CO₂ may delay the next natural glaciation by 50+ kyr in moderate-emission scenarios (Ganopolski et al. 2016 ) or up to ~100 kyr in high-emission scenarios. Not modelled here.

Methodology highlights

Reproducibility scripts: 3d repository .

Data sources

• LR04 stack — Lisiecki & Raymo 2005 , Paleoceanography 20, PA1003. 57 benthic δ¹⁸O records, 0–5,320 kyr BP. Primary spectral analysis; per-regime L1+L2+L3 fits.
• CENOGRID composite — Westerhold et al. 2020 , Science 369, 1383. Cenozoic δ¹⁸O + δ¹³C global composite, 0–67 Myr BP. Deep-time L3 step generalization.
• EPICA Dome C CO₂ — Bereiter et al. 2015 , Geophys. Res. Lett. 42, 542. Antarctic ice-core atmospheric CO₂, 0–805 kyr BP. Cross-proxy validation.
• CenCO2PIP CO₂ — Consortium 2023 , Science 382, eadi5177. Bayesian multi-proxy CO₂ synthesis, 0–66 Ma. Deep-time CO₂ validation.
• Cheng 2016 Asian speleothem composite — Cheng et al. 2016 , Nature 534, 640. U-Th-dated Asian Monsoon δ¹⁸O, 0–640 kyr BP. Non-tuned chronology-bias control.

See also

Climate cluster siblings: Climate Summary · L1 Attribution · Insolation Null Test · Related Work

Elsewhere: Eigenfrequencies · Obliquity & Inclination · Eccentricity · Why Earth Is Special · Supporting Evidence §1 (100-kyr problem + ice-core dating methodology) · 8H period table