Interglacials of the Quaternary defined by northern hemispheric land ice distribution outside of Greenland

The study addresses how interglacials of the Quaternary should be defined and paced, particularly across the Mid-Pleistocene Transition (MPT) when dominant glacial cycle periodicity shifted from ~41 kyr to ~100 kyr. Traditional orbital theory and recent rules (e.g., caloric summer half-year insolation at 65°N) link deglaciations to obliquity maxima, but these approaches largely infer ice-volume changes from benthic δ18O and may conflate temperature and ice-volume signals. The authors propose redefining interglacials by the absence of substantial northern hemispheric land ice outside Greenland, using land-ice distributions inferred from model-based deconvolution of the LR04 benthic δ18O stack. The goal is to test whether interglacial onsets are regularly paced by obliquity throughout the last 2.6 Myr and to reassess the MPT dynamics with a land-ice-based criterion that better reflects ice-sheet physics and albedo feedbacks.

Prior work established orbital forcing (Milankovitch) as pacing glacial cycles, with early Quaternary 41-kyr obliquity dominance and later quasi-100-kyr cycles. Explanations for the MPT include internal feedbacks: ice-sheet dynamics and regolith removal, and gradual declines in glacial CO2. Tzedakis et al. (2017, T17) proposed a simple insolation-based rule (caloric summer half-year at 65°N) to predict which cycles yield interglacials, defining interglacials via thresholds in a detrended benthic δ18O stack, aligning with a community synthesis for the last 800 kyr. Other glaciological studies deconvolved benthic δ18O into sea level (ice volume) and temperature using 3-D ice-sheet models (e.g., ANICE), implicating North American ice-sheet mechanics (e.g., merging domes, basal sliding at pressure-melting conditions) in post-MPT longer cycles. Simplified insolation metrics often emphasize a single latitude (65°N), potentially underrepresenting latitudinal ice-distribution effects and the full-year energy budget. The literature also discusses CO2’s role as feedback vs. trigger, with mixed proxy evidence and modeling suggesting increased glacial-interglacial CO2 amplitude across the MPT due to ocean circulation and iron fertilization feedbacks.

• Data basis: Use the ANICE 3-D ice-sheet model-based deconvolution of the LR04 benthic δ18O stack to separate deep-ocean temperature and global mean sea-level (ice volume) and to simulate spatial-temporal ice distribution for four ice-sheet regions (North America, Eurasia, Greenland, Antarctica) at 40×40 km (Greenland 20×20 km). Results cover the last 5 Myr; analysis restricted to the Quaternary. Model output temporal resolution: 2 kyr, interpolated to 1 kyr for onset detection.
• Interglacial definition: An interglacial is defined by the absence of substantial northern hemisphere land ice outside Greenland, operationalized via thresholds on the integrated ice volume of North America + Eurasia. Two thresholds are used: Threshold 1 (lower) must be crossed downward to flag a new interglacial; Threshold 2 (upper) must be crossed upward to separate two interglacials by an intervening glacial. If Threshold 1 is not crossed within an obliquity cycle (minima-to-minima), that cycle is a skipped termination; if Threshold 2 is not crossed, it is a continued interglacial.
• Threshold calibration: Chosen to (i) reproduce the interglacial detections during the last 800 kyr from a community review, and (ii) minimize the number of obliquity cycles without new interglacials earlier. Standard, time-varying thresholds maintain a constant difference between the two thresholds. Numerical values (in 10^15 m^3): after 1650 kyr BP, n1=7 and n2=10; before 1750 kyr BP, n1=3 and n2=6; linearly interpolated between 1750–1650 kyr BP. Alternative sensitivity tests use time-invariant high or low threshold pairs.
• Cycle assignment: Onsets are assigned to obliquity cycles defined by consecutive obliquity minima, without requiring phase alignment to insolation maxima, thus reducing sensitivity to age uncertainties.
• Additional analyses: Latitudinal distribution metrics of NA+EA ice volume (e.g., latitude of 50% and 10–90% extent), fraction of global ice-volume change attributable to NA+EA, and fraction of benthic δ18O attributable to ice-volume (δ18Osw). Frequency analyses performed with Analyseries. CO2 radiative forcing calculated as ΔR = 5.35 W m^-2 ln(CO2/278 ppm).

• Interglacial realization fractions (new interglacial per obliquity cycle): pre-MPT 67% (16/24, MIS 55–101), MPT 88% (14/16, MIS 27–53), post-MPT 52% (12/23, MIS 1–25). This indicates irregular interglacial occurrence throughout the Quaternary, not strictly obliquity-paced.
• Mean interglacial return times: pre-MPT 60 ± 22 kyr, MPT 47 ± 19 kyr, post-MPT 79 ± 24 kyr, illustrating misleading average periodicities due to skipped and continued cycles.
• Early Quaternary glacial-interglacial amplitudes in NA+EA land ice are small relative to post-MPT, despite prominent benthic δ18O variability; prior to the MPT, only ~20% of the 41-kyr running-mean δ18O signal is due to sea level/ice volume, rising to ~50% after the MPT during glacials.
• The fraction of global ice-volume change due to NA+EA declines to ~50% during long early interglacials but remains ~85% during the last 1 Myr, explaining discrepancies with δ18O-based interglacial identification.
• Multiple early Quaternary obliquity cycles lack new interglacials under the land-ice criterion (e.g., skipped terminations around MIS 61, 65, 71, 77; continued interglacials around MIS 55, 79, 93, 101), differing from T17 in 4–5 robust cases.
• Latitudinal ice distribution shifted southward over time: at glacial maxima, the NA+EA ice 50% line moved from ~70°N (pre-1.8 Ma) to 62°N (LGM). The 10% southern extent moved from 55–63°N early to ~50°N after the MPT.
• The climate system shifted from an interglacial-dominated early Quaternary to a glacial-dominated late Quaternary when using constant thresholds.
• CO2 likely acts predominantly as a feedback across the MPT: modeling reproduces the 41→100 kyr transition without externally forced CO2 declines, though CO2 amplitude increases likely contributed to longer cycles; glacial/interglacial transitions can occur with fixed CO2 in models, but without lengthened cycles.
• Results converge with independent skipped-termination analyses for much of the last ~1.6 Myr, but not for the early Quaternary, highlighting definition sensitivity.

By defining interglacials via minimal NA+EA land ice rather than benthic δ18O minima, the study shows that terminations are irregularly realized across obliquity cycles throughout the last 2.6 Myr. This challenges simple obliquity-paced rules and indicates that δ18O minima before the MPT often reflect deep-ocean temperature changes and ice outside NA+EA rather than substantial deglaciation of northern hemisphere continents. The findings underscore the dominant role of ice-sheet dynamics (including basal conditions and possible regolith removal) in setting cycle length after the MPT and suggest CO2’s role is primarily as a feedback amplifying cycle length rather than the initial trigger. The southward shift and expansion of NA+EA ice through the Quaternary altered albedo feedbacks and sensitivity to insolation across latitudes, implying that latitude-specific insolation metrics are of secondary importance compared to actual ice-distribution changes. Differences with δ18O-based interglacial catalogs (e.g., T17) arise mainly in the early Quaternary where ice-volume signals are small, reinforcing the notion that binary glacial/interglacial classifications depend strongly on chosen variables and thresholds. The study advocates focusing on gradual climate indices and the physical mechanisms (ice dynamics, carbon cycle feedbacks) rather than rigid binary definitions, particularly for early Quaternary intervals.

The paper introduces a land-ice-based definition of interglacials using model-derived NA+EA ice volumes, revealing that interglacial onsets are irregular across the Quaternary, with realization fractions of 67% (pre-MPT), 88% (MPT), and 52% (post-MPT). It demonstrates that early Quaternary δ18O variability overstates NA+EA ice-volume changes, explaining discrepancies with δ18O-based interglacial catalogs. The work highlights a transition from interglacial-dominated to glacial-dominated climate states and a southward shift of ice sheets across the Quaternary, supporting a primary role for ice-sheet physics and a feedback role for CO2 in the MPT. Future research should (i) develop independent constraints on NA+EA ice volumes and their latitudinal distribution, (ii) expand global δ18Osw/sea-level stacks with broader spatial coverage, (iii) improve ice-sheet model physics and parameterizations, (iv) refine CO2 reconstructions and age models across the MPT, and (v) emphasize continuous climate indices over binary classifications to better capture the mechanisms driving Quaternary climate variability.

• Dependence on the LR04 benthic δ18O stack: any biases (stacking method, spatial heterogeneity, undersampling, age model) propagate into the deconvolution and inferred ice volumes.
• Model-based uncertainties in ice-sheet deconvolution and parameterizations; estimated average uncertainty for NA+EA ice volume ~0.5×10^15 m^3 (~2%).
• Coarse temporal resolution (2 kyr, interpolated to 1 kyr) limits assessment of millennial-scale variability and precise termination phasing.
• Threshold selection introduces arbitrariness; early Quaternary small amplitude changes make classification sensitive to threshold choice.
• Limited independent geological constraints on NA+EA minimum ice extents; notable mismatch with evidence for an early Pleistocene Laurentide advance to ~38–39°N around 2.5 Ma.
• Alternative δ18O stacks (e.g., probabilistic stacks) differ around 1.8–1.9 Ma (MIS 65–71), potentially affecting classifications in those cycles.
• Sea-level and δ18Osw reconstructions used for indirect comparison may be regionally biased or insufficiently sampled to represent global means.