Skip to Content
📄 Fibonacci Laws — Read the paper
The ModelClimate Formula

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.

Orbital forcing is not climate. L1 captures orbital forcing; L2 adds the silicate-weathering 405-kyr carbon-cycle thermostat; L3 adds Heaviside step terms at known Cenozoic transitions. The formula does not model ice-sheet hysteresis dynamics, CO₂ amplification feedbacks beyond the L2 thermostat, internal variability (Heinrich/D-O), regional asymmetries, or anthropogenic CO₂. Orbital cycles set the timing of glacial/interglacial transitions; the magnitude of the observed signal is dominated by Earth-system feedbacks not modelled here.

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 .

Climate Formula on LR04 δ¹⁸O over the past 700 kyr (post-MPT regime). Model curve (red) tracks the observed LR04 data (black) through seven glacial-interglacial cycles. Per-regime R² = 0.87.

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. 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

Computational reference card: For a copy-paste-ready formula card with definitions, validity windows, denormalization step, JSON coefficient locations, and a worked numerical walk-through, see Formulas → §8 Canonical Climate Formula. This page is the scientific narrative; the reference page is the math card.

C(t)  =  c0  +  nNL1[ancos ⁣(2πnt8H)+bnsin ⁣(2πnt8H)]L1: orbital forcing  +  pPL2[αpcos ⁣(2πtp)+βpsin ⁣(2πtp)]L2: carbon thermostat  +  iTL3γiH(tti)L3: boundary-condition stepsC(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.

RegimeRecord + windowR² L1ΔR² L2ΔR² L3Total R²
post-mptLR04 0–1 Ma0.870+0.003+0.0000.874
inhg-mptLR04 1.0–2.7 Ma0.722+0.007+0.0000.729
pre-inhgLR04 2.7–5.3 Ma0.381+0.049+0.0000.430
lr04-fullLR04 0–5.3 Ma0.239+0.009+0.0080.255
epica-co2EPICA 0–800 kyr0.834+0.012+0.0000.845
cenco2pipCenCO2PIP 0–66 Ma0.161+0.000+0.6020.763
cenogrid (δ¹⁸O)CENOGRID 0–67 Ma0.001+0.004+0.6130.618
cenogrid (δ¹³C)CENOGRID 0–67 Ma0.001+0.008+0.3420.351

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.

Climate Formula on EPICA Dome C atmospheric CO₂ over 0–800 kyr BP. Same canonical 3-layer formula and the same 32 integer-divisor L1 lattice, fitted to a completely different proxy (atmospheric CO₂ from Antarctic ice cores rather than benthic δ¹⁸O). R² = 0.85 (L1 alone = 0.83).

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:

nperiod (kyr)physical interpretation
9298.1g₂−g₇ eccentricity beat / Mercury Axial = AscNode = 8H/9
22121.9s₂−s₄ nodal / g₄−g₂ — strongest L2-signature lattice integer
25107.3s₁−s₄ nodal eigenmode beat (100-kyr centroid)
2895.8g₄−g₅ eccentricity beat (Berger 95k)
3968.8Jupiter ecliptic perihelion = 8H/39 (≈ Laskar s₃ eigenmode)
6541.27k+s₃ obliquity beat (Berger 41k) = Saturn ecliptic = Jupiter ICRF — Law 6 resonance lock; dominant peak
6839.5k+s₄ obliquity sub-peak = Mars ICRF perihelion = 8H/68
11323.74k+g₅ climatic precession sub-peak (Berger 23.7k)
12022.35k+g₂ climatic precession sub-peak = H/15 (Berger 22.4k)

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.

Try it interactively. The 3D simulation  ships with a Climate Formula Explorer modal (Tools menu) that plots the formula reconstruction directly on top of four proxy records across eight tabs: CenCO2PIP · 66M, CENOGRID · 67M, LR04 · 5.3M, LR04 · 1.2M, EPICA · 800k, LR04 · 700k, LR04 · 200k, and LR04 forward (–250 → +250 kyr). Three layer toggles (Total, L1 alone, L2 alone) and a per-regime R² breakdown panel.


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

PlanetDirect L1 matchesExclusiveNotes
Jupiter4n=39 (Peri_ecl)n=16 Obliq, n=21 Axial, n=39 Peri_ecl, n=65 ICRF
Mars3n=68 (ICRF)n=16 Axial, n=21 Obliq, n=68 ICRF
Mercury2 (same n)n=9 (Axial = AscNode)Cassini-state lock
Venus2 (same n)n=110 (ICRF = Obliq)Two cycles co-located
Saturn1n=65 shared with Jupiter ICRF — the Law 6 lock
Uranus1n=16 shared with Mars + Jupiter Obliq
Earth0All six Earth cycles sit off the L1 lattice (Fibonacci H/n anchors are deliberately off-lattice)
Neptune0All Neptune cycles off L1

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

PlanetPublished periodReferenceModel H/nDeviation
Mercury895,000 yrBills & Comstock 20058H/3 = 894,179 yr+0.09%
Earth41,000 yrLaskar 2004; Berger 1978H/8 = 41,915 yr+2.2%
Mars124,800 yrWard 1973; Laskar 20048H/21 = 127,740 yr+2.4%

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:

TestEccentricity attributionInclination-side framework
405-kyr term presence (Berger’s strongest predicted eccentricity beat)Predicts dominant; observed essentially absent in post-MPT LR04 (amplitude ratio 0.12)Not predicted as an orbital signal here — modelled as Layer-2 carbon-cycle internal resonance (§Framework validation)
Eccentricity-beat phase coupling (95k+125k bispectrum)Predicts strong; not detected (bicoherence 0.507 < null-95 0.555 — replicates Muller & MacDonald 1997)No contradiction
100-kyr centroid identificationPredicts 95–99 kyr eccentricity beats; observed centroid at 107 kyr is a nodal beats₁−s₄ nodal eigenmode at n=25 = 107 kyr — matches
Berger climatic-precession 6-peak spectrumNot predicted (these are k+g_j sub-peaks, not eccentricity-beat periods)All 6 peaks match H/n divisors to <0.4%

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.


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

Climate Formula on full LR04 (0–5,320 kyr BP). Three-regime stitched fit: pre-iNHG, iNHG-MPT, post-MPT. Transition markers at iNHG (~2.7 Ma) and MPT (~1 Ma) visible as dashed vertical lines. Amplitude grows from past (left) to present (right), reflecting climate sensitivity changes — not orbital forcing changes.

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:

BandPre→post ratioDirection
41-kyr obliquity0.72×shrank
100-kyr band (n=25, 28)1.64×grew
23.7-kyr climatic precession (n=113)2.19×grew most
405-kyr empirical line (L2 carbon-cycle resonance)0.34×shrank

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.


Forward Projection

Climate Formula forward projection from −250 to +250 kyr. Past = formula + LR04 data; future = formula extrapolation only. Markers: next glacial onset at ~60,500 AD, warmest interglacial at ~166,000 AD, strongest glacial at ~198,500 AD.

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:

Years from nowAD dateC(t) normalisedStrength
58,500~60,500 AD+2.27next natural glaciation onset
106,000~108,000 AD+0.24mild
130,500~132,500 AD+0.04very mild
153,000~155,000 AD−0.99(local max, interglacial-range)
196,500~198,500 AD+2.48strongest in window

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

FrameworkNext glacial onset (kyr)Method
Berger & Loutre 2002 ~50LLN-2D astronomical-insolation model; current eccentricity minimum (MIS 11 analogue)
Loutre & Berger 2003 ~50–100Same + low-CO₂ scenarios
Tzedakis et al. 2012 ~50 (analog-bound)MIS 19c astronomical analogue + pre-industrial CO₂ inception threshold
Ganopolski et al. 2016 ~50 / ~100 / ≥500CLIMBER-2 climate-ice-sheet model (no / moderate / high anthropogenic CO₂)
Holistic Climate Formula~5832-integer 8H lattice + 3-line L2 carbon thermostat + 6-step L3 transitions, ridge fit on LR04 + CENOGRID + EPICA + CenCO2PIP

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

Formula predictionActual paleoclimate eventDiscrepancy
Glacial maximum at 29 kyr BPLast Glacial Maximum ≈ 20 kyr BP+9 kyr early
Glacial maximum at 138 kyr BPMIS 6 ≈ 140 kyr BP−2 kyr

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:

MetricLate Pliocene (2.43–2.68 Ma BC)Post-MPT (0–250 kyr BP)
LR04 correlation with formula C(t) [r]0.7600.945
LR04 raw detrended std (‰)0.2570.464 (1.8× Pliocene)
Glacial-peak intervals (kyr)mean 39.5mean 29.3
Dominant LR04 spectral period41.3 kyr (obliquity)125.5 kyr (100-kyr band)

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

TopicDetail
Spectral methodsMTM (Thomson multitaper, NW=3, K=5); Lomb-Scargle for Cheng2016’s irregular sampling; Hinich bispectrum for eccentricity-beat coupling; multi-component OLS for amplitude fit; single-component OLS for integer-divisor scan
Rayleigh resolutionΔP ≈ P²/T. At T = 5,320 kyr (full LR04) and P = 41 kyr: ΔP ≈ 0.32 kyr (Berger sub-peaks resolvable). At T = 1.2 Myr and P = 110 kyr: ΔP ≈ 10 kyr (95k/100k/112k spectrally collinear).
Pre-registrationEvery empirical test had data source, method, parameters, and verdict rules locked in writing before any data analysis
Bootstrap1,000 resampling iterations for amplitude confidence intervals

Reproducibility scripts: 3d repository .


Data sources

  • LR04 stackLisiecki & Raymo 2005 , Paleoceanography 20, PA1003. 57 benthic δ¹⁸O records, 0–5,320 kyr BP. Primary spectral analysis; per-regime L1+L2+L3 fits.
  • CENOGRID compositeWesterhold 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 compositeCheng 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

Last updated on: