Global Evidence of Obliquity Damping in Climate Proxies and Sea-Level Record during the Last 1.2 Ma: A Missing Link for the Mid-Pleistocene Transition

Preprint

Article

Global Evidence of Obliquity Damping in Climate Proxies and Sea-Level Record during the Last 1.2 Ma: A Missing Link for the Mid-Pleistocene Transition

Altmetrics

Downloads

142

Views

Comments

A peer-reviewed article of this preprint also exists.

Paolo Viaggi^*

This version is not peer-reviewed

Submitted:

06 September 2023

Posted:

08 September 2023

You are already at the latest version

Alerts

Abstract

A recent research has identified an inverse amplitude link between obliquity damping and short eccentricity amplification during Mid-Late Pleistocene based on LR04 δ18O and equatorial Pacific Site 846 sea surface temperature records, which is associated with the Earth's long-term cooling. In the present study, new evidence of this anticorrelation is presented from Antarctic D-CO2-CH4 records, global benthic-planktic δ18O, and regional (Atlantic, Pacific, Mediterranean, and Indian) climate-related proxies. Based on a critical review of theoretical constrains (Earth's oblateness changes and ice-volume phase lag in the obliquity band <5.0 kyr), this widespread and symmetric (bipolar) obliquity response damping has been interpreted as an effect of the obliquity-oblateness feedback, which could be the latent physical mechanism at the origin of the Mid-Pleistocene Transition (MPT). Indeed, results and considerations of the present work suggest that fast and positive/negative net variation of the Earth's oblateness in the obliquity band was controlled by dominant glacio-eustatic water mass component and, assuming a rapid response of the ice-volume to surface temperature changes, the mean obliquity lag response is estimated to be <5.0 kyr over the past 800 kyr. These elements may explain the interglacial/glacial damping observed in the obliquity response. The consolidation of the Earth's long-term icy-state in the subtrend IV, culminated with the post-MPT obliquity damping, might have contributed to the strengthening of the short eccentricity response by mitigating the obliquity ‘ice-killing’ during obliquity maxima (interglacials), favouring the obliquity-cycle skipping and a feedback amplified ice-growth in the short eccentricity band (obliquity damping hypothesis). This could suggest a different impact of the climate friction than what is generally believed, which is presumably the latent physical mechanism that triggers transient ‘competitive’ interaction between obliquity and short eccentricity started early during the Piacenzian.

Keywords:

Subject: Environmental and Earth Sciences - Geophysics and Geology

1. Introduction

Understanding the Mid-Pleistocene Transition has remained one of the most fascinating and unresolved question in Paleoclimatology over the course of several decades of research. MPT is marked by a gradual shift in periodicity from 41 kyr to ~100 kyr and a progressive increase in the amplitude and ‘saw-tooth’ shape of climate oscillations. This scientifically challenging phenomenon represents a key element to constrain the mechanisms that governed the climate system in the past, and it remains a fundamental learning lesson also for present-day climate and the associated changes. The MPT began at approximately 1.2 Ma (Head and Gibbard, 2005; Clark et al., 2006; Head et al., 2008) and marks the beginning of the classical Pleistocene ‘Ice Ages’ when major glaciations developed over North America, northern Europe, and the Alps. The final stage of the MPT at approximately 0.7 Ma (Clark et al., 2006) marks the start of the dominant ~100-kyr climate response, a paleo-climatological paradox considering the extremely low energy of the eccentricity insolation changes (Imbrie et al., 1993; Berger et al., 1999; Clark et al., 2006; Berger and Loutre, 2010; Lisiecki, 2010). Indeed, new insights into the nonlinear dynamics of the climate indicate short eccentricity response amplification over nominal solution of +400% (global δ¹⁸O) and +180% (equatorial Pacific sea surface temperature, SST) during Mid-Late Pleistocene, which are linked with the long-term mean icy-state (δ¹⁸O and SST exponential trends) (Viaggi, 2018a). This shift in periodicity, amplitude, and asymmetry of climate responses towards the Late Pleistocene evinces the fundamental changes in the dynamics of climate system mentioned in the literature, however, the cause of these features still remains poorly understood (Clark et al., 2006; Bintanja and van de Wal, 2008; Ellis and Palmer, 2016; Chalk et al., 2017; Köhler and van de Wal, 2020). In the following lines, we critically examine some of the many models and hypotheses related to the MPT, referring to the cited literature for discussion and further details.

One way to address MPT questions is to replicate the paleoclimatic records using a variety of glacial cycle models of different complexity, referred to as conceptual models (Imbrie et al., 2011; Crucifix, 2012; Ditlevsen and Ashwin, 2018; Quinn et al., 2018; Mukhin et al., 2019; Nyman and Ditlevsen, 2019; Berends et al., 2021). On the other hand, the numerous mechanisms through which the 100-kyr cycles can be generated make it difficult to determine the best model to describe the source of the 100-kyr variability (Huybers, 2011; Imbrie et al., 2011). Moreover, since models cannot accurately reproduce the observed record, they lack information on certain physical processes (Berends et al., 2021). In some models, the MPT is represented either by an increase in the nonlinearity of the climate system response to obliquity forcing (Huybers, 2007; Liu et al., 2008; Nyman and Ditlevsen, 2019) or by an enhanced sensitivity of climate response to eccentricity owing to modulation of the orbital precession (Imbrie et al., 1993; Raymo, 1997; Shackleton, 2000; Lisiecki, 2010). Based on obliquity-pacing results, Huybers (2007) proposed that the 100-kyr variability could arise from the skipping of one or two obliquity beats, corresponding to 80-kyr or 120-kyr glacial cycles, which provide an average periodicity of ~100 kyr, resulting from the long-term growth of ice-sheet (obliquity-cycle skipping). However, this model produces approximately constant power across the MPT for the 41-kyr cycle, in disagreement with observational data (Lisiecki and Raymo, 2007; Viaggi, 2018a). The idea of skipped obliquity cycles appears also in a phase-space model that simulates changes in Pleistocene ice volume based on Earth’s orbital parameters (Imbrie et al., 2011). Their model demonstrates 40 kyr epochs with termination in each obliquity cycle and 100-kyr epochs where the model skips certain obliquity cycles when it fails to reach the threshold for terminations. Models of obliquity skipping are intriguing mainly because they introduce the idea of weakening obliquity rather than explaining the 100 kyr cycle as multiples of the 40 kyr cycle, and the physical mechanism underpinning the obliquity-cycle skipping remains unclear. Perhaps, the most interesting model is that proposed by Abe-Ouchi et al. (2013) in which the hysteresis loop of the North American ice sheet is such that after inception of the ice sheet, its mass balance remains mostly positive through several precession cycles, whose amplitudes decrease towards an eccentricity minimum. Their model simulates the sawtooth characteristic of glacial cycles, the timing of terminations and the amplitude of the Northern Hemisphere ice-volume at the Last Glacial Maximum with a dominant 100-kyr spectral peak. However, obliquity is not resulted the driver of the 100-kyr cycle, although it helps to amplify the ice-volume changes from glacial to interglacial states.

Most of the hypotheses put forward for the origin of the MPT invoke certain internal changes of the climate system in response to long-term cooling, possibly induced by a decrease in atmospheric CO₂ (Clark et al., 2006; Zachos et al., 2001, 2008; Kender et al., 2018; Shoenfelt et al., 2018; Farmer et al., 2019; Hasenfratz et al., 2019) along with the concurrent Northward geodynamic migration of Greenland and North America, favouring a net accumulation of continental Arctic ice (Steinberger et al., 2015; Daradich et al., 2017). Viaggi (2018a) reinforced the notion of a geo-tectonic control on the long-term composition of the atmospheric greenhouse gases (GHG) (seafloor spreading rate, explosive volcanic activity, orography and erosion, paleogeographic configuration, oceanic paleocirculation, and ocean fertilisation), and highlights four phases (subtrends I, II, III, IV) of progressive lowering of the average temperatures and ice-volume growth, triggered by long-term changes in the atmospheric composition. However, these interesting hypotheses might explain the origin of the Plio-Pleistocene long-term cooling, but it remains unclear what kind of ‘response’ to this cooling has to deal with the MPT origin and its observational features. Indeed, in the context of increasing power of obliquity nominal forcing and decrease in eccentricity nominal power (Laskar et al., 2004; Laskar et al., 2011) during MPT and post-MPT period, it remains unclear why obliquity cycles which dominated the pre-MPT time with an incisive ‘ice-killing’ action (Berger et al., 1999; Lisiecki and Raymo, 2007; Liu et al., 2008; Viaggi, 2018a), gradually became less effective in preventing the prolonged survival of the ice-sheet through possible obliquity-cycle skipping (Huybers, 2007; Liu et al., 2008) and obliquity response damping, associated to strong short eccentricity response (Viaggi, 2018a, 2021a).

In this study, we present new observational evidence of damping of the obliquity response and a physical mechanism (obliquity-oblateness feedback) that may explain its origin, which could be related with the MPT origin.

1.1. Obliquity-oblateness feedback

Viaggi (2018a) quantitatively investigated the strength relationships of orbital climate responses to nominal forcings during the Plio-Pleistocene, and showed a sharp decline in the obliquity response sensitivity of both global benthic δ¹⁸O stack (Lisiecki and Raymo, 2005) and equatorial Pacific Site 846 SST (Herbert et al., 2010) during the Mid-Late Pleistocene, coupled with a strong amplification of short eccentricity response, linked to the Earth’s long-term cooling. This impressive coupling between obliquity-damped and short eccentricity-amplified responses could be an understudied feature in the MPT debate. On the other hand, Viaggi (2018a) hypothesised an attenuation mechanism on the obliquity forcing by obliquity-oblateness feedback (Rubincam, 1993, 1995; Bills, 1994, 1998; Williams et al., 1998; Levrard and Laskar, 2003; Laskar et al., 2004) and a reduction of feedback amplification processes during obliquity maxima (interglacials). This might have led to weakening of the obliquity ‘ice-killing’, thereby, favouring a ~100-kyr long-life feedback-induced ice growth in the short eccentricity band (here referred as obliquity damping hypothesis, ODH). However, this ‘early stage’ idea requires more comprehensive development, which is the purpose of the present work.

Climate friction by obliquity-oblateness feedback is dissipative feedback between obliquity variations and climate, which may cause a secular drift of the spin axis (Laskar et al., 2004). Glacial and interglacial conditions drive the redistribution of the ice/water mass and the isostatic adjustment to surface loading, affecting the dynamical ellipticity of the Earth (oblateness). Because both the Earth’s ice load history and viscoelastic structure are not strongly constrained, it is difficult to produce accurate predictions of the coupled response of the entire system, which depends by both oblateness changes and phase lag estimates of the ice-volume in the obliquity band (Williams et al., 1998; Levrard and Laskar, 2003; Laskar et al., 2004; Skinner and Shackleton, 2005; Lisiecki and Raymo, 2009). Therefore, simplified assumptions have to be made to evaluate the magnitude and the direction (positive or negative) of the obliquity secular change. Specifically, based on the analysis of benthic δ¹⁸O records, Levrard and Laskar (2003) estimated a positive mean secular obliquity change of 0.01° Ma^-1 during the glaciation periods of the last 3 Ma. However, ODH requires a negative secular change of obliquity during the last 1.2 Ma to explains the obliquity reduction on its maxima (interglacials), i.e., during the most critical phase to promote long-lasting feedback-induced ice growth in the short eccentricity band. The obliquity-oblateness feedback is still an understudied phenomenon that needs further investigation (Viaggi, 2018a), while the anomalous decline in obliquity response power during the Mid-Late Pleistocene (Lisiecki and Raymo, 2007; Liu et al., 2008; Lisiecki, 2014; Viaggi, 2018a) and the detailed analysis of the Plio-Pleistocene long-term cooling trend (Viaggi, 2018a) are new insights that need to be placed in the context of the MPT origin. The obliquity-oblateness feedback could be the latent physical mechanism connecting obliquity damping with the short eccentricity amplification.

The aim of the present work is to expand early stage idea of ODH by the following topics:

1): searching worldwide new observational evidences for the link between obliquity damping and short eccentricity amplification from global and regional (Antarctic, Pacific, Atlantic, Mediterranean, Indian) climate-related proxies (Sect. 3.1);
2): discuting the role of the long-term cooling trend in the MPT debate and the relationships between orbital forcings and proxies (Sect. 3.2 and 3.3);
3): by critically review the requisite theoretical constrains of ODH to establish that the obliquity-oblateness feedback could be the driving mechanism of the interglacial/glacial damping observed in Mid-Late Pleistocene obliquity responses (Sect. 3.4);
4): refreshing by new cross-spectral data the role of the short eccentricity forcing (Sect. 3.5).

The implications of ODH are that the onset of the ~100 kyr cycle associated with the MPT would be the climate system’s reaction to the weakening of the obliquity and the related strengthening, mediated by feedback mechanisms, of the phase-locked short eccentricity response.

1.2. Key role of Obliquity Forcing on the Earth’s Climate System

According to the Milankovitch theory, obliquity is a key component in high latitude insolation. High obliquity significantly amplifies the seasonality, especially in the polar areas, thereby creating cold winters and hot summers with maximal melting. This prevents high-latitude ice accumulation (Levrard and Laskar, 2003). Huybers (2006) argued that the integrated summer (June-July-August, JJA) insolation between 30° and 70 °N, which is dominated by the 41-kyr obliquity signal, is more representative of summer melt compared to the 65 °N JJA insolation, which is dominated by precession. Raymo and Nisancioglu (2003) proposed that the strong obliquity signal imprinted on the Ice Age record is generated by the exertion of controlled meridional temperature gradients on the poleward transport of moisture. As obliquity decreases, cooling at high latitudes occurs and the gradient of insolation heating between high- and low-latitude increases, both of which promote ice sheet growth. Similar results were obtained by Mukhin et al. (2019). They construct a Bayesian data- driven model from LR04 δ¹⁸O record that could account for the main factors which may potentially impact Pleistocene climate. The only insolation forcing that matters in the Pleistocene climate is that the meridional gradient of insolation is dominantly affected by obliquity oscillations. The insolation gradient was regarded as the driving force for the 41-kyr climate cycles during the pre-MPT epoch (from 3 Ma to 0.8 Ma), through the atmospheric meridional heat and moisture fluxes that were modulated (Mukhin et al., 2019). Moreover, an assessment of eleven radiometrically dated terminations suggests that obliquity exerted a persistent influence on ice age terminations since the MPT (Bajo et al., 2020). The key role of obliquity on the Earth’s climate system is also attested by its largest quantitative impact on the global LR04 δ¹⁸O record during the Plio-Pleistocene, whose variance is estimated to be 9.9% with only 2.0% precession (Viaggi, 2018a). Similar results were obtained for the SST record of equatorial Site 846, where obliquity carries 5.5% variance and precession of 1.6%. The relevant action of obliquity in preventing long life of ice-sheet during Plio-Pleistocene pre-MPT period was referred as obliquity’s ‘ice-killing’ by Viaggi (2018a).

Since obliquity plays a key role on the climate variability, any process or mechanism that can weaken obliquity forcing might have promoted the ‘skipping’ of obliquity-cycle, especially in terms of climate response to interglacial stages during the icy-state background. Thus, a hypothetical weakened obliquity during interglacials might have reduced seasonality and refreshed the polar summers, favouring high-latitude snow precipitation (poleward transport of moisture) and ice preservation (high latitudes cooling). This implies that the obliquity damping may have mitigated the obliquity’s ‘ice-killing’, favouring obliquity cycle skipping. Indeed, the skipping of obliquity cycle is more frequent during the Late Pleistocene, when most deglacial/glacial events are separated by low-amplitude obliquity peaks (Huybers, 2007; Viaggi, 2018a; this study).

2. Materials and Methods

2.1. Materials

This work contributes to the research by considering a variety of global and regional (Antarctic, Atlantic, Pacific, Mediterranean, and Indian) climate-related records from literature. The information reported on each record is more detailed for those records that were statistically processed in the present work:

Global proxies

The δ¹⁸O values of the Plio-Pleistocene LR04 benthic stack was obtained from 57 oceanic series, which are globally distributed over a wide latitudinal range (Atlantic, Pacific, and Indian Ocean) (Lisiecki and Raymo, 2005). The 1-kyr resampled orbital components isolated by singular spectrum analysis (SSA) (Viaggi, 2018a) were applied in the present work. The benthic δ¹⁸O measures changes in global ice-volume and deep-water temperature, which are controlled by high-latitude surface temperatures (Skinner and Shackleton, 2005; Lisiecki and Raymo, 2007). LR04 is a global average signal, from which regional variability was smoothed. Although LR04 is not an area-weighted stack and is biased to the Atlantic Ocean and eastern equatorial Pacific region (Elderfield et al., 2012), this record is actually one of the best syntheses of the benthic δ¹⁸O signal to study some features of the global climate system. The LR04 stack was found to be orbitally tuned to an ice-model, driven by the insolation that occurs on 21st June at 65 °N in La93(1,1) orbital solutions (Laskar et al., 1993), with corrected sedimentation rate (Lisiecki and Raymo, 2005). There was a temporal increase in the errors between La93 and La2004-La2010 orbital solutions, which becomes noticeable only after 10 Ma (Laskar et al., 2004; Laskar et al., 2011). In particular, within the 5.3 Ma of the LR04 record, the insolation phase and amplitude errors between La93 and La2004 orbital solutions are very low (Viaggi, 2018a) and can be considered negligible for the last 800 kyr.

The δ¹⁸O of benthic-planktic stack over the past 2.0 Ma (Huybers, 2007).

Antarctica proxies

The δD, CO₂, CH₄ records acquired by the European Project for Ice Coring in Antarctica (EPICA) for the last 800 kyr (Jouzel et al., 2007; Luthi et al., 2008; Loulergue et al., 2008; Bereiter et al., 2015). The EPICA time-series provided by Past Interglacials Working Group (PIWG 2016) were recalibrated based on the robust AICC2012 age model, which was built by combining glaciological inputs and data constraints, including a wide range of relative and absolute gas and ice stratigraphic markers (Bazin et al., 2013). The 0.5-kyr resampled EPICA orbital components isolated by SSA (Viaggi, 2021a) and recalibrated on the AICC2012 age model (PIWG, 2016) were applied in the present work. The EPICA stack is the mean signal among the highly covariant standardised δD, CO₂ and CH₄ SSA-filtered records (Viaggi, 2021a), and the EPICA smoothed stack (stack_s) is the Savitzky-Golay filtered signal by a 4-order least squares polynomial fitting across a moving window of 15. This is to compensate for the higher resolution of the original EPICA record. The δD proxy depends on Antarctic site temperatures derived from the ratios of water isotopes in ice cores (Jouzel et al., 1997). CO₂ and CH₄ concentrations were determined directly from measurements in the air within ice cores (Luthi et al., 2008; Loulergue et al., 2008). The CH₄ content results from changes in wetlands in the tropics and high northern latitudes (Loulergue et al., 2008; PIWG, 2016). Thus, ice core CH₄ may be viewed as a signal integrating parts of the northern hemisphere terrestrial biosphere (Loulergue et al., 2008; Luthi et al., 2008; PIWG, 2016).

Atlantic proxies

The Si/Sr and Ca/Sr ratios were used as proxies of ice-rafted debris (IRD) during the last 1.4 Ma in the IODP Site U1308 (eastern North Atlantic). The age model of Site U1308 was constructed using radiocarbon dates and oxygen isotope stratigraphy (Hodell et al., 2008).

Ca/Ti record over the past 1.4 Ma reflects the relative changes of biogenic carbonate and detrital sediment at Site U1385 (Shackleton Site) on the SW Iberian Margin. Biogenic carbonate increases during interglacial and interstadial climate states, and decreases during glacial and stadial periods (Hodell et al., 2015).

SST records from the high latitude North Atlantic ODP Site 982 located on top of the Rockall Plateau over the past 4 Ma (Lawrence et al., 2009) were studied. The sample resolution is 4.0±3.4 kyr.

SST records at North Atlantic DSDP Site 607 for the last 3.2 Ma were also considered. The age model was obtained directly from the LR04 stack (Lawrence et al., 2010).

The Atlantic benthic δ¹⁸O stacks for the last 800 kyr were studied. This is a new stack of 20 Atlantic sites of the LR04 including the new records ODP 926 and 928. Synchronisation is based on the adjusted LR04 age model to incorporate the differences between Atlantic and Pacific termination durations (Lisiecki and Raymo, 2009).

Pacific proxies

The ODP Site 846 SST record (Herbert et al., 2010) and its 1-kyr resampled orbital SSA-components (Viaggi, 2018a) were considered. This is a Plio-Pleistocene (last 4.8 Ma) alkenone-based SST record located in the Eastern Equatorial Pacific region. The record is synchronised via the δ¹⁸O values of benthic foraminifera measured in the same sediments and aligned to the LR04 stack (Herbert et al., 2010).

SST variability over the Pleistocene (last 1.5 Ma) estimated from a southern Coral Sea sediment core at site MD06-3018 (subtropical southwest Pacific) was compared with those of the sites MD97-2140 (Western Equatorial Pacific) and ODP 846 (Eastern Equatorial Pacific) (Russon et al., 2011).

The planktonic δ¹⁸O (Cheng X. et al., 2004), benthic δ¹⁸O (Ao et al., 2011), and SST (Li et al., 2011) records covering the last 4 Ma were studied at ODP Site 1143 located on the southern continental margin of the South China Sea (tropical western Pacific). The resolution of the records are: 2.4±1.8 kyr (planktonic δ¹⁸O), 2.6±1.3 kyr (benthic δ¹⁸O), 2.5±1.9 kyr (SST).

The tropical SST changes at IODP Site 1146 (last 2.2 Ma) located in the South China Sea with an age model synchronised to the LR04 δ¹⁸O stack (Herbert et al., 2010).

Pacific benthic δ¹⁸O stack for the last 800 kyr was considered. This is a new stack of 14 Pacific sites of the LR04, including the new site PC18. The synchronisation is based on the LR04 age model with adjustments to include the differences between Atlantic and Pacific termination durations (Lisiecki and Raymo, 2009).

Mediterranean proxies

The benthic δ¹⁸O records (1-kyr resolution) from the eastern Mediterranean were studied for the last 1.2 Ma from ODP Sites 967 and 968. The age model is based on tuning the elemental ratio Ti/Al against insolation, the former being dominated by the precession-related changes in northern African climate, i.e., monsoonal forcing (Konijnendijk et al., 2015).

Indian proxies

The robust dates and the highly resolved (mean 0.26±0.32 kyr) record of the Red Sea relative sea level (RSL), which approximates the global (eustatic) sea level changes (ESL), were analysed for the last 500 kyr (Grant et al., 2014). The age model is derived by synchronisation of the Red Sea dust and RSL records with the U/Th-dated δ¹⁸O of Asian speleothems. This record was partitioned in orbital SSA-components in the present work.

The tropical SST changes at ODP Site 722 covering the last 3.3 Ma located in the Arabian Sea with an age model synchronised to the LR04 δ¹⁸O stack (Herbert et al., 2010). The sample resolution is 2.0±0.9 kyr.

The planktic δ¹⁸O and SST composite records from the Sites RC11-120 and E49-18 (southern Indian Ocean) during the last 450 kyr (Hays et al., 1976) were studied.

2.2. Statistical methods

Several statistical methods have contributed significantly to the recognition of key aspects of ODH. In this work, resampling by cubic spline interpolation (CSI) at 1.0 kyr (0.5 kyr when the reference record is EPICA) was applied on time series that are processed with methodologies for which evenly-spaced records are necessary or preferable (response sensitivity, SSA, cross-spectral, Morlet wavelet filtering). CSI is a robust procedure for both up-sampling and down-sampling because the range and the number of output points are specified with data constrained from algorithm that prevents the generation of artefacts. The procedure fits a cubic spline interpolant where the knots of the spline are the X (times) of the magnitude of the data. The spline passes exactly through each data point. Endpoint conditions correspond to the not-a-knot Akima condition to prevent the wild swings observed when noisy data are fitted to cubic splines. Although resampling can generate oversampling of some records, the robustness of the CSI algorithm in preventing artefacts and the fact that the target of orbital frequencies used in this study lies well above the Nyquist’s original limits do not impose any problems related to the resolution of the original records.

The response sensitivity (R_s) to orbital forcings is a convenient way to compare the amplitude of climate proxies and nominal solutions following the work of Viaggi (2018a). Here, R_s has been applied to the EPICA and LR04 δ¹⁸O time series for the last 800 kyr. It was computed using the following equation:

Rs = (σ2resp ‒ σ2forc) / σ2forc * 100

where σ2resp and σ2forc are the variance of the standardised (0-mean; 1-standard deviation, SD) orbital SSA components (comp) of EPICA δD, CO₂, CH₄ (Viaggi, 2021a) and LR04 δ¹⁸O (Viaggi, 2018a) and the forcing nominal solutions, i.e., eccentricity La2010 (Laskar et al., 2011), obliquity and precession La2004 (Laskar et al., 2004), calculated by arbitrary time segments binned at 80 kyr. The δ¹⁸O SSA-components and the nominal solutions were resampled by CSI at the same spacings of EPICA, i.e., at 0.5 kyr (Viaggi, 2021a) to prevent data density distortion in descriptive statistics. For the purposes of this exercise, the differences in age model between EPICA and LR04 records should be negligible on large time bins. A self-sustained climate system, which is only paced by orbital forcings, would generate near zero values of R_s because the variance of the response would match the variance in forcing, whereas values that are significantly larger than zero would indicate a nonlinear reinforcement of the response signal. Negative R_s values (up to the physical limit of -100%) show a response variance lower than orbital forcing, indicating a damping mechanism of the climate system.

Factor analysis by principal component analysis (PCA) was applied to R_s data to investigate correlation patterns. PCA is a statistical methodology that attempts to identify the underlying variables (latent factors) that explain structured patterns of linear correlations within a multivariate dataset, thereby, reducing the complex data matrix into a hierarchical set of low-dimensional embedment of data points, ordered by variance. One benefit of this procedure is that it may capture some hidden meaning from the data, such as underlying physical (Steinhilber et al., 2012; Zharkova et al., 2012) or chemical (Viaggi et al., 2019) processes that are yet to be identified.

Cross-spectral analysis conducted with Hann window on evenly-spaced records were applied to investigate coherency and phase relationships between La2004 nominal solution and EPICA/SST/benthic δ¹⁸O signals in the obliquity band, and between La2004-La2010 nominal forcings and the related SSA-components from the Red Sea RSL/benthic δ¹⁸O. Analysis on detrended standardised records (0-mean, 1-SD), and δ¹⁸O/nominal precession signals were inverted to obtain the same climate polarity.

The Red Sea RSL record were partitioned into orbital components using SSA with the aim of assessing the obliquity damping on a time-series of 500 kyr long only, and to evaluate the spectral coherence of the RSL orbital components with global LR04 δ¹⁸O and nominal forcings. The SSA is an advanced method for time-series analysis used to isolate independent components based on signal strength (variance) and to the adaptive basis generated by the time series itself (Vautard and Ghil, 1989; Elsner and Tsonis, 1996; Ghil et al., 2002; Hassani, 2007; Viaggi, 2018a). Indeed, one of the most important features of the SSA is the implicit approach towards both the quantitative estimation of the signal strength and component extraction, which are fundamental task in paleoclimatology and complex signal processing (Viaggi, 2018a; 2021a). SSA requires evenly spaced time-series, and therefore, the RSL record was resampled at constant time intervals of 1.0 kyr throughout CSI in accordance to the resolution of both nominal forcings (Laskar et al., 2004; Laskar et al., 2011) and LR04 components (Viaggi, 2018a). SSA was achieved by the computation of 200-order forward-backward covariance data matrix using the singular-value decomposition procedure. Fourteen SSA orbital-components were identified spanning from short eccentricity to precession that cumulatively explain 86.1% of the original RSL variance. Signal components with similar spectral frequencies were merged into a single data component marked with the rank number of the basic components (e.g., component-1-2; 7-10 represents the sum of components 1 and 2, and from 7 to 10). The final result of this grouping process generated three RSL orbital components. The SSA-components were subsequently investigated using Fourier frequency spectrum (FFS) with data window to analyse the frequency power and test the significance of the components. The data tapering used in this study was the cs2-Hann window and the Prime Factor fast Fourier transform algorithm. The SSA-component time-series of the Red Sea RSL can be found online as supplementary material.

3. Results and Discussion

3.1. Evidence of post-MPT Obliquity Damping

3.1.1. EPICA record

R_s to orbital forcing - Quantitative analysis of the R_s to orbital forcings is demonstrated based on EPICA δD, CO₂, CH₄ and LR04 δ¹⁸O records for the last 800,000 years. Although the length of the current EPICA time series does not capture the complete temporal evolution of the MPT, some interesting features can still be observed. Viaggi (2018a) reported a sharp decrease in global δ¹⁸O and equatorial SST to near-zero values in obliquity R_s after 533 kyr. A similar post-MPT feature is observed in all EPICA obliquity δD, CO₂, and CH₄ signals since the 560–640 kyr time bin (Figure 1).

In fact, the obliquity R_s exhibits a nonlinear amplification of δD by 125%, and it reached values of ~340% and ~400% for CO₂ and CH₄, respectively at the 560-640 kyr time bin. After this amplification episode, a sharp decline in generalised obliquity R_s to near zero values (balance line) or up to ‒80%, –60% is evinced since 560 kyr. The obliquity amplification and damping effects are more significant in the CO₂ and CH₄ responses. Also the δ¹⁸O obliquity R_s shows nonlinear amplification of up to 170%, followed by a decline to near zero value or to –60%, ‒50% damped signal. The short eccentricity R_s exhibits a similar shape among signals with a dominantly strong nonlinear amplification pattern within the 720-800 kyr time bin and after 480 kyr. Indeed, the R_s in short eccentricity amplification reaches its maxima up to ~400%, 600% with an attenuation episode at the 560-640-kyr time bin, thereby, reaching values of –40%, ‒80% which correlates with the obliquity amplification maxima (Figure 1). Two amplification maxima in short eccentricity R_s may be highlighted at 400-480 kyr (δD=510%; CO₂=440%; CH₄=590%; and δ¹⁸O=620%) and 80-160 kyr time bin (δD =390%; CO₂=390%; CH₄=150%; δ¹⁸O=470%), the former being the highest expression of the climate response amplification. The precession R_s oscillates quite regularly between damped cycles till up to –90%, and amplified responses up to 100%, 200%.

In general, we observe an inverse correlation between the R_s of short eccentricity vs. obliquity. These results are in agreement with the notion of a more intense MISs after 450 kyr, compared to that from 800 to 450 kyr, which is consistent with a transition between two distinct climate states (Lang and Wolff, 2011; PIWG, 2016; Barth et al., 2018). The increase in amplitude of the 100-kyr climate cycles began ~430 kyr ago, and is known as the Mid-Brunhes Event (MBE) (Berger and Yin, 2012). It occurred within the Mid-Brunhes Oscillation (MBO) of global δ¹⁸O (Viaggi, 2018a) and EPICA δD, CO₂, and CH₄ records (Viaggi, 2021a). Indeed, the MBE corresponds to the maximum short eccentricity R_s amplification of the last 800 kyr. In summary, the analysis of the Antarctic R_s reveals that since 480 kyr, a strong amplification of the short eccentricity signals has occurred (up to ~400%, 600%) coupled with damping of obliquity responses (up to ‒80%, –60%). Table 1 shows the maximum amplitude amplification (R_s max) and the maximum amplitude damping (R_s min) values of the signals during the last 800 kyr, normalised in % for kyr. Interestingly, the R_s data show a marked asymmetry of the climate responses with a very robust amplitude magnifications and moderate damping, similar to that documented from global δ¹⁸O, and equatorial Pacific SST signals (Viaggi, 2018a).

Factor model of orbital R_s - A PCA model of EPICA R_s data was obtained to investigate certain correlation patterns that are recognisable from Figure 1. The model converges on principal component 1 (PC-1), explaining 43.8% of the total variance, whose loadings on δD, CO₂, CH₄ and δ¹⁸O R_s express a positive correlation with short eccentricity and precession (which reflects the eccentricity modulation of precession), and an anticorrelation with obliquity (Figure 2). It is interesting to note that PC-1 tends to exhibit a positive correlation with the mean of the δ¹⁸O exponential trend, supporting a link between orbitals R_s and the long-term ice volume (Viaggi, 2018a).

In particular, Figure 2b exhibits the PC-1 regression factor score and the δ¹⁸O exponential trend varying according to time bins, which reflects the increasing tendency of PC-1 factor score (i.e., short eccentricity/precession amplification vs. obliquity damping) together with the δ¹⁸O enrichment towards the Late Pleistocene (post-MPT). However, the most pronounced δ¹⁸O enrichment rate occurred earlier in subtrend IV, which includes the MPT (Viaggi, 2018a; present work). The turning point (PC-1 = 0) between the mean values of negative (high obliquity R_s vs. low short eccentricity R_s) and positive (low obliquity R_s vs. high short eccentricity R_s) PC-1 is within the time range 560-400 kyr (Figure 2b), marking the transition to MBE at ~430 kyr. These results suggest that PC-1 is a latent factor underpinned by a post-MPT anticorrelation among obliquity and short eccentricity/precession R_s, which is related to the long-term growth of the cryosphere volume.

Rescaled quantitative impact - Viaggi (2021a) estimated the quantitative impact (% variance) of orbital and millennial-scale suborbital components using EPICA δD, CO₂, and CH₄ records for the last 800,000 years. The signal strength obtained by variance is a relative concept, and is inversely proportional to the total variance of the record being investigated, and thus to its temporal extension. Therefore, for a short time series, the same component exhibits more relevant strength compared to long time series because of the former ‘defect’ in total variance owing to the lack of high order cycles or long-term trend (Viaggi, 2021a). In the present paper, the EPICA stack data were rescaled for the Plio-Pleistocene LR04 δ¹⁸O/SST orbital/suborbital variance (Viaggi, 2018a) to provide a rough estimate on a long time-frame and compare the results. Because the extremely small fraction of δ¹⁸O (0.5%) and SST (0.4%) ‘noise’ components of Viaggi (2018a) partly includes suborbital signals below the Nyquist’s condition, this fraction was included in the mean calculation base of 23.6% for rescaling EPICA stack signals. As expected in a post-MPT record, we observed that the short eccentricity is primarily responsible for the variance of the EPICA components, which was estimated to have a stack mean of 51.6% (Viaggi, 2021a). This relevant strength of the short eccentricity signal reflects the robust nonlinear amplification pattern shown by R_s analysis (Figure 1 and Figure 2). By rescaling this value on the Plio-Pleistocene variance, we obtain an estimate of 12.2%, which is very high compared to the global δ¹⁸O of 6.5% (Viaggi, 2018a). This difference should be apparent because the EPICA short eccentricity response is a post-MPT strong amplified signal, whereas the Plio-Pleistocene LR04 δ¹⁸O contains a long pre-MPT record of weak eccentricity response (Lisiecki and Raymo, 2007; Lisiecki, 2010; Viaggi, 2018a). The EPICA obliquity band accounts for a stack mean variance of 19.0%, which is almost equivalent to a Plio-Pleistocene rescaled variance of 4.5%. As opposed to the short eccentricity, this value is clearly too low compared to the LR04 δ¹⁸O value of 9.9%. This may be because of the dominance of damped signal in the EPICA post-MPT obliquity record, which is in agreement with the R_s analysis, whereas the δ¹⁸O includes a long pre-MPT record of strong obliquity response (Lisiecki and Raymo, 2007). The EPICA stack mean variance driven by precession can be estimated to be 8.4%. On the Plio-Pleistocene rescaled variance, the Antarctic precession stack correspond to 2.0%, which is close to the global δ¹⁸O value of 2.0% because even the LR04 precession signal exhibits an amplitude increase during MPT and post-MPT times (Lisiecki and Raymo, 2007; Viaggi, 2018a).

Wavelet power spectra and SSA signal structural observations - Figure 3 shows the Morlet wavelet power spectra of EPICA short eccentricity, obliquity, and precession SSA-stacks (Viaggi, 2021a). We observed a decrease in the obliquity response after ~550 kyr associated with an evident increase in both short eccentricity and precession power.

An Antarctic SSA signal structural observation related to the eccentricity-obliquity link can be seen in Figure 4. δD, CO₂, and CH₄ component-2 signals exhibit two or three low-amplitude 41-kyr obliquity peaks (glacial/interglacial) embedded in a weak ~93/75-kyr framework (Viaggi, 2021a). In contrast, the obliquity maxima associated with glacial terminations appear more pronounced. These features are similar to the shapes shown by Viaggi (2018a) for the global δ¹⁸O and equatorial SST components-3-4 during the MPT, and may be further evidence of an obliquity attenuation phenomenon. This suggests a connection between short eccentricity and obliquity, a possible observational reminiscent of the ‘obliquity-cycle skipping’ model (Huybers, 2007; Liu et al., 2008).

3.1.2. Sea-level Record of the Red Sea

Fourteen SSA orbital-components were identified in the present work from the RSL record of the Red Sea (Grant et al., 2014). The components spanned from short eccentricity to precession, and they cumulatively explain 86.1% of the original variance (Table 2). The reconstructed component of 100- kyr short eccentricity (comp-1-2;7-10) explain 61.4% of the variance. The 41-31 kyr obliquity component (comp-3-4;13-14) describes 13.1% of variance, followed by 11.7% of the 23-19 kyr precession signal (comp-5-6;11-12). The strength of these orbital RSL signals for the last 500 kyr is very similar to the corresponding bands of the EPICA records for the last 800 kyr (Viaggi, 2021a).

As in the case of EPICA, the RSL data were rescaled with respect to the Plio-Pleistocene variance to provide a rough estimate on a long-time frame and compare the results (calculation base of 23.6%), which were found to be largely similar to that of EPICA (Table 3). Indeed, as expected in a post-MPT record, the RSL short eccentricity is primarily responsible for the variance of the Red Sea record, being rescaled to 14.5%. Again, this strength is very high compared to the global δ¹⁸O of 6.5% (Viaggi, 2018a). The RSL obliquity signal is approximately equivalent to a Plio-Pleistocene rescaled variance of 3.1%. In contrast to the short eccentricity, this value is clearly too low compared to the LR04 δ¹⁸O of 9.9%, corroborating the RSL obliquity damping evidence, similar to the EPICA record. The Red Sea RSL precession signal should correspond to 2.8% rescaled variance, which is slightly amplified compared to δ¹⁸O.

3.1.3. LR04 δ¹⁸O and equatorial ODP Site 846 SST

R_s analysis – A remarkable indication of inverse coupling between obliquity damping and short eccentricity/precession amplification during post-MPT time is observed from LR04 δ¹⁸O and Site 846 SST records (Viaggi, 2018a). In this section, these observations are further developed. Figure 5 plots the Plio-Pleistocene orbital R_s (short eccentricity, obliquity and precession) of global LR04 δ¹⁸O (a, b) and equatorial Site 846 SST (c, d) as a function of the long-term mean climate state (averages of exponential-fit trend and long-term SSA-component expressing the climate background).

The values were calculated by arbitrary time segments binned at 532 kyr (data from Viaggi 2018b)1. R_s was calculated with respect to nominal forcing (eccentricity La2010, Laskar et al., 2011; obliquity and precession La2004, Laskar et al., 2004). Figure 5a-b exhibits a δ¹⁸O nonlinear amplification of the short eccentricity R_s up to ~440% towards the Mid-Late Pleistocene, similar in shape with that of precession response but with a magnitude of up to ~200%. The δ¹⁸O obliquity R_s shows a nonlinear amplification of up to 180% towards the end of the MPT, followed by a strong depletion in variance to near zero value at post-MPT (0-532 kyr bin). This R_s pattern is similar to the regional equatorial SST record (Figure 5c–d), although the values are more scattered with smaller positive magnitude of short eccentricity (it does not reach 200%) and precession amplification of up to ~100%. The SST obliquity R_s also shows strong variance depletion of up to ‒20% post-MPT.

R_s orbital factor model - PCA factor model was obtained to investigate the correlation patterns recognisable from Figure 5 by integrating the temporal R_s data of both δ¹⁸O and SST records with the long-term mean climate state (Figure 6). For this Plio-Pleistocene dataset, the model solution identified two principal components that can explain 94.6% of the total variance. PC-1 (57.7%) (Figure 6a) exhibits a strong positive loading on both δ¹⁸O and SST short eccentricity/precession R_s, which is strictly related to the long-term δ¹⁸O enrichment and long-term SST depletion. PC-2 (36.9%) (Figure 6b) shows robust positive loadings on both δ¹⁸O and SST obliquity R_s, which is related to the long-term δ¹⁸O enrichment/SST depletion.

Noteworthy, these Plio-Pleistocene factors (Figure 6c) corroborate the latent link between the increasing amplitude of all orbital climate responses and the anticorrelation between obliquity and short eccentricity with the progressive development of the Earth’s icy-state (Viaggi, 2018a), which is in agreement to the factor model of EPICA. In detail, as shown in Figure 6d, the spread of PC factor scores (namely the difference PC-2 ‒ PC-1) highlights the post-MPT anomalous depletion of the obliquity’s R_s which appears as the strongest final stage of similar and less marked PCs anticorrelation patterns since ONHG and INHG events (red bar). Indeed, the absolute factor spread exhibits two patterns of increasing value and inverse coupling (anticorrelation) among climate responses to forcings through time. They include positive spread (blue bar), showing high-obliquity associated to low-short eccentricity/precession R_s and negative spread (red bar), showing low-obliquity linked to high-short eccentricity/precession R_s. These response configurations are associated with three transition patterns of positive to negative spreads including ONHG (Transition-1, TRA-1), INHG (Transition-2, TRA-2), and MBE (Transition-3, TRA-3), the latter being characterised by high magnitude and incorporation of the MPT. Thus, the strengthening of the short eccentricity signal appears to be associated with the obliquity damped response, whose maximum expression is observed in post-MPT, and is preceded by an earlier and weaker expression straddling the ONHG and INHG. These novel results suggest the traditional notion of MPT to be the final, nonlinear transitional stage of a complex competing interaction between obliquity vs. short eccentricity forcing. It occurred under the influence of the long-term cooling that was already started during the Piacenzian (3.60–2.58 Ma), when the first weak spectral patterns of short eccentricity response are observed (Meyers and Hinnov, 2010; Viaggi, 2018a; present work). This early evidence of short eccentricity response is sometimes known as ‘premature 100-kyr cycle’ (McClymont et al., 2013) or ‘early onset of the 100-kyr-like cycle’ (Liu et al., 2008), however, it is referred by an accelerated cooling of SSTs from 1.2 to 0.9-0.6 Ma, which roughly coincides only with subtrend IV (Viaggi, 2018a) and MPT.

The R_s analysis documented by EPICA δD-CO₂-CH₄, global δ¹⁸O, and equatorial Pacific SST records support the orbital climate responses to be interpreted mainly as nonlinear feedback signals, paced by orbitals (Hays et al., 1976; Shackleton, 2000; Huybers, 2007; Lisiecki, 2010; Feng and Bailer-Jones, 2015) and damped or amplified asymmetrically in a range from −100% to +400%, +600% the forcing (i.e., as stacked signal among the primary perturbation of the forcing and the feedback responses), depending on the long-term mean climate state and cycle periodicity (Viaggi, 2018a, 2021a, 2021c). The nonlinearity is linked to the scale-dependent net balance of positive and negative feedback processes associated with carbon/water cycles and albedo mechanisms (e.g., ice-albedo, dust-albedo, vegetation-albedo) under the influence of variable initial conditions (Imbrie et al., 1993; Berger and Loutre, 1997; Shackleton, 2000; Ruddiman, 2003; Archer et al., 2004; Rial et al., 2004; Brovkin et al., 2007; Hansen et al., 2007; Herbert et al., 2010; Kohler et al., 2010; Lisiecki, 2010; Hansen et al., 2011; Russon et al., 2011; Abe-Ouchi et al., 2013; Ellis and Palmer, 2016; Shoenfelt et al., 2018; Köhler and van de Wal 2020; Viaggi, 2018a, 2021a, 2021c). In other words, orbitals pace the frequency beat of the climate response and the feedback mechanisms non-linearly transfer most of the system energy, with the least powerful climate mitigation processes. This implies that the natural climate system is dominated by positive feedbacks capable of making extremely strong amplification with responses of up to 5-7 times the forcing. This leads to the exponential behaviour of the climate system (see Figure 10 in Viaggi, 2018a), which has crucial quantitative implications also for present day climate change (Miller et al., 2010; Box et al., 2022). This interpretation may explain the energy excess ‘paradox’ of the astronomically paced signals compared to the small energy of orbital forcing, especially with respect to the eccentricity bands (Imbrie et al., 1993; Wunsch, 2004; Berger et al., 2005; Clark et al., 2006; Berger et al., 2012; Ellis and Palmer, 2016). Indeed, the maximum Plio-Pleistocene exponential growth of the LR04 δ¹⁸O short eccentricity variance (e^2.93x) suggests an enhanced sensitivity of short eccentricity forcing to nonlinear amplifying feedback mechanisms on long duration cycle during the Mid-Late Pleistocene icy-state (Viaggi, 2018a). The δ¹⁸O exponential coefficient progressively decreases for short length cycles from obliquity (e^1.88x), precession (e^1.78x), to half-precession (e^1.05x).

SSA-structural signal observations: A SSA-structural evidence of the latent link between obliquity and short eccentricity comes from LR04 δ¹⁸O component-2 (transition obliquity to short eccentricity), where the transitional change in the nonlinear dynamic can be configured as a sort of ‘switching’ of the orbital power between the 41-kyr obliquity and 73/93-kyr eccentricity cycles (Viaggi, 2018a) (Figure 7). In particular, the increasing amplitude of the short eccentricity response is coupled with an obliquity amplitude reduction during MPT and post-MPT time. Variations in amplitude are basically gradual, although nonlinear. Similar SSA-features are also recorded in the equatorial Site 846 SST (Figure 8). The transition nature from obliquity to short eccentricity is less evinced in this SST record likely because the low latitude of this site. The increasing amplitude of the short eccentricity is clearer in subcomponent-1-4. The decrease in the obliquity response is highly evident.

Further observations of SSA-structural signal for a special link between obliquity and short eccentricity come from Figure 9. Here, we can see a particular shape of the δ¹⁸O and SST component- 3-4, where two or three low-amplitude 41-kyr obliquity peaks (glacial/interglacial) are embedded in a weak ~74/92-kyr framework. Similar to those shown in the EPICA records, these shapes seem reminiscent of the ‘obliquity-cycle skipping’ model (Huybers, 2007) and corroborate the observations of Liu et al. (2008).

3.1.4. Atlantic, Pacific, Mediterranean, and Indian proxies

New additional evidence of the coupling between obliquity damping vs. short eccentricity amplification is highlighted in Figure 10 at Site U1308 (eastern North Atlantic) from the work of Hodell et al. (2008). Here, the authors depict the 41- and 100-kyr filtered components of the Si/Sr signal, which show a clear decrease in obliquity amplitude coupled with an increase in short eccentricity amplitude during the MPT and post-MPT time (Figure 10d). The onset of IRD deposition, enriched in detrital carbonate (Ca/Sr peaks), during the MIS-16 marks the reversal of the anticorrelation pattern between obliquity and short eccentricity (Figure 10b). Hodell et al. (2008) speculate that the thickness of ice-sheet and the duration surpassed the critical threshold during MIS-16 and activated the dynamical processes responsible for Laurentide ice-sheet instability in the Hudson Strait region. It should be noted that MIS-16 is at the top of the global δ¹⁸O subtrend IV, indicating the largest ice volume was reached by the Plio-Pleistocene long-term trend of the ice-sheet growth (Viaggi, 2018a).

Hodell et al. (2015) generated a continuous time series of Ca/Ti over the past 1.4 Ma, which reflect relative changes of biogenic carbonate and detrital sediment at Site U1385 (Shackleton Site) on the SW Iberian Margin (Figure 11). Noteworthy, on this site as well the orbital bandpass filters show the distinctive inverse amplitude coupling between obliquity and short eccentricity/precession components. Lisiecki and Raymo (2007) have reported that the change in δ¹⁸O glacial dynamics at ~1.4 Ma is associated with an abrupt decline in the 41-kyr power and a decrease in modulation sensitivity to obliquity. Moreover, Lisiecki (2014) found that the obliquity power of the middle deep Atlantic (2300-4010 m) benthic δ¹³C decreased by 50% at 0.6 Ma, indicating a restricted change in obliquity-driven overturning circulation. This is concurrent with a dramatic increase in the 100-kyr power of benthic LR04 δ¹⁸O. Thus, MPT likely contributed to a ∼50% decrease in the obliquity power coupled to a 100-kyr power of relevant increase in the middle deep Atlantic δ¹³C (Lisiecki, 2014). Figure 12 shows the global δ¹⁸O glacial variability over the last 2 Ma in two different multi- site stacks, i.e., benthic-planktic δ¹⁸O (Huybers, 2007) and benthic δ¹⁸O (Lisiecki and Raymo, 2005). Figure 12b exhibits a strong post-MPT increase in 100-kyr (and more moderate 22-kyr) spectral power coupled to a relevant decrease in 41-kyr power, as evinced from the benthic-planktic δ¹⁸O stack (Huybers 2007). These spectral features are similar to those of the LR04 record (Figure 12c) estimated from the work of Imbrie et al. (2011).

The reconstructions of subtropical southwest Pacific climate variability over the Pleistocene were derived by Russon et al. (2011) from coupled planktic foraminiferal δ¹⁸O-Mg/Ca measurements acquired from a southern Coral Sea sediment core at site MD06-3018. A clear shift from ∼40 kyr to ∼100 kyr modes of reconstructed glacial-interglacial SST variability was observed over the MPT when compared with those of the sites MD97-2140 (Western Equatorial Pacific) and ODP 846 (Eastern Equatorial Pacific). The wavelet power spectra in these sites (Figure 13) again highlight the same relationship between SST obliquity damped responses and the subsequent increase in short eccentricity power.

Konijnendijk et al. (2015) provided a ~1.2 Ma long benthic δ¹⁸O record from the eastern Mediterranean based on ODP Sites 967/968. Figure 14a exhibits the wavelet power spectra of the δ¹⁸O data provided by Konijnendijk et al. (2015) and Figure 14b-c shows the power spectra and time-series of the wavelet filter centred on 41-kyr cycle. Although it is less evident in the Mediterranean site, there is a slight reduction in the amplitude of δ¹⁸O obliquity, associated with an increase in power at both responses of short eccentricity (strong) and precession.

Lawrence et al. (2009) provided the first continuous, orbital-resolution SST record from the North Atlantic ODP Site 982 that is located on the Rockall plateau, which is a critical region to understand the origin of the Plio-Pleistocene ice ages because of its high sensitivity to changes of North Atlantic atmospheric pressure systems and North Atlantic Drift. This site is particularly interesting because it exhibits high-amplitude variations derived from unusual obliquity (Lawrence et al., 2009). Figure 15a shows the wavelet power spectra of the SST 41-kyr filtered component from the 1-kyr resampled data of Site 982. The plot demonstrates a relevant decrease in the 41-kyr spectral power and signal amplitude during MPT and post-MPT. Herbert et al. (2010) published the estimates of SST at four tropical sites on the basis of alkenone paleotemperature determinations. They interpreted the similarity of tropical SST changes, in dynamically dissimilar regions, to reflect ‘top-down’ forcing through atmosphere circulation and strong GHGs feedback amplification, which connects the fate of Northern Hemisphere ice sheets with global ocean temperatures. More significantly, Figure 15b exhibits the SST obliquity-damping feature by wavelet filtering and power spectra from the ODP Site 722 located in the Arabian Sea, Indian Ocean (Herbert et al., 2010).

ODP Site 1143 is located within a basin on the southern continental margin of the South China Sea (tropical western Pacific). It provides planktonic δ¹⁸O (Cheng X. et al., 2004), benthic δ¹⁸O (Ao et al., 2011) and SST (Li et al., 2011) records covering the last 4 Ma. Figure 16 shows the wavelet power spectra of the 41-kyr filtered components over the last 2 Ma from the 1-kyr resampled data of Site 1143. A significant decrease in 41-kyr spectral power and signal amplitude is evinced during MPT and post-MPT period.

Finally, it is very interesting to observe evidence of obliquity damping from the historical work of Hays et al. (1976), which is based on composite records from Sites RC11-120 and E49-18 of the southern Indian Ocean (Figure 17). Here, the 40-kyr components of both planktic δ¹⁸O and SST records exhibit a clear damping trend in amplitude during the last 450 kyr.

These relevant observations from global and regional records (Antarctic, North Atlantic, Atlantic - Iberian margin, Mediterranean, western and eastern tropical Pacific, Pacific - South China Sea, Indian

- Red Sea, Indian - Arabian Sea, southern Indian) of different climate-related proxies show that the coupling between obliquity-damping and short eccentricity-amplification is a widespread feature of the Mid-Late Pleistocene climate system, which could be the missing link in the MPT origin. There is a substantial variability in the intensity of obliquity damping at some single-site records (e.g., low intensity in Mediterranean Site 967/968) likely due to the different sensitivity to obliquity of ocean basins and site locations (e.g., latitude, circulation). Understanding this phenomenon is crucial because the obliquity attenuation could have contributed to the strengthening of the short eccentricity response, favouring a long-lived feedback amplified ice-growth.

3.2. Long-Term Cooling sets Boundary Conditions for Glacial/Interglacial cycles

This section discusses the role of the long-term trend in the MPT debate in relation to the results highlighted above. Viaggi (2018a) investigated the records of the global LR04 δ¹⁸O stack (Lisiecki and Raymo, 2005), the equatorial ODP Site 846 SST (Herbert et al., 2010), and the global ΔSST stack (Martinez-Boti et al., 2015) using SSA. The study found similar trend components in all records, which explain ~76% of the Plio-Pleistocene variance and appear to be related to the long-term pCO₂ proxies (Hönisch et al., 2009; Seki et al., 2010; Bartoli et al., 2011; Martinez-Boti et al., 2015). Statistical analysis suggests that such long-term trends significantly modified the mean climate state with respect to orbitals through time (Viaggi, 2018a), thereby hinting the background forcing to be relevant in setting boundary conditions for orbital climate responses (Berger et al., 1999; Bintanja and van de Wal, 2008; Huybers and Tziperman, 2008; McClymont et al., 2013). These evidences of non-regional cooling trends are well supported by the work of McClymont et al. (2013) on Pleistocene global SST records. The study showed that the variability of SSTs is superimposed upon a long-term cooling trend in oceanographic systems spanning the low- to high-latitudes and is accompanied by evolving pCO₂, abyssal ocean ventilation, atmospheric circulation, and/or dust inputs to the Southern Ocean. The ODP Site 982 (North Atlantic) and the ODP Site 1143 (southern South China Sea) both exhibit a remarkable trend of SST/δ¹⁸O cooling/enrichment over the last 4 Ma (Lawrence et al., 2009; Li et al., 2011). The close link between long-term cooling and the coupling between amplification of the short eccentricity vs. damping of the obliquity responses is corroborated by the results of the Plio-Pleistocene R_s orbital factor model shown in Fig.s 5 and 6. The long-term cooling documented by the global LR04 δ¹⁸O stack is characterised by four step-wise subtrends (Viaggi, 2018a), where terminal accelerations configure a mild curvilinear shape that is broken by slope changes representing four thresholds (TH1 to TH4) of mean climate state variation (Figure 18). The subtrend I occurs in the time interval of 5.33-3.30 Ma (Zanclean-Early Piacenzian) and culminates with the ONHG through a mild rate change of δ¹⁸O (linear equation Y_I). This subtrend is associated with the TRA-1 (Figure 6d), where the first weak response of the short eccentricity occurs straddling to the ONHG (Viaggi, 2018a). Subtrend II is Piacenzian in age, and exhibits the first strong acceleration of the changing rate, which terminates with the INHG (linear equation Y_II). It is followed by subtrend III with intermediate rate of variation till up to 1.4 Ma (Gelasian-Calabrian p.p.). Subtrends II and III include TRA-2, where the second early manifestation and the more consistent short eccentricity response are observed across the INHG. The beginning of the MPT plots within subtrend IV (1.4-0.65 Ma), which marks the second highest rate of change (linear equation Y_IV) during the Plio-Pleistocene long-term cooling trend and ice- volume growth. The end of the subtrend IV (MIS-16) indicates that the largest global ice-volume was reached during the Cenozoic cooling (Zachos et al., 2001, 2008), and marks the final stage of the MPT (~0.7 Ma). Interestingly, subtrend IV is associated to accelerated rate of erosion/weathering cycles of the cratons surrounding the North Atlantic (Yehudai et al., 2021). After MIS-16, the cooling trend was broken by a wide swing of global δ¹⁸O depletion and temperature recovery (SST, EPICA δD), which is characterised by a relative stasis of the ice growth and temperature cooling during the last 600 kyr (MBO). These are associated with a mild midterm GHGs recovery (Viaggi, 2018a, 2021a). TRA-3, which includes MPT, occurs across subtrends III/IV and MBO, and shows the most relevant competing interaction between obliquity vs. short eccentricity responses under the influence of the long-term cooling (Figure 6). The strengthening of the short eccentricity signal appears to be associated to the obliquity damped response (early manifestations across ONHG and INHG), whose maximum expression is observed during post-MPT (i.e., MBO time).

Although the proxy-based evidences for long-term cooling that are associated with increased global ice volume appear to be reasonably credible, it is unclear how this background climate affected the MPT. Changes in diverse environmental parameters suggest that glacial climate boundary conditions evolved across the MPT, and may have altered climate sensitivity to orbital forcing by placing pre-existing ice-sheets closer to some threshold of climate-ice sheet response (McClymont et al., 2013). Berger et al. (1999) modelled the Northern Hemisphere ice-sheet volume to insolation forcing over the last 3 Ma. Under a warmer climate induced by higher atmospheric pCO₂, ice sheets never grow large enough during insolation minima such that they cannot survive any subsequent moderate insolation maxima, resulting in obliquity-related dominant cycles. Under cooler climate conditions, ice sheets survive through moderate insolation maxima and deglaciate entirely only under maximum insolation forcing (high eccentricity, high obliquity, boreal summer at the perihelion). Modelling results by Bintanja and van de Wal (2008) suggest that the gradual emergence of the 100-kyr cycles can be attributed to the increased ability of the merged North American ice sheets to survive insolation maxima and reach a continental-scaled size. An important feature is the increased inability of the system to reach full interglacial levels at the timescales of obliquity and precession, in particular in North America, where continuing cooling enabled ice sheets to overcome insolation maxima at ever-greater volumes (Bintanja and van de Wal, 2008). However, the larger the ice sheet grows and extends towards lower latitudes, the smaller is the insolation required to make the mass balance negative (Abe-Ouchi et al., 2013). Cooling of climate is expected to be analogous to increasing the threshold for ablation in the insolation integrated over the summer period (summer energy) (Huybers and Tziperman, 2008). Summer energy suggests that an ice sheet will be the most sensitive to obliquity when its ablation zone exists at high latitudes with a low insolation threshold (warm climate background). However, the model produces mostly linear obliquity response, and fails to produce trends towards greater amplitude (Huybers and Tziperman, 2008). This is a common problem with models that fail to accurately reproduce the observed records (Huybers, 2011; Imbrie et al., 2011; Berends et al., 2021).

The Plio-Pleistocene long-term cooling, and the peculiar geographic setting of Greenland and North America, have certainly led to more favourable boundary conditions for snow preservation at high latitudes. Likely, the key role of the long-term cooling in MPT is to increase the scaling effect of feedback mechanisms in response to orbital forcing, and the resulting increase in the amplitude of climate responses. However, the long-term cooling cannot explain alone the obliquity-damping observations from climate proxies in the context of an increase of the obliquity nominal forcing during Mid-Late Pleistocene, as shown in the next section.

3.3. Amplitude relationships between Orbital Forcings and Proxies

A common feature of the anticorrelation among climate proxies is the unexpected obliquity amplitude reduction coupled with an increasing amplitude of the short eccentricity response because this characteristic is not found in the nominal solutions over the last 1.2 Ma (Laskar et al., 2004; Laskar et al., 2011) (Figure 19). In fact, the amplitude trends and the power spectra of the orbital solutions appear to have inverse relationship with those of the related climate proxies: obliquity nominal-increase vs. proxy-decrease and short eccentricity/precession nominal-decrease vs. proxy-increase. These relationships are reflected in R_s data (Viaggi, 2018a; this study). On the other hand, there is no univocal relationship between the strength of the astronomical forcing and the intensity of the related climate response (PIWG, 2016). Although the amplitude of the climate response should be theoretically consistent with that of forcing, the amplification or damping in relation to the net balance of positive and negative feedback mechanisms, and the effect of the initial conditions, may alter this linear connection (Shackleton, 2000; Archer et al., 2004; Brovkin et al., 2007; Bintanja and van de Wal, 2008; Kohler et al., 2010; Lisiecki, 2010; Hansen et al., 2011; Ellis and Palmer, 2016; Shoenfelt et al., 2018; Viaggi, 2018a, 2021a). Since obliquity is a relevant forcing to high latitudes ‘ice-killing’, it is difficult to understand what ordinary negative feedback mechanism might have dampened the climate response of obliquity in the context of increasing amplitude of nominal forcing during the last 1.2 Ma.

In the following section, we have critically reviewed the theoretical constrains required by ODH for obliquity-oblateness feedback relevant to the uncertainties in the obliquity phase lag estimates and oblateness changes.

3.4. Why does Obliquity’s Response Damping?

A crucial characteristic of the damping of climate proxies in the obliquity band is its symmetry, which suggests attenuation of the forcing during both interglacial (high obliquity) and glacial (low obliquity) stages. ODH to achieve this bipolar damping by obliquity–oblateness feedback, theoretically requires that a negative secular change occurs during the obliquity maxima (interglacial) and a positive change takes place during obliquity minima (glacial), based on the model of Levrard and Laskar (2003) (Figure 20). Specifically, mandatory constrains are low obliquity lag response of the ice-sheet (<44°; <5.0 kyr), and positive (interglacial damping)/negative (glacial damping) changes of oblateness.

3.4.1. Remarks on Obliquity phase lag and Temperature/Ice-volume proxies

Response phase lags are difficult to estimate (e.g., age models) and should be compared within the same time interval with preferably evenly spaced records and the same spectral analysis parameters. Conceptually, the benthic δ¹⁸O measures changes in global ice-volume and deep-water temperature, which are controlled by high-latitude surface temperature (Skinner and Shackleton, 2005; Lisiecki and Raymo, 2007). δ¹⁸O changes are globally assumed to be synchronous within the mixing time of the ocean (~1 kyr) because the δ¹⁸O composition of foraminiferal tests is primarily controlled by the storage of ¹⁸O-depleted water in ice sheets (Lisiecki and Raymo, 2009). However, few studies have revealed that diachronous changes in deep water temperature between Atlantic and Pacific basins produce significant lags in benthic δ¹⁸O (Skinner and Shackleton, 2005; Waelbroeck et al., 2006; Lisiecki and Raymo, 2009). Moreover, the benthic δ¹⁸O records from different depths in the same ocean might have experienced diachronous responses. Because the benthic δ¹⁸O signal is a composite of these two different processes, the net phase of the δ¹⁸O signal reflects the variable response time of the ice-volume and deep-water temperature components across ocean basins and site locations (Skinner and Shackleton, 2005; Lisiecki and Raymo, 2009; Herbert et al., 2010; Lawrence et al., 2010). This may result in a lag bias in the benthic δ¹⁸O data, affecting phase shift estimates of the ice-volume. To overcome this problem, other proxies may be considered to be not affected by deep-temperature bias, and the benthic δ¹⁸O deep-temperature bias may be estimated. The EPICA stack_s better approximates the global benthic δ¹⁸O signal (Figure 21), the latter being intrinsically related to the global average surface temperature (Lisiecki and Raymo, 2007). There is a striking similarity in the glacial-interglacial variability between EPICA stack_s and LR04 δ¹⁸O, even though amplitude offsets exist possibly resulting from the deep-temperature component of benthic δ¹⁸O. Indeed, the cross-correlation function (CCF) indicates a Pearson correlation at 0-lag of 0.87, which is the highest among the single EPICA δD, CO₂ and CH₄ filtered records that are in the range 0.80-0.83. The CCF exhibits the highest coefficient (0.91) at ‒5 lag number, equivalent to a δ¹⁸O ‘bulk’ lags of ~2.5 kyr. Modelling results of Yin and Berger (2012) indicate that the variations of the annual mean temperatures averaged over the Earth and over the southern high latitudes were dominantly controlled by GHGs, whereas over the northern high latitudes, insolation played a key role. These differences between the Arctic and Southern Ocean are mainly due to their different geographical configuration, which leads to more climate responding to local and seasonal insolation forcing over Arctic and more to global and annual GHGs forcing over the Southern Ocean (Yin and Berger, 2012). Moreover, simulated changes of polar (Arctic and Antarctic) temperatures are strongly related to changes in simulated global temperatures, confirming that ice-core-based reconstructions provide quantitative insights on global-scale temperature variations (Masson-Delmotte et al., 2006). Hence, the EPICA stack may be considered a proxy of the global atmospheric temperature averaged by GHGs and related to the NH temperature that approximates the global ice-volume and deep-water temperature δ¹⁸O signal.

Studies on Greenland and Antarctica indicate a fast response of the cryosphere mass balance and sea-level rise to recent changes in atmospheric and sea surface temperature (Cuffey and Clow, 1997; Lemke et al., 2007; Hanna et al., 2008; Lawrence et al., 2008; Ettema et al., 2009; Tarnocai et al., 2009; Miller et al., 2010; Yin et al., 2011; Pattyn et al., 2018; King et al., 2020; Box et al., 2022). The cryosphere (which consists of snow, river and lake ice, sea ice, glaciers and ice caps, ice shelves and ice sheets, and frozen ground) is intricately linked to the surface energy budget, the water cycle, sea level changes, and the surface gas exchange (Lemke et al., 2007). The recent decline in ice mass are correlated with rising surface air temperatures, especially in the north of 65 °N (Lemke et al., 2007). The current change in Arctic temperature has been consistently exceeding the average temperature of the Northern hemisphere by a factor of 3-4, suggesting that Arctic warming will continue to significantly exceed the global average with concomitant reductions in terrestrial ice masses and increasing rate of sea level rise (Miller et al., 2010; Box et al., 2022). Since 1990, the Greenland Ice Sheet surface mass balance has decreased rapidly, which has been mainly generated by increased melting and runoff induced by global warming (Hanna et al., 2008; Ettema et al., 2009). The inverse correlation between accumulation rate and temperature in the mid and late Holocene suggests that the Greenland Ice Sheet is more prone to volume reduction in a warmed climate than previously thought (Cuffey and Clow, 1997). Also Yin et al. (2011) reported that owing to a mid-range increase in the atmospheric GHG concentrations, the subsurface oceans surrounding the Greenland and Antarctica ice sheets have warmed substantially at depths of 200-500 m compared with the changes observed till date. The observed acceleration of the outlet glaciers and ice flows in Greenland and Antarctica is closely linked to ocean warming in the vicinity of the ice sheets (Yin et al., 2011). The late 20th- century glacier wastage was likely a response to post-1970 warming. The strongest mass losses per unit area of mountain glaciers have been observed in Patagonia, Alaska, northwest USA, and southwest Canada (Lemke et al., 2007). Also the temperature at the top of the permafrost layer has increased by up to 3 °C since the 1980s in the Arctic region (Lemke et al., 2007). The permafrost degradation has led to changes in the characteristics of land surface and drainage systems, and most importantly, it represents an extremely large carbon pool stored in the form of peat and methane in circumpolar areas (Tarnocai et al., 2009). Lawrence et al. (2008) suggest that a rapid melting of Arctic Sea ice may initiate a feedback loop that rapidly melts Arctic permafrost, triggering further warming feedback owing to the release of methane and carbon dioxide. Modelling results on Last Glacial Maximum suggest that, once a large ice sheet is established, a moderate increase in insolation is sufficient to trigger a negative mass balance, leading to an almost complete retreat of the ice sheet (and sea-level equivalent of ice-volume changes) within several thousand years. This fast retreat of ice-volume is governed mainly by rapid ablation due to the lowered surface elevation resulting from delayed isostatic rebound (Abe-Ouchi et al., 2013). Although current climate changes are probably more rapid than those observed in the geological record, these studies demonstrate the high resilience of the climate system to temperature changes (further insights into this topic are shown in the next section). Considering this documented rapid response of the cryosphere volume and sea-level to surface temperature changes, the EPICA stack may approximate the ice-volume changes without the lag bias of the benthic δ¹⁸O, the latter likely represents the most lagged signal in the climate chain (forcing > surface temperature > GHGs > ice- volume > ESL > bottom water temperature). It follows that phase lags averaged among surface temperature proxies (EPICA stack, SSTs) is likely to be a better approximate of the lag of the ice-volume.

Next, we consider phase lag data from literature. According to Levrard and Laskar (2003), the climate friction impact was positive and most likely negligible over the last 3 Ma. However, it is based on the model assumption that the ice sheet obliquity response lag from benthic δ¹⁸O records is 8 kyr (70°) or within the range of 6-10 kyr (50°-90°) considered plausible for Plio-Pleistocene glaciation, and positive change of oblateness (Figure 20). Lisiecki and Raymo (2005) have reported that the benthic δ¹⁸O lag in obliquity reaches ~6.8 kyr during the Pleistocene. Viaggi (2018a) estimated a benthic δ¹⁸O obliquity lag of 7.4 kyr post-MPT, which is close to the value considered by Levrard and Laskar (2003). Also Imbrie et al. (1993) estimated a high obliquity lag (9.1 kyr) of the planktonic δ¹⁸O stack over nominal forcing during the last 400 kyr. In contrast, the obliquity lag estimated by Hilgen et al. (1993) provided a lower value of 5.6 kyr with respect to the Late Pleistocene if the phase relation with precession is kept constant. Konijnendijk et al. (2015) provided spectral cross-correlation data between their 1.2 Ma long benthic δ¹⁸O record from ODP Sites 967/968 and 65 °N summer insolation, and suggested an average lag of 5.5±0.8 kyr for the obliquity. In cores V3097/607 and RC11-120/E49-18 (44 °S), the planktic δ¹⁸O stack is reported to lag SST records by 1.5-1.6 kyr over the last 400 kyr (Imbrie et al., 1993). The SST obliquity component of the equatorial Pacific Site 846 exhibits a lag of 2.7 kyr versus nominal solution post-MPT start (Viaggi, 2018a). Viaggi (2021a) provided an average EPICA obliquity lag of carbon GHGs at 2.6-kyr (CO₂ 2.8-kyr; CH₄ 2.4-kyr) with respect to the atmospheric temperature. Lisiecki (2014) reported obliquity phase data of climate responses and circulation proxies for the last 1.5 Ma. While the obliquity phase lag appears to be high for both benthic δ¹⁸O and Middle deep Atlantic δ¹³C (6.7 kyr), the SST proxies from North Atlantic Site 607 and U1313 (Site 607 redrilled as Integrated Ocean Drilling Program U1313) exhibit values <4.9 kyr, with an average of 3.5±1.4 kyr. Lawrence et al. (2010) interpreted the slight obliquity phase lead at North Atlantic Sites 982 (~20°; ~2.3 kyr) and 607 (~10°; ~1.1 kyr) of SST records relative to benthic δ¹⁸O throughout the late Pleistocene as an indication of the relative response times of these components to insolation forcing, where the high inertia of the ice sheets lag the rapid response of the ocean surface. These phase data, although not referred to nominal solutions, agree with Herbert et al. (2010), whose results from four tropical ODP Sites indicate that SST leads by approximately 2 to 5 kyr the glacial cycles as recorded by benthic δ¹⁸O. Nevertheless, a surficial proxy that leads benthic δ¹⁸O may not necessarily indicate that the climate response occurred before ice-volume change because of bias in benthic δ¹⁸O deep-temperature lag (Lisiecki and Raymo, 2009).

New cross-spectral data in the obliquity band acquired from different climate proxies (EPICA, SSTs and benthic δ¹⁸O) and obliquity nominal solution (Laskar et al., 2004) for the last 800 kyr are shown in Table 4. Despite potential age-model discrepancies, the EPICA stack exhibits a lagging pattern of 38° (4.2 kyr), lower than the theoretical threshold of 44° (5.0 kyr). Incorporating the SSTs data, the average delay can be estimated to be 33°±15°, which is equivalent to 3.7±1.7 kyr. Also considering the benthic δ¹⁸O, the mean phase lag is 48°±6° (5.3±0.6 kyr), which is very close to the theoretical threshold and significantly lower than the range of 50°-90° (6-10 kyr) that has been considered for Plio-Pleistocene glaciations. Therefore, assuming a rapid response of the ice-volume to surface temperature changes, the lowering of the obliquity phase lags compared to the ODH theoretical threshold seem to be more credible.

3.4.2. Observations of Orbital phase lags between δ¹⁸O and Red Sea RSL records

In this section, phase lags are explored further by considering the orbital relationships between the global benthic δ¹⁸O and the Red Sea RSL vs. nominal forcings during the last 500 kyr. According to Grant et al. (2014), the RSL approximates global (eustatic) sea level and is broadly similar to benthic δ¹⁸O variations; however, phase and amplitude offsets exist because of the deep temperature component from benthic δ¹⁸O (Skinner and Shackleton, 2005; Lisiecki and Raymo, 2009; Lawrence et al., 2010; Elderfield et al., 2012). Marine δ¹⁸O record contains a globally synchronous glacio-eustatic component and provides a useful proxy for sea level changes (Shackleton, 1967; Skinner and Shackleton, 2005). The SSA decomposition of the RSL record carried out in the present work (Table 2) allows us to evaluate the phase relationships and the coherency of the signals using cross-spectral analysis in the orbital bands (Figure 22, Table 5). Assuming a certain margin of error due to different age models, all signals exhibit good coherency (0.75-0.94) and phase lags with respect to nominal solutions that range from 2.7 to 7.65 kyr. The coherency between δ¹⁸O and RSL in all orbital bands is very high (0.83-0.94). This high cross-coherency among orbital forcings, global δ¹⁸O, and Red Sea RSL records corroborate that the RSL approximates very well the glacio-eustatic sea-level fluctuations linked to the ice-volume and paced by orbital forcings, also in the short eccentricity band. In terms of the phase relationships, the RSL lags nominal obliquity by 5.65 kyr, whereas δ¹⁸O lags obliquity by 7.65 kyr, resulting a benthic δ¹⁸O delays of 2 kyr compared to that of RSL. In the precession band, the RSL lags by 2.92 kyr and the δ¹⁸O lags by 3.96 kyr, resulting in benthic δ¹⁸O delays of ~1 kyr with respect to RSL. The RSL and the δ¹⁸O lags short eccentricity by 3.5 kyr and 2.72 kyr, respectively, resulting a δ¹⁸O lead of 0.78 kyr compared to RSL.

Notably, the global benthic δ¹⁸O lags RSL obliquity by 2.0 kyr and precession by 1.0 kyr, which are counterintuitive because ESL follows the changes in ice-volume (Skinner and Shackleton, 2005; Miller et al., 2010; Box et al., 2022). In case there is no systematic errors in the age models, the data suggests a delay bias that could be attributed to the benthic δ¹⁸O deep-water temperature component (Skinner and Shackleton, 2005; Waelbroeck et al., 2006; Lisiecki and Raymo, 2009; Lawrence et al., 2010). In contrast, the short eccentricity exhibits opposite relationship, where benthic δ¹⁸O leads RSL by 0.78 kyr, which seems less credible considering the long periodicity and the bias of the deep temperature component. This result could be an artefact owing to the averaging of the lead-lag phase variations over a small number of short eccentricity peaks (9-10) during the last 500 kyr (Figure 22a). If the estimate of the δ¹⁸O obliquity deep-water temperature bias is correct, then the benthic δ¹⁸O phase lags in the obliquity bands (Table 4) should be reduced by ~2 kyr. Accordingly, the benthic δ¹⁸O obliquity mean phase lag of 5.3 kyr could be corrected to 3.3 kyr, which is close to the mean EPICA stack/SST of 3.7 kyr. The slight advance of the average ‘unbiased’ benthic δ¹⁸O over surface temperature is likely to be within the error margins of these data. Finally, the obliquity lag of the ESL of 5.65 kyr during the last 500 kyr corroborates previous considerations of a rapid ice-volume response to temperature variations also at the time scale of the geological record.

In the next section, we have considered changes in Earth’s oblateness from present-day satellite data that are very useful, albeit worrying, in evaluating the sensitivity and the direct effect of mass redistribution on the Earth’s oblateness.

3.4.3. Changes in Earth’s Oblateness

Here, we have discussed how polar ice-loading and global water-mass redistribution might affect the changes in the Earth’s dynamic oblateness (J₂). Oblateness refers to the Earth’s ellipsoidal shape, i.e., its equatorial radius is ∼21 km larger than its polar radius, and most of this oblateness arises due to Earth’s rotation, which pushes its mass out towards the equator (Chao, 2006; Nerem and Wahr, 2011).

In particular, J₂ changes are determined by both glacial isostatic adjustment (GIA) and the gravitational redistribution of water in the oceans and atmosphere (W). GIA describes the response of the solid Earth, the gravitational field, and the oceans to mass changes of the global ice sheets (Whitehouse, 2018). A commonly studied component of GIA is ‘postglacial rebound’, which specifically relates to the uplifting of the land surface following ice melt, resulting in a decrease in J₂ (Mitrovica and Peltier, 1993; Cheng and Tapley, 2004; Cheng et al., 2013; Seo et al., 2015; Nakada and Okuno, 2017; Whitehouse, 2018). A negative trend indicates that Earth is becoming less oblate. During deglaciation stage, water-ice masses are transferred from the polar regions to the global oceans, thereby increasing the polar moment of inertia and hence the J₂ (Levrard and Laskar, 2003; Chao, 2006). Simulations conducted by Bills (1994) suggest that the mass transport from the oceans to the polar ice sheets during Late Pleistocene glaciations can significantly change the gravitational oblateness of the Earth. The author concludes that it may be plausible that the MPT involves an obliquity–oblateness feedback, reflecting the evolution of the orbital-rotational-climatic system away from obliquity dominated oscillations. The present-day variations in temporal gravity estimated by Earth satellite observations offer the opportunity to measure the ongoing changes in J₂. Nerem and Wahr (2011) use present-day gravity variations from the gravity recovery and climate experiment (GRACE) satellite mission to investigate changes in a 34-year time series of the Earth’s oblateness observed by satellite laser ranging. The authors observed that from 1975 to 1997, the oblateness variability exhibits a secular negative trend, which is believed to be due to postglacial rebound. Noteworthy, their results suggest that since 2002, the ice loss from Greenland and Antarctica has remained the dominant contributor to the current GIA-corrected positive J₂ trend, which apparently began sometimes during the 1990s. Similar results are reported by Seo et al. (2015) over the period 1979- 2010, which established that the variations in ice mass in Greenland and Antarctica have remained the dominant cause of the observed positive J₂ trend since 1998. Interestingly, the residual J₂ variations show a clear oscillation near 10-11 yr, suggesting some connection with solar activity (Seo et al., 2015), and thus high sensitivity of the oblateness changes. Cheng et al. (2013) reported an evident deceleration in the rate of the GIA-induced decrease in J₂ over the last few decades, and this is likely due to changes in the rate of the global mass redistribution from melting of the glaciers and ice sheets, as well as mass changes in the atmosphere and ocean. Also Nakada and Okuno (2017) estimated the total positive J₂ change of Greenland-Antarctic ice-sheets and mountain glaciers, and attested the change to recent melting models. Noteworthy, Sun et al. (2019) and Sun and Riva (2020) proposed an approach to simultaneously estimate J₂ temporal variations based on gravity observations from the GRACE satellite mission. One of the main advantages of their methodology is that it allows to separate the quantitative contributions of individual sources to the total signal. They found that the secular trend is mainly generated by 1) ice sheet melt (primarily controlled by Greenland, mountain glaciers and Antarctica, respectively), which has a positive J₂ change, and 2) the ongoing solid Earth response to past glaciation (GIA component), which has a negative trend. In particular, the trend due to ice sheet melting during 2002-2017 has a larger absolute value, such that the overall trend is rising, and the Earth is currently becoming flatter. More specifically, by using the J₂ estimates from Table 2 of Sun and Riva (2020), we calculated the contribution of the J₂ water layer to be ∼70% over the total absolute J₂ signal in 13 years only (2003-2016). A recent study on the last interglacial (MIS-5e) found that colder North Atlantic surface water and weaker Atlantic overturning characterised the early last interglacial (126 kyr) compared to the late last interglacial episode (122 kyr) (Govin et al., 2012). Also the δ¹⁸O values from the Greenland ice cores shows a multimillennial lagged response to summer insolation during the early part of MIS-1 and MIS-5e interglacials, which may be explained by changes in ice sheet elevation and the impact of meltwater on ocean circulation (Masson-Delmotte et al., 2011). These results from past interglacial stages suggest that both insolation and ice sheet melting have to be considered to reproduce the climatic pattern, and especially points out a relevant role of meltwater flux in interglacial dynamics. In particular, the remarkable result of the cross-spectral analysis of the ESL of the Red Sea (Figure 22, Table 5) suggests a rapid and coherent oscillation in the obliquity band even of the water layer component of the Earth’s oblateness.

These studies are outstanding because they indicate a high sensitivity of the Earth to variations in oblateness that is induced by fast response of the cryosphere mass balance to surface temperature changes. In addition, they highlight a concurrent slow-moderate negative (GIA rebound) and fast-robust positive (water-mass redistribution) J₂ change in the recent postglacial and global warming context, resulting in a J₂ positive net change. The secular obliquity change depends only on the amplitude and on the phase lags of the oblateness variations in the obliquity band (Levrard and Laskar, 2003). Thus, it can be speculatively affirmed that during the Plio-Pleistocene long-term glacial trend, it was the rapid and quantitatively dominant water-mass redistribution component that controlled the J₂ net-change in the obliquity band. This water-mass redistribution component, which is reasonably related to the glacio-eustatic variations of the sea level, could be a crucial element in determining J₂ positive net change during obliquity maxima. In this sense, the strong coherence between the Red Sea ESL and the global δ¹⁸O in the obliquity band (Table 5) is very significant. Also the good congruence between the ESL and EPICA in terms of the rescaled variance during the last 500 kyr (Table 3) suggests similarities in amplitude patterns. It is therefore plausible to assume that the chaining ‘forcing > surface temperature > GHGs > ice-volume > ESL > J₂ water mass redistribution’ in the obliquity band during the Plio-Pleistocene was dominant over the solid Earth’s response with an amplitude pattern similar to that of the SST-δ¹⁸O, characterised by transient phases of progressive amplitude increase linked to the icy-state (Figs. 6 and 18). Following the work of Sun and Riva (2020) and assuming obliquity phase lag <5 kyr, we can schematically approximate the effect of the obliquity-oblateness feedback as follows:

High ε (interglacial): ΔJ₂ = –GIA₃₀ +W₇₀ ≈ +W ⟹ –Δε ⟹ interglacial damping,

Low ε (glacial): ΔJ₂ = +GIA₃₀ –W₇₀ ≈ –W ⟹ +Δε ⟹ glacial damping,

where, ε refers to obliquity; ΔJ2 is the net change of oblateness; GIA₃₀ reflects slow-moderate GIA component of oblateness; W₇₀ is the fast-robust water component of oblateness; and Δε is the secular change of obliquity. It should be considered that a strong variation of the obliquity is probably not necessary owing to the exponential behaviour of climate responses caused by feedback mechanisms (Shackleton, 2000; Ruddiman, 2003; Herbert et al., 2010; Kohler et al., 2010; Lisiecki, 2010; Russon et al., 2011; Shoenfelt et al., 2018; Viaggi, 2018a, 2021a; Köhler and van de Wal, 2020). The pre-MPT amplification of the obliquity R_s (δ¹⁸O, SST) with the development of the icy-state (Figure 6d, blue bars) might have led to the increase in the amplitude of both the glacio-eustatic oscillation and the oblateness change related to the dominant and fast water mass redistribution component. This possible step-wise increase in amplitude of the obliquity-paced J₂ net change within each transition interval (TRA-1, 2, 3), along with obliquity phase lags which may also be below the theoretical threshold of 44°, might have reached J₂ sensitive values. This might have induced secular changes of obliquity, each culminating in a ‘reversal’ of obliquity weakening that favoured the strengthening of the short eccentricity response by feedbacks amplification (Figure 6d, red bars). In an icy-state background, the most critical stage of this phenomenon would be during obliquity maxima (interglacials) with the consequent weakening of the obliquity ice-killing and strengthening of the ~100-kyr long-life feedbacks-induced ice growth. The maximum expression of this mechanism would occur during the TRA-3, which is associated with the maximum development of the icy-state (subtrend IV) and strong amplification of the obliquity response of up to the MPT end (TH4). TRA-1 and TRA-2 would be early precursor stages of this complex mechanism that was initiated as early as during the ONHG.

3.5. Earth’s eccentricity and the 100,000-year issue

Milankovitch’s theory describes climate cycles controlled by high-latitude NH summer insolation which results from the collective effects of changes in the Earth’s movements (eccentricity, obliquity, precession) over thousands of years. Although broadly accepted, some specific observations are not explained by the theory. For a discussion, see Ellis and Palmer (2016) and PIWG (2016). Here we focus our discussion on the role of the short eccentricity which is closely related to the MPT conundrum. Since the historical work of Hays et al. (1976), the 100,000-year issue derives by the observation that 65°N mean monthly insolation (21 june-20 july) is dominated by precession (24-22-19 kyr, 85% power) and secondarily by obliquity (41 kyr, 15% power), while eccentricity is negligible (Figure 23a). The frequency spectrum does not change significantly if we consider the 85°N mean monthly insolation where we observe an increase in the obliquity power (23%) over precession (77%) (Figure 23b).

Indeed, the eccentricity controls primarily the amplitude of the precession cycle (Berger et al., 1993; Ellis and Palmer, 2016). However the eccentricity, through altering the mean distance from the Earth to the Sun, is the only orbital parameter which change the Earth’s mean annual insolation (MAI), although the total energy received is very small quantitatively (Berger et al., 1993; Laskar et al., 2004; Ellis and Palmer, 2016) (Figure 23c). As is well known, variance analysis of climate proxies such as EPICA shows an opposite picture respect to the nominal forcings with the short eccentricity dominating (51.6%), followed by obliquity (19%) and lastly, precession (8.4%) (Viaggi, 2021a). Hence, the main paradox of the Milankovitch theory concern 1) the origin of the short eccentricity glacial periodicity (which is negligible on high-latitude NH summer insolation), and 2) the strength of the short eccentricity component in climate proxies with respect to the forcing. Concerning the latter, Hays et al. (1976) wrote ≪… an explanation of the correlation between climate and eccentricity probably requires an assumption of nonlinearity≫. R_s analysis by Viaggi (2018a), and in the present work (Table 1 and Table 5), which suggests that the climate system is not self-sustaining because the variance of the response does not match the variance of the forcing, clearly corroborates this statement. Likely, orbitally-paced climate responses result by nonlinear changes trought the combined action of internal feedback mechanisms (albedo, GHGs, etc.) which are dependent on the mean climate state and are linked to the time scale of orbital cycles (Viaggi, 2018a), operating a strong energy transfer (Hays et al., 1976; Shackleton, 2000; Archer et al., 2004; Brovkin et al., 2007; Bintanja and van de Wal, 2008; Kohler et al., 2010; Lisiecki, 2010; Hansen et al., 2011; Shoenfelt et al., 2018; Viaggi, 2021a, 2021c). According to Ellis and Palmer (2016) glacial cycles are forced by orbital cycles and Milankovitch insolation, but regulated by ice-albedo and dust-albedo feedbacks. A corollary of Milankovitch’s theory is that the climatic effect of orbital forcings must be in phase and consequent to the cause because the climate system is inertial. As already argued by Hays et al. (1976) ≪the dominant 100,000-year climatic component has an average period close to, and is in phase with, orbital eccentricity≫, but this result has over time been marginalised mainly because of the low insolation fraction in the eccentricity band. Here it is the case to highlight how the new cross-spectrum analysis data between nominal forcings and the respective components of global δ¹⁸O, Red Sea RSL (Table 5) and EPICA stacks (Table 6, Figure 24) confirm the high cospectral coherence also of the 90-100 kyr cycle, and a general delayed pattern of proxies consistent with an inertial climate system (with the only exception of the EPICA short eccentricity). In particular, it should be noted that also the weak short eccentricity component of the MAI has a coherence of 0.72 (Table 6), and that the La2010 nominal short eccentricity signal acquired from Morlet wavelet filtering has a high cross-coherency with the associated components of global benthic δ¹⁸O (0.75) and the Red Sea RSL (0.78) during the last 500 kyr (Table 5). Assuming a certain margin of error due to different age models, the phase delay of EPICA obliquity and precession stacks respect to nominal forcings (4.2-kyr and ~4.0-kyr, respectively) is in agreement with the putative cause-effect relationship, and support the action of feedback mechanisms which are lagged by definition (Viaggi, 2018a, 2021a). On the other hand, the EPICA short eccentricity stack shows an advance of 4.3-4.5-kyr, which seems uncredible. Instead, the phase relationships of the short eccentricity within the EPICA record (i.e., with the same age model) are consistent with an inertial climate system affected by feedbacks since CO₂ and CH₄ lags δD (Viaggi, 2021a). The EPICA ‘anomalous’ result with the eccentricity nominal solution could be an artefact owing to the averaging of the lead-lag phase variability over a small number of short eccentricity peaks (16) during the last 800 kyr, perhaps in relation to possible inaccuracies in the age model and error propagation on long periodicity.

Nevertheless, these short eccentricity cospectral data cannot be undervalued simply because the energy involved is too low. Presumably there are other, different and complex mechanisms at stake, which we have probably underestimated or ignored. On the other hand, the weaker short eccentricity external forcing should not implies any problem if the ~100-kyr response is driven primarily by internal feedbacks which transfer most of the climate system energy (Hays et al., 1976; Rial, 1999; Shackleton, 2000; Liu et al., 2008; Lisiecki, 2010; Viaggi, 2018a, 2021a, 2021c). In particular, Shackleton (2000) believes that the ~100-kyr cycle is not generated from ice sheet dynamics, but it is probably the effect of orbital forcing eccentricity, with the changing atmospheric CO₂ concentration influencing the global carbon cycle and controlling the timing of glaciations (Raymo, 1997), and Lisiecki (2010) pointed towards an internally driven climate response phase-locked to eccentricity.

In this context, the ODH comes into play as it offers a possible physical mechanism that mitigates the obliquity forcing which is the main summer ice-melt agent at NH high latitudes as previously discussed. This mechanism, by favouring a longer life of the ice-sheet through the obliquity-cycle damping and ‘skipping’, combined with the Earth’s long-term cooling, may have increased the sensitivity of the polar climate system to the minima of MAI by triggering a strong feedback-related nonlinear energy transfer in the short eccentricity band.

4. Summary and conclusions

4.1. Main results

The evidences of obliquity damping based on proxy records can be summarised as follows:

The Antarctic orbital R_s demonstrates that since 560 kyr, a strong amplification of the short eccentricity signals has been occurring (up to 400%, 600%) coupled with damping of obliquity responses (up to ‒80%, –60%), confirming the marked asymmetry of the climate responses to orbital forcing. The PCA model of R_s data (EPICA, LR04 δ¹⁸O) suggests PC-1 to be a latent factor, indicating a post-MPT anticorrelation among obliquity and short eccentricity/precession R_s, which is related to the long-term growth of the cryosphere volume.
The PCA model integrating Plio-Pleistocene orbital R_s data and the long-term components of both LR04 δ¹⁸O and Site 846 SST records identifies two PCs that are strictly related to the δ¹⁸O/SST short eccentricity/precession (PC-1) and obliquity (PC-2) amplification. Both are linked to the long-term δ¹⁸O enrichment and SST reduction. PC-2 factor highlights the post-MPT anomalous depletion of the obliquity R_s. These factors corroborate the latent link among the increasing amplitude of all orbital climate responses, the obliquity damping and the Earth’s icy-state developed through four stages of step-wise growth (subtrend I to IV).
The spread of Plio-Pleistocene PCs factor exhibits two anticorrelation patterns of increasing absolute magnitude among forcing responses through time: positive spread showing high-obliquity associated with low-short eccentricity/precession R_s, and negative spread showing low-obliquity linked to high-short eccentricity/precession R_s. These response configurations are associated with three transition patterns of positive to negative spread including ONHG (Transition-1), INHG (Transition-2), and MBE (Transition-3), the latter being characterised by extremely high-magnitude and containing the MPT.
EPICA orbital SSA-stacks, which are rescaled on Plio-Pleistocene variance, exhibit a short eccentricity estimate of 12.2%, which is very high compared to the global δ¹⁸O value of 6.5%. In addition, an obliquity estimate of 4.5% is very low compared to the δ¹⁸O value of 9.9%. This is in agreement with the results of the R_s analysis, indicating a post-MPT short eccentricity amplification vs. obliquity damping.
Orbital SSA-components from the RSL record of the Red Sea exhibit a Plio-Pleistocene rescaled variance consistent with that of EPICA: short eccentricity (14.5%), obliquity (3.1%), and precession (2.8%). These data suggest post-MPT short eccentricity amplification vs. obliquity damping even in ESL fluctuations.
Antarctic SSA-structural signal observation of two or three low-amplitude 41-kyr obliquity peaks (glacial/interglacial) embedded in a weak ~93/75-kyr short eccentricity framework determined from dD, CO₂, and CH₄ component-2s, are very similar in shape to the global LR04 δ¹⁸O and Site 846 equatorial Pacific SST component-3-4s during the MPT. These shapes may be further evidence of an obliquity attenuation phenomenon linked to the short eccentricity, and seem observational reminiscent of the ‘obliquity-cycle skipping’ model.
Additional evidences of MPT and post-MPT anticorrelation between obliquity damping and short eccentricity amplification are highlighted from a variety of global and regional (Atlantic, Pacific, Mediterranean, Indian) climate-related records, hinting them to be a widespread feature of the Pleistocene climate system. However, this feature does not seem to be consistent with the nominal solutions.

The main results and considerations concerning the theoretical constrains required by ODH can be summarised as follows:

Studies on Greenland and Antarctica indicate a fast response of the cryosphere mass balance to recent atmospheric and SST changes. The EPICA stacks better approximates the global benthic δ¹⁸O, the latter probably represents the most lagged signal in the climate chain, and may be considered a proxy of the global atmospheric temperature averaged by GHGs. Thus, it is likely that the phase lags averaged among surface temperature proxies approximate the phase lag of the ice-volume better, overcoming the deep temperature lag bias of the benthic δ¹⁸O.
Assuming a fast response of the ice-volume to surface temperature changes (EPICA stack, SST), the mean value of obliquity lag <5.0 kyr (3.7±1.7 kyr) has been documented during the last 800 kyr. Also considering the benthic δ¹⁸O records, the obliquity mean phase lag is 5.3±0.6 kyr, which is very close to the theoretical threshold of 5.0 kyr and is significantly lower than the range of 6-10 kyr considered in previous studies.
Cross-coherency data demonstrate that the Red Sea RSL record approximates extremely well the glacio-eustatic sea-level fluctuations linked to ice-volume and paced by orbital forcings, also in the short eccentricity band. Global benthic δ¹⁸O lags Red Sea ESL by ~2.0 kyr in the obliquity band, suggesting a delay bias that could be attributed to benthic δ¹⁸O deep-water temperature signal. Thus, the benthic δ¹⁸O obliquity mean phase lag of 5.3 kyr could be corrected to 3.3 kyr, close to the EPICA stack/SST mean of 3.7 kyr, which is significantly lower than the theoretical threshold of 5-kyr.
Recent studies on variations in present-day satellite temporal gravity suggest the high sensitivity of the Earth to oblateness and highlight a concurrent slow-moderate negative (GIA rebound) and fast-robust positive (water-mass redistribution) J₂ changes during the recent postglacial and global warming context, resulting in a J₂ positive net change. Cross-spectral results from the Red Sea RSL over the last 500 kyr suggest a rapid and coherent oscillation in the obliquity band of the water layer component of the Earth’s oblateness.
It is hypothesised that the fast and robust J₂ water-mass redistribution component in the obliquity band, which is reasonably related to the glacio-eustatic variation of the sea level, could be a crucial element in determining +J2 ⟹ –during obliquity maxima (interglacial damping), and –J2 ⟹ +during obliquity minima (glacial damping). This would explain the fact that obliquity damping in climate proxies is basically symmetrical.

4.2. Obliquity damping hypothesis

The proposed hypothesis can be outlined as follows:

The widespread evidence from proxy records of anticorrelation between obliquity and short eccentricity/precession has been interpreted as an effect of the obliquity-oblateness feedback by critically reviewing its theoretical constrains, which support negative/positive secular change of obliquity for both low ice-volume phase lag (<5.0 kyr) and positive/negative net change of oblateness, respectively, in the obliquity band that are likely dominated by the J₂ water-mass redistribution component.
Obliquity damping during the interglacial stages might have contributed to the strengthening of the short eccentricity response by mitigating the obliquity’s ice-killing, favouring the obliquity-cycle skipping, and a ~100-kyr long-life feedback amplified ice-growth in the context of global icy-state.
Orbitals, including short eccentricity, may pace the frequency beat of the climate response. The phase-locked feedback mechanisms might have non-linearly transferred most of the system energy depending on the long-term climate state and the cycle duration, overcoming the energy excess ‘paradox’ especially with respect to the eccentricity band.
The observed transition patterns (TRA-1,-2,-3) of the orbitals R_s and the early onset of the short eccentricity response suggest the traditional notion of the MPT to be the final, high-magnitude, nonlinear transitional stage of a complex competing interaction between obliquity vs. short eccentricity forcing under the influence of both long-term cooling and obliquity-oblateness feedback, that had already started during the Piacenzian. The maximum expression of this mechanism would occur during TRA-3 that is associated with the maximum ice-volume development (subtrend IV) and a strong amplification of the obliquity response till the termination of the MPT (TH4).
The Plio-Pleistocene long-term cooling is a relevant background forcing in setting boundary conditions to orbital climate responses, and is characterised by four step-wise subtrends (I to IV), where a mild curvilinear shape is broken by slope changes representing four thresholds (TH1 to TH4) of mean climate state change.
The role of the long-term cooling is outlined as follows:

a.

The non-linearity of the orbital responses is increased by the scale effect of the ice-sheet growth on feedback mechanisms, which are sensitive to long orbital periods. Specifically, the amplitude increases in the obliquity band of the ice-volume changes. It is hypothesised that the glacio-eustatic oscillations in the obliquity band inducing oblateness variations by dominant/fast water mass redistribution component could have led to the overcoming of the thresholds sensitive to secular change of obliquity. This mechanism would explain why the 100-kyr cycle reached its maximum expression post-MPT, albeit after a period of early/weak manifestations as early as in Piacenzian.

b.

The icy-state may have increased the sensitivity of the polar climate system to the minima of MAI by triggering a strong nonlinear energy transfer in the short eccentricity band by feedback mechanisms (primarily ice-albedo and GHGs).

The coupling between obliquity damping and short eccentricity amplification through the long-term cooling could be a missing link in the Pleistocene climate system and MPT debate. ODH could explain the rise of a dominant ~100-kyr response during the Mid-Late Pleistocene as its enhanced sensitivity to short eccentricity by a competing interaction with obliquity through an obliquity-oblateness feedback mechanism during the development of the Earth’s icy-state. This could suggest a different impact of the climate friction than previously thought.

4.3. Outlook

ODH needs independent corroboration studies. Here we make some recommendations for future research:

Unbiased estimates of the ice-volume lag in the obliquity band.
Modelling the changes in Earth’s oblateness in the obliquity band and its effects on the forcing, determined especially by considering the dominant/fast water mass redistribution component. This could make the difficult and uncertain solid Earth’s oblateness component less stringent among the theoretical constrains.
Modelling the hypothesis with full observational evidence control.
Role of the MAI and feedback mechanisms in a context of the Earth’s long-term icy-state.

Funding

This research did not receive any funding from funding agencies of the public, commercial, or not-for-profit sectors.

Data Availability Statement

Supplementary data to this article can be found online at link

Acknowledgements

The author is grateful to the anonymous Colleague for the useful suggestions and enhancements to the paper.

Conflicts of Interest

The author declares that there are no competing interests.

Abbreviations

CCF: cross-correlation function; CSI: cubic spline interpolation; comp: component; EPICA: European Project for Ice Coring in Antarctica; ESL: global (eustatic) sea level; FFS: Fourier frequency spectrum; GHG: green-house gas; GIA: glacial isostatic adjustment; GRACE: gravity recovery and climate experiment; INHG: Intensification of Northern Hemisphere Glaciation; IRD: ice-rafted debris; J2: Earth’s dynamic oblateness; MAI: mean annual insolation; MBE: Mid- Brunhes Event; MBO: Mid-Brunhes Oscillation; MIS: marine isotope stage; MPT: Mid-Pleistocene transition; ODH: obliquity damping hypothesis; ONHG: Onset of the Northern Hemisphere Glaciation; PCA: principal component analysis; Rs: response sensitivity to orbital forcings; RSL: relative sea level; SSA: singular spectrum analysis; SST: sea surface temperature; subcomp: subcomponent; W: gravitational redistribution of water in the oceans and atmosphere

Note

Reconstruction of orbital SSA-components as follows (sub = subcomponent):

δ¹⁸O short eccentricity = (detrend comp-1 sub-3-11 * 0.95 + comp-2 sub-1-4;9-15 * 5.5) / 6.45;

δ¹⁸O obliquity = (comp-3-4 * 8.1 + comp-2 sub-5-8 * 1.8) / 9.9;

δ¹⁸O precession = comp-5-7;

δ¹⁸O lomg-term trend = comp-1 sub-1;

SST short eccentricity = (comp-2 sub-1-4 * 3.65 + comp-2 sub-5-8;11-15 * 0.5 + detr comp-1 sub-3-11 * 0.93) / 5.08; SST obliquity = (comp-3-4 * 4.5 + comp-2 sub-9-10 * 0.27 + comp-2 sub-5-8;11-15 * 0.76) / 5.53;

SST precession = comp-5-7;

SST lomg-term trend = comp-1 sub-1;

References

Abe-Ouchi A, Saito F, Kawamura K, Raymo ME, Okuno J, Takahashi K, Blatter H (2013) Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet volume. Nature 500. [CrossRef]
Ao H, Dekkers MJ, Qin L, Xiao G (2011) Age models of ODP Site 184-1143. PANGAEA. [CrossRef]
Archer D, Martin P, Buffett B, Brovkin V, Rahmstorf S, Ganopolski A (2004) The importance of ocean temperature to global biogeochemistry. Earth Planet Sci. Lett. 222(2004):333–348. [CrossRef]
Bajo P, Drysdale RN, Woodhead JD, Hellstrom JC, Hodell D, Ferretti P, Voelker AHL, Zanchetta G, Rodrigues T, Wolff E, Tyler J, Frisia S, Spötl C, Fallick AE (2020) Persistent influence of obliquity on ice age terminations since the Middle Pleistocene transition. Science 367, 1235–1239 (2020). [CrossRef]
Barth AM, Clark PU, Bill NS, He F, Pisias NG (2018) Climate evolution across the Mid-Brunhes Transition. Clim. Past, 14, 2071–2087, 2018. [CrossRef]
Bartoli G, Hönisch B, Zeebe RE (2011) Atmospheric CO2 decline during the Pliocene intensification of northern hemisphere glaciations. Paleoceanography 26, PA4213. [CrossRef]
Bazin L, Landais A, Lemieux-Dudon B, Toyé Mahamadou Kele H, Veres D, Parrenin F, Martinerie P, Ritz C, Capron E, Lipenkov V, Loutre M-F, Raynaud D, Vinther B, Svensson A, Rasmussen SO, Severi M, Blunier T, Leuenberger M, Fischer H, Masson-Delmotte V, Chappellaz J, Wolff E (2013) An optimized multi-proxy, multi-site Antarctic ice and gas orbital chronology (AICC2012): 120–800 ka. Clim. Past, 9, 1715–1731, 2013. [CrossRef]
Bereiter B, Eggleston S, Schmitt J, Nehrbass-Ahles C, Stocker TF, Fischer H, Kipfstuhl S, Chappellaz J (2015) Revision of the EPICA Dome C CO2 record from 800 to 600 kyr before present. Geophys. Res. Lett., 42, 1-8. [CrossRef]
Berends CJ, Köhler P, Lourens LJ, van de Wal RSW (2021) On the cause of the mid-Pleistocene transition. Reviews of Geophysics, 59, 1-20. [CrossRef]
Berger A, Loutre MF, Tricot C (1993) Insolation and Earth’s orbital periods. Journal of Geophysical Research Atmospheres, 98, D6, 93. 19 June. [CrossRef]
Berger A, Loutre MF (1997) Long-term variations in insolation and their effects on climate, the LLN experiments. Surv. Geophys 18, 147–161. [CrossRef]
Berger A, Li XS, Loutre MF (1999) Modeling northern hemisphere ice volume over the last 3 Ma. Quaternary Science Reviews, 18, 1–11.
Berger A, Yin Q (2012) Modeling the interglacials of the last 1 million years. In: Berger et al. (eds) Climate change: inferences from paleoclimate and regional aspects. Springer-Verlag, Wien 2012.
Berger A, Mélice JL, Loutre MF (2005) On the origin of the 100-kyr cycles in the astronomical forcing, Paleoceanogr., 20:PA4019. [CrossRef]
Berger A, Yin QZ, Herold N (2012) MIS-11 and MIS-19, analogs of our Holocene interglacial. In: 3rd Int. Conf. on Earth Syst. Model., Hamburg. https://meetingorganizer.copernicus.org/3ICESM/3ICESM-11.pdf.
Bills BG (1994) Obliquity-oblateness feedback: Are climatically sensitive values of obliquity dynamically unstable? Geophys. Res. Lett., 21, 3, 177-180, , 1994. 1 February.
Bills BG (1998) An oblique view of climate. Nature, 396, 405-406, 3 December.
Bintanja R, van de Wal RSW (2008) North American ice-sheet dynamics and the onset of 100,000-year glacial cycles. Nature, Vol 454, , 869- 872. 14 August. [CrossRef]
Box JE, Hubbard A, Bahr DB, Colgan WT, Fettweis X, Mankoff KD, Wehrlé A, Noël B, van den Broeke MR, Wouters B, Bjørk AA, Fausto RS (2022) Greenland ice sheet climate disequilibrium and committed sea-level rise. Nature Climate Change, Vol 12, 22, pp. 808–813. 20 September. [CrossRef]
Brovkin V, Ganopolski A, Archer D, Rahmstorf S (2007) Lowering of glacial atmospheric CO2 in response to changes in oceanic circulation and marine biogeochemistry. Paleoceanogr 22, PA4202. [CrossRef]
Chalk TB, Hain MP, Foster GL, Rohling EJ, Sexton PF, Badger MPS, Cherry SG, Hasenfratz AP, Haug GH, Jaccard SL, Martínez-García A, Pälike H, Pancost RD, Wilson PA (2017) Causes of ice age intensification across the Mid-Pleistocene Transition. Proceedings of the National Academy of Sciences, 17, 1-6. 20 November. [CrossRef]
Chao BF (2006) Earth’s oblateness and its temporal variations. C. R. Geoscience, 338 (2006), 1123–1129. [CrossRef]
Cheng M, Tapley BD (2004) Variations in the Earth’s oblateness during the past 28 years. J. of Geophys. Res.,109, B09402. [CrossRef]
Cheng M, Tapley BD, Ries JC (2013) Deceleration in the Earth’s oblateness. J. of Geophys. Res.: Solid Earth, 118, 740–747. [CrossRef]
Cheng X, Tian J, Wang P (2004) Stable foraminifer isotope composition of ODP Site 184-1143, South China Sea. PANGAEA. [CrossRef]
Clark PU, Archer D, Pollard D, Blum JD, Rial JA, Brovkin V, Mix AC, Pisias NG, Roy M (2006) The middle Pleistocene transition: characteristics, mechanisms, and implications for long-term changes in atmospheric pCO2. Quat. Sci. Rev. [CrossRef]
Crucifix M (2012) - Oscillators and relaxation phenomena in Pleistocene climate theory. Phil. Trans. R. Soc. A (2012) 370, 1140–1165. [CrossRef]
Cuffey KM, Clow GD (1997) - Temperature, accumulation, and ice sheet elevation in central Greenland through the last deglacial transition. Journal of Geophysical Research, vol. 102, No. C12, pages 26,383-26,396, november 30, 1997.
Daradich A, Huybers P, Mitrovica JX, Chan N-H, Austermann J (2017) The influence of true polar wander on glacial inception in North America. Earth and Planetary Science Letters, 461 (2017), 96–104. [CrossRef]
Ditlevsen PD, Ashwin P (2018) Complex Climate Response to Astronomical Forcing: The Middle-Pleistocene Transition in Glacial Cycles and Changes in Frequency Locking. Front. Phys 6, 62. [CrossRef]
Elderfield H, Ferretti P, Greaves M, Crowhurst S, McCave IN, Hodell D, Piotrowski AM (2012) Evolution of Ocean Temperature and Ice Volume Through the Mid-Pleistocene Climate Transition. Science 337, 704–709 (2012). [CrossRef]
Ellis R, Palmer M (2016) Modulation of ice ages via precession and dust-albedo feedbacks. Geoscience Frontiers, v. 7 (2016), pp. 891-909. [CrossRef]
Elsner JB, Tsonis AA (1996) Singular spectrum analysis: a new tool in time series analysis. Springer, 1996.
Ettema J, van den Broeke MR, van Meijgaard E, van de Berg WJ, Bamber JL, Box JE, Bales RC (2009), Higher surface mass balance of the Greenland ice sheet revealed by highresolution climate modeling, Geophys. Res. Lett., 36, L12501. [CrossRef]
Farmer JR, Hönisch B, Haynes LL, Kroon D, Jung S, Ford HL, Raymo ME, Jaume-Seguí M, Bell DB, Goldstein SL, Pena LD, Yehudai M, Kim J (2019) Deep Atlantic Ocean carbon storage and the rise of 100,000-year glacial cycles. Nature Geosci. [CrossRef]
Feng F, Bailer-Jones CAL (2015) Obliquity and precession as pacemakers of Pleistocene deglaciations. Quaternary Science Reviews, v. 122, , pp. 166-179. 15 August. [CrossRef]
Ghil M, Allen RM, Dettinger MD, Ide K, Kondrashov D, Mann ME, Robertson A, Saunders A, Tian Y, Varadi F, Yiou P (2002) Advanced spectral methods for climatic time series. Rev. Geophys. 40(1):3.1–3.41. [CrossRef]
Govin A, Braconnot P, Capron E, Cortijo E, Duplessy J-C, Jansen E, Labeyrie L, Landais A, Marti O, Michel E, Mosquet E, Risebrobakken B, Swingedouw D, Waelbroeck C (2012) Persistent influence of ice sheet melting on high northern latitude climate during the early Last Interglacial. Clim. Past, 8, 483–507, 2012. [CrossRef]
Grant KM, Rohling EJ, Bronk Ramsey C, Cheng H, Edwards RL, Florindo F, Heslop D, Marra F, Roberts AP, Tamisiea ME, Williams F (2014) Sea-level variability over five glacial cycles. Nat. Commun., 5, 5076. [CrossRef]
Hanna E, Huybrechts P, Konrad Steffen K, Cappelen J, Huff R, Shuman C, Irvine-Fynn T, Wise S, Griffiths M (2008) Increased Runoff from Melt from the Greenland Ice Sheet: A Response to Global Warming. American Meteorological Society, 21, 331-341. [CrossRef]
Hansen J, Sato M, Kharecha P, Russell G, Lea DW, Siddall M (2007) Climate change and trace gases. Phil Trans. R. Soc. A 365, 1925–1954. [CrossRef]
Hansen J, Sato M, Kharecha P, von Schuckmann K (2011) Earth’s energy imbalance and implications. Atmos. Chem. Phys 11, 13421–13449. [CrossRef]
Hasenfratz AP, Jaccard SL, Martínez-García A, Sigman DM, Hodell DA, Vance D, Bernasconi SM, Kleiven HF, Haumann FA, Haug GH (2019) The residence time of Southern Ocean surface waters and the 100,000-year ice age cycle. Sci. 363, 1080-1084 (2019), 1-5. [CrossRef]
Hassani H (2007) Singular spectrum analysis: methodology and comparison. J. Data Sci. 5(2007):239–257. 2: 5(2007).
Hays JD, Imbrie J, Shackleton NJ (1976) Variations in the Earth’s Orbit: Pacemaker of the Ice Ages. Science, Vol. 194, No. 4270 (Dec. 10, 1976), pp. 1121-1132. https://www.jstor.org/stable/1743620.
Head MJ, Gibbard PL (2005) Early–Middle Pleistocene transitions. An overview and recommendation for the defining boundary. From: Head MJ & Gibbard PL (eds). Early–middle Pleistocene transitions: the land–ocean evidence. Geological Society, London, Special Publications, 247, 1–18. 0305–8719/05/$15© The Geological Society of London 2005.
Head MJ, Pillans B, Farquhar SA (2008) The Early–Middle Pleistocene transition: characterization and proposed guide for the defining boundary. Episodes 31(2):255–259. 2: Episodes 31(2).
Herbert TD, Peterson LC, Lawrence KT, Liu Z (2010) Tropical ocean temperatures over the past 3.5 million years. Sci. 2010, 328, 1530–1534. [CrossRef] [PubMed]
Hilgen FJ, Lourens LJ, Berger A, Loutre MF (1993) Evaluation of the astronomically calibrated time scale for the Late Piocene and the earliest Pleistocene. Paleoceanography 8, 549–565. [CrossRef]
Hodell DA, Channell JET, Curtis JH, Romero OE, Rohl U (2008) Onset of ‘‘Hudson Strait’’ Heinrich events in the eastern North Atlantic at the end of the middle Pleistocene transition (~640 ka)? Paleoceanogr., 23, PA4218. [CrossRef]
Hodell D, Lourens L, Crowhurst S, Konijnendijk T, Tjallingii R, Jiménez-Espejo F, Skinner L, Tzedakis PC, Shackleton Site Project Members (2015) A reference time scale for Site U1385 (Shackleton Site) on the SW Iberian Margin. Global and Planetary Change, 133 (2015), 49–64. [CrossRef]
Hönisch B, Hemming NG, Archer D, Siddall M, McManus JF (2009) Atmospheric carbon dioxide concentration across the mid-Pleistocene transition. Science 324, 1551. [CrossRef]
Huybers P (2006) Early Pleistocene glacial cycles and the integrated summer insolation forcing. Science, 313. [CrossRef]
Huybers P (2007) Glacial variability over the last two million years: an extended depth-derived age model, continuous obliquity pacing, and the Pleistocene progression. Quat. Sci. Rev. 26 (2007):37–55. [CrossRef]
Huybers P, Tziperman E (2008) Integrated summer insolation forcing and 40,000-year glacial cycles: The perspective from an ice-sheet/energy-balance model. Paleoceanography, vol. 23, 1-18. [CrossRef]
Imbrie J, Berger A, Boyle EA, Clemens SC, Duffy A, Howard WR, Kukla G, Kutzbach J, Martinson DG, Mcintyre A, Mix AC, Molfino B, Morley, JJ, Peterson LC, Pisias NG, Prell WL, Raymo ME, Shackleton NJ, Toggweile JR (1993) On the structure and origin of major glaciation cycles 2. The 100,000-year cycle. Paleoceanogr. 8(6):699–735.
Imbrie J, Imbrie-Moore A, Lisiecki LE (2011) A phase-space model for Pleistocene ice volume. Earth Planet Sci Lett 307, 94–102. [CrossRef]
Jouzel J et al. (1997) Validity of the temperature reconstruction from water isotopes in ice cores, J. Geophys. Res., 102(C12), 26,471–26,487.
Jouzel J et al. (2007) Orbital and millennial Antarctic climate variability over the past 800,000 years. Sci 317, 793–796. [CrossRef]
Kender S, Ravelo AC, Worne S, Swann GEA, Leng MJ, Asahi H, Becker J, Detlef H, Aiello IW, Andreasen D, Hall IR (2018) Closure of the Bering Strait caused Mid-Pleistocene Transition cooling. Nature Commun. (2018) 9, 5386. [CrossRef]
King MD, Howat IM, Candela SG, Noh MJ, Jeong S, Noël BPY, van den Broeke MR, Wouters B, Negrete A (2020) Dynamic ice loss from the Greenland Ice Sheet driven by sustained glacier retreat. Communications Earth & Environment (2020) 1, 1. [CrossRef]
Köhler P, Bintanja R, Fischer H, Joos F, Knutti R, Lohmann G, Masson-Delmotte V (2010) What caused Earth’s temperature variations during the last 800,000 years? Data-based evidence on radiative forcing and constraints on climate sensitivity. Quat. Sci. Reviews 29 (2010) 129–145. [CrossRef]
Köhler P, van de Wal RSW (2020) Interglacials of the Quaternary defined by northern hemispheric land ice distribution outside of Greenland. Nature Commun. (2020) 11, 5124. [CrossRef]
Konijnendijk TYM, Ziegler M, Lourens LJ (2015) On the timing and forcing mechanisms of late Pleistocene glacial terminations: Insights from a new high- resolution benthic stable oxygen isotope record of the eastern Mediterranean. Quaternary Science Reviews, 129 (2015), 308-320. [CrossRef]
Lang N, Wolff EW (2011) Interglacial and glacial variability from the last 800 ka in marine, ice and terrestrial archives. Clim. Past, 7, 361–380, 2011. [CrossRef]
Laskar J, Joutel F, Boudin F (1993) - Orbital, precessional and insolation quantities for the Earth from -20 Myr to +10 Myr. Astronomy & Astrophysics, Vol. 270, pp. 522– 533.
Laskar J, Robutel P, Joutel F, Gastineau M, Correia ACM, Levrard B (2004) A long term numerical solution for the insolation quantities of the Earth. Astron. Astrophys 428, La2004.
Laskar J, Fienga A, Gastineau M, Manche H (2011) La2010: a new orbital solution for the long-term motion of the Earth. Astron. Astrophys 532, A89. [CrossRef]
Lawrence DM, Slater AG, Tomas RA, Holland MM, Deser C (2008) Accelerated Arctic land warming and permafrost degradation during rapid sea ice loss. Geophys Res Lett 35, L11506. [CrossRef]
Lawrence KT, Herbert TD, Brown CM, Raymo ME, Haywood AM (2009) High-amplitude variations in North Atlantic sea surface temperature during the early Pliocene warm period. Paleoceanography, vol. 24, PA2218. [CrossRef]
Lawrence KT, Sosdian S, White HE, Rosenthal Y (2010) North Atlantic climate evolution through the Plio-Pleistocene climate transitions. Earth and Planetary Science Letters, 300 (2010), 329–342. [CrossRef]
Lemke P, Ren J, Alley RB, Allison I, Carrasco J, Flato G, Fujii Y, Kaser G, Mote P, Thomas RH, Zhang T (2007) Observations: Changes in Snow, Ice and Frozen Ground. In: Climate Change 2007.
Levrard B, Laskar J (2003) Climate friction and the Earth’s obliquity. Geophys. J. Int 154, 970–990. [CrossRef]
Li L, Qianyu L, Jun T, Pinxian W, Hui W, Zhonghui L. (2011) A 4-Ma record of thermal evolution in the tropical western Pacific and its implications on climate change. Earth and Planetary Science Letters, 309 (2011), 10–20. [CrossRef]
Lisiecki LE, Raymo ME (2005) A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records. Paleoceanogr 20, PA1003. [CrossRef]
Lisiecki LE, Raymo ME (2007) Plio–Pleistocene climate evolution; trends and transitions in glacial cycle dynamics. Quat. Sci. Rev. 26(2007):56– 69. [CrossRef]
Lisiecki LE, Raymo ME (2009) Diachronous benthic δ18O responses during late Pleistocene terminations, Paleoceanography, 24, PA3210. d.
Lisiecki LE (2010) Links between eccentricity forcing and the 100,000-year glacial cycle. Nature Geoscience, vol 3, may 2010, 349-352. [CrossRef]
Lisiecki LE (2014) Atlantic overturning responses to obliquity and precession over the last 3 Ma. Paleoceanogr., 29, 71–86. [CrossRef]
Liu Z, Cleaveland LC, Herbert TD (2008) Early onset and origin of 100-kyr cycles in Pleistocene tropical SST records. Earth and Planetary Science Letters 265 (2008), 703–715. [CrossRef]
Loulergue L, Schilt A, Spahni R, Masson-Delmotte V, Blunier T, Lemieux B, Barnola JM, Raynaud D, Stocker TF, Chappellaz J (2008) Orbital and millennial-scale features of atmospheric CH4 over the past 800,000 years. Nature 453. [CrossRef]
Luthi D, Le Floch M, Bereiter B, Blunier T, Barnola JM, Siegenthaler U, Raynaud D, Jouzel J, Fischer H, Kawamura K, Stocker TF (2008) High-resolution carbon dioxide concentration record 650,000–800,000 years before present. Nature 453. [CrossRef]
Martinez-Boti MA, Foster GL, Chalk TB, Rohling EJ, Sexton PF, Lunt DJ, Pancost RD, Badger MPS, Schmidt DN (2015) Plio-Pleistocene climate sensitivity evaluated using high-resolution CO2 records. Nature 518. [CrossRef]
Masson-Delmotte V et al. (2006) Past and future polar amplification of climate change: climate model intercomparisons and ice-core constraints. Climate Dynamics (2006) 26, 513–529. [CrossRef]
Masson-Delmotte V, Braconnot P, Hoffmann G, Jouzel J, Kageyama M, Landais A, Lejeune Q, Risi C, Sime L, Sjolte J, Swingedouw D, Vinther B (2011) Sensitivity of interglacial Greenland temperature and δ18O: ice core data, orbital and increased CO2 climate simulations. Clim. Past, 7, 1041–1059, 2011. [CrossRef]
McClymont EL, Sosdian SM, Rosell-Melé A, Rosenthal Y (2013) Pleistocene sea-surface temperature evolution: Early cooling, delayed glacial intensification, and implications for the mid-Pleistocene climate transition. Earth-Science Reviews, Volume 123, 173-193. [CrossRef]
Meyers SR and Hinnov LA (2010), Northern Hemisphere glaciation and the evolution of Plio-Pleistocene climate noise, Paleoceanography, 25, PA3207. d.
Miller GH, Alley RB, Brigham-Grette J, Fitzpatrick JJ, Polyak L, Serreze MC, White JWC (2010) Arctic amplification: can the past constrain the future? Quaternary Science Reviews, 29 (2010), 1779–1790. [CrossRef]
Mitrovica JX, Peltier WR (1993) - Present-day secular variations in the zonal harmonics of the Earth’s geopotential. J. Geophys. Res., 98 (B3), 4509–4526. [CrossRef]
Mukhin D, Gavrilov A, Loskutov E, Kurths J, Feigin A (2019) Bayesian Data Analysis for Revealing Causes of the Middle Pleistocene Transition. Scientific Reports (2019), 9, 7328. [CrossRef]
Nakada M, Okuno J (2017) Secular variations in zonal harmonics of Earth’s geopotential and their implications for mantle viscosity and Antarctic melting history due to the last deglaciation. Geophys. J. Int. (2017) 209, 1660–1676. [CrossRef]
Nerem RS, Wahr J (2011) Recent changes in the Earth’s oblateness driven by Greenland and Antarctic ice mass loss. Geophys. Res. Lett. 38, L13501. [CrossRef]
Nyman KHM, Ditlevsen PD (2019) The middle Pleistocene transition by frequency locking and slow ramping of internal period. Climate Dynamics. [CrossRef]
Past interglacials working group of PAGES (2016), Interglacials of the last 800,000 years. Rev. Geophys 54, 162–219. [CrossRef]
Pattyn F, Ritz C, Hanna E, Asay-Davis X, DeConto R, Durand G, Favier L, Fettweis X, Goelzer H, Golledge NR, Munneke PK, Lenaerts JTM, Nowicki S, Payne AJ, Robinson A, Seroussi H, Trusel LD, van den Broeke M (2018) The Greenland and Antarctic ice sheets under 1. 5 °C global warming. Nature Climate Change. [CrossRef]
Quinn C, Sieber J, von der Heydt AS, Lenton TM (2018) The Mid-Pleistocene Transition induced by delayed feedback and bistability. Dynamics and Statistics of the Climate System, 2018, 1–17. [CrossRef]
Raymo ME (1997) The timing of major climate terminations. Paleoceanography, Vol. 12, No. 4, Pages 577-585, August 1997.
Raymo ME, Nisancioglu K (2003) The 41 kyr world: Milankovitch’s other unsolved mystery. Paleoceanography, vol. 18, no. 1. [CrossRef]
Rial JA (1999) Pacemaking the Ice Ages by frequency modulation of Earth’s orbital eccentricity. Science 285, 564–568.
Rial JA, Pielke RA Sr, Beniston M, Claussen M, Canadell J, Cox P, Held H, deNoblet-Ducoudré N, Prinn R, Reynolds JF, Salas JD (2004) Nonlinearities, feedbacks and critical thresholds within the Earth’s climate system. Clim. Chang. 65(1–2):11–38, 2004.
Rubincam DP (1993) The obliquity of Mars and climate friction. J. Geophys. Res 98, 10827–10832.
Ruddiman WF (2003) Orbital insolation, ice volume, and greenhouse gases. Quat. Sci. Rev. 22, 1597–1629.
Russon T, Elliot M, Sadekov A, Cabioch G, Corrège T, De Deckker P (2011), The mid-Pleistocene transition in the subtropical southwest Pacific, Paleoceanography, 26, PA1211.
Seki O, Foster GL, Schmidt DN, Mackensen A, Kawamura K, Pancost RD (2010) Alkenone and boron-based Pliocene pCO2 records. Earth Planet Sci Lett 292(2010):201–211. [CrossRef]
Seo K-W, Chen J, Wilson CR, Lee C-K (2015) Decadal and quadratic variations of Earth’s oblateness and polar ice mass balance from 1979 to 2010. Geophys. J. Int. (2015) 203, 475–481. [CrossRef]
Shackleton NJ (1967) Oxygen isotope analysis and pleistocene temperature reassessed. Nature 215, 15–17.
Shackleton NJ (2000) The 100,000-year ice-age cycle identified and found to lag temperature, carbon dioxide, and orbital eccentricity. Sci 289, 1897–1902.
Shoenfelt EM, Winckler G, Lamy F, Anderson RF, Bostick BC (2018) Highly bioavailable dust-borne iron delivered to the Southern Ocean during glacial periods. PNAS, 1-6,. [CrossRef]
Skinner LC, Shackleton NJ (2005), An Atlantic lead over Pacific deep-water change across Termination I: Implications for the application of the marine isotope stage stratigraphy, Quat. Sci. Rev., 24, 571–580. [CrossRef]
Steinberger B, Spakman W, Japsen P, Torsvik TH (2015) The key role of global solid-Earth processes in preconditioning Greenland’s glaciation since the Pliocene. Terra Nova, 27, 1–8, 2015. [CrossRef]
Steinhilber, F., Abreu, J.A., Beer, J., Brunner, I., Christl, M., Fischer, H., Heikkil¨a, U., Kubik, P.W., Mann, M., McCracken, K.G., Miller, H., Miyahara, H., Oerter, H., Wilhelms, F., 2012. 9,400 years of cosmic radiation and solar activity from ice cores and tree rings. PNAS, April 109 (16), 5967–5971. 10 April. [CrossRef]
Sun Y, Riva R, Ditmar P, Rietbroek R (2019) Using GRACE to explain variations in the Earth’s oblateness. Geophysical Research Letters, 46. [CrossRef]
Sun Y, Riva R (2020) A global semi-empirical glacial isostatic adjustment (GIA) model based on Gravity Recovery and Climate Experiment (GRACE) data. Earth Syst. Dynam., 11, 129–137, 2020. [CrossRef]
Tarnocai C, Canadell JG, Schuur EAG, Kuhry P, Mazhitova G, Zimov S (2009) Soil organic carbon pools in the northern circumpolar permafrost region. Global Biogeochem. Cycles 23, GB2023. [CrossRef]
Vautard R, Ghil M (1989) Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series. Physica D 35(1989):395–424 North-Holland, Amsterdam.
Viaggi P (2018a) δ18O and SST signal decomposition and dynamic of the Pliocene-Pleistocene climate system: new insights on orbital nonlinear behavior vs. long-term trend. Prog. in Earth and Planet. Sci., 2018, 5:81, 7 December 2018. [CrossRef]
Viaggi P (2018b) δ18O and SST signal decomposition and dynamic of the Pliocene-Pleistocene climate system: new insights on orbital nonlinear behavior vs. long-term trend. Prog. in Earth and Planet. Sci., 2018, 5:81, 7 December 2018, Supplementary data set: https://static-content.springer.com/esm/art%3A10.1186%2Fs40645-018-0236-z/MediaObjects/40645_2018_236_MOESM2_ESM.xls. 7 December.
Viaggi P, Scotti P, Previde Massara E, Knezaurek G, Menichetti E, Piva A, Torricelli S, Gambacorta G (2019) Paleoceanographic Evolution of mid-Cretaceous Paleobioproductivity and Paleoredox Chemometric Signals. Link to Cretaceous Long-term Eustatic Climax from a Central Atlantic Well Sequence. Conference presentation, STRATI 2019, 3rd International Congress on Stratigraphy (2-, Milan). 5 July. [CrossRef]
Viaggi P (2021a) Quantitative impact of astronomical and sun-related cycles on the Pleistocene climate system from Antarctica records. Quaternary Science Advances, Volume 4, October 2021. [CrossRef]
Viaggi P (2021b) Quantitative impact of astronomical and sun-related cycles on the Pleistocene climate system from Antarctica records. Quaternary Science Advances, Volume 4, October 2021, Supplementary data set: https://ars.els-cdn.com/content/image/1-s2.0-S2666033421000162-mmc2.xlsx.
Viaggi P (2021c) Evidence of Climate Mitigation Feedbacks by Global Forest Cover Expansion during the Pliocene. Int J Environ Sci Nat Res, 2021, 28(5), 1-5. [CrossRef]
Waelbroeck C, Levi C, Duplessy J-C, Labeyrie L, Michel E, Cortijo E, Bassinot F, Guichard F (2006) Distant origin of circulation changes in the Indian Ocean during the last deglaciation. Earth Planet. Sci. Lett., 243, 244–251.
Whitehouse PL (2018) Glacial isostatic adjustment modelling: historical perspectives, recent advances, and future directions. Earth Surf. Dynam., 6, 401–429, 2018. [CrossRef]
Williams DM, Kasting JF, Frakes LA (1998) Low-latitude glaciations and rapid changes in the Earth’s obliquity explained by obliquity–oblateness feedback. Nature 396, 453–455. [CrossRef]
Wunsch C (2004) Quantitative estimate of the Milankovitch-forced contribution to observed Quaternary climate change. Quat. Sci. Rev. 23(2004):1001–1012. [CrossRef]
Yehudai M, Kim J, Penac LD, Jaume-Seguì M, Knudson KP, Bolge L, Malinverno A, Bickert T, Goldstein SL (2021) Evidence for a Northern Hemispheric trigger of the 100,000-y glacial cyclicity. PNAS 2021, Vol. 118 No. 46, pp. 1-11. [CrossRef]
Yin J, Overpeck JT, Griffies SM, Hu A, Russell JL, Stouffer RJ (2011) Different magnitudes of projected subsurface ocean warming around Greenland and Antarctica. Nature Geoscience, 1-5. [CrossRef]
Yin QZ, Berger A (2012) Individual contribution of insolation and CO2 to the interglacial climates of the past 800,000 years. Clim Dyn (2012) 38, 709–724. [CrossRef]
Zachos JC, Pagani M, Sloan L, Thomas E, Billups K (2001) Trends, rhythms, and aberrations in global climate 65 Myr to present. Science 292, pp.686-693. Zachos JC, Dickens GR, Zeebe RE (2008) An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics. NATURE, Vol. 451, 17 January 2008. [CrossRef]
Zharkova V, Shepherd SJ, Zharkov SI (2012) Principal component analysis of background and sunspot magnetic field variations during solar cycles 21–23. Mon. Not. R. Astron. Soc. 424, 2943–2953 (2012). [CrossRef]

Figure 1. Analysis of the orbital response sensitivity during the last 800 kyr (fuchsia: short eccentricity; brown: obliquity; green: precession) according to temporal segments binned at 80 kyr from EPICA: (a) δD; (b) CH₄; (c) CO₂; and (d) LR04 δ¹⁸O records. The balance line (0%) indicates that the variance of the response matches the variance in forcing, indicating a self-sustained climate system that is paced only by orbital forcings. Values significantly larger than zero indicate a nonlinear reinforcement of the response signal. Negative R_s indicates a response variance lower than the orbital forcing, suggesting a damping mechanism for the climate system. Note the impressive asymmetry of the climate responses with extremely robust amplitude magnifications and moderate damping. Scale of the y-axis up to -200% for graphical requirements. MBE = Mid-Brunhes Event; MPT = Mid-Pleistocene transition. SSA-components are acquired from EPICA (Viaggi, 2021b), LR04 (Viaggi, 2018b). Nominal solutions are from Laskar et al. (2004) and Laskar et al. (2011).

Figure 2. Standardised PCA model of EPICA δD, CO₂, CH₄ and LR04 δ¹⁸O orbital R_s during the last 800 kyr, binned at 80-kyr: (a) PC-1 loadings; (b) PC-1 regression factor score (blue line), and δ¹⁸O exponential trend mean (green line). The asterisk indicates the transition to MBE. MBE = Mid-Brunhes Event; MPT = Mid-Pleistocene transition.

Figure 3. Morlet wavelet power spectra of EPICA SSA-stacks: (a) short eccentricity; (b) obliquity; (c) precession. MBE = Mid-Brunhes Event; MPT = Mid-Pleistocene transition. Data from Viaggi (2021b).

Figure 4. (a) EPICA dominant short eccentricity stack (51.6% variance) and related cycles defined on maxima (black horizontal lines) and numbered from 1st to 8th (left), compared with δD, CO₂, CH₄ SSA-component-2 time-series (obliquity signal embedded in a weak short eccentricity). Asterisks on δD comp-2 (not shown on the other comp-2s) exhibit two or three low-amplitude 41-kyr obliquity peaks (glacial/interglacial) embedded in a weak ~93/75-kyr framework after the MPT, a possible observational reminiscent of the ‘obliquity-cycle skipping’ (Huybers, 2007); (b) FFS of the δD comp-2, very similar to the CO₂ and CH₄ signals (not shown), showing a dominant 41-kyr cycle embedded in a weak ~93/75-kyr framework. Coloured lines are first-order autoregressive AR(1) curves (50%, 90%, 95%, 99%, and 99.9% significance levels). MBE = Mid-Brunhes Event; MPT = Mid-Pleistocene transition. Data from Viaggi (2021b).

Figure 5. Sensitivity of Pliocene-Pleistocene orbital response (blue: short eccentricity, brown: obliquity, green: precession) of (a, b) global LR04 δ¹⁸O and (c, d) equatorial ODP Site 846 SST as a function of the long-term mean climate state (left x-axis: exponential trend; right x-axis: long-term trend). The data are calculated from arbitrary time segments binned at 532 kyr. The balance line (0 %) indicates that the variance of response matches the variance in forcing, suggesting a self- sustained climate system, which is only paced by orbital cycles. Values much larger than zero would indicate a nonlinear reinforcement of the signal response. Note the different magnitude of amplification between the global LR04 δ¹⁸O (up to ~440%) and the equatorial Site 846 SST (up to ~180%) records. Negative R_s shows a response variance lower than orbital forcing, suggesting a damping mechanism for the climate system. The red box highlights the strong inverse coupling between obliquity damping and short eccentricity/precession amplification during post-MPT time. Data from Viaggi (2018b).

Figure 6. Pliocene-Pleistocene standardised PCA model for LR04 δ¹⁸O/Site 846 SST orbitals R_s and the mean climate state by time binned at 532 kyr: (a) PC-1 loadings; (b) PC-2 loadings; (c) PCs regression factor score (57.7% and 36.9% variance, respectively); (d) PCs factor score spread (PC-2 ‒ PC-1). Positive spread (blue bar): high-obliquity vs. low-short eccentricity/precession R_s. Negative spread (red bar): low-obliquity vs. high-short eccentricity/precession R_s. ONHG (Onset of the Northern Hemisphere Glaciation), INHG (Intensification of Northern Hemisphere Glaciation) and MBE label time bins containing these events. Blue arrows: transition patterns (TRA-1, 2, 3) of positive to negative factor spread including ONHG (TRA-1), INHG (TRA-2) and MBE (TRA-3), the latter containing the MPT. Data from Viaggi (2018b).

Figure 7. Time-series panel of δ¹⁸O SSA-component-2 (transition obliquity to short eccentricity) and its subcomponents and FFS of comp-2 for the three time intervals, highlighted by red labelled boxes (a, b, c). Note the transition in periodicity from obliquity (41-54 kyr) to short eccentricity (73-93 kyr) with a gradual upwards strengthening of the short eccentricity power. Interestingly, the amplitudes of the eccentricity in comp-2 subcomp-1-4, 9-15 increases (blue lines), whereas the amplitude of the obliquity in comp-2 subcomp-5-8 decreases during the MPT. X-axis is scale variable (expressed in ‰). The y-axis scale of the frequency spectrum increases upwards. ONHG = Onset of the Northern Hemisphere Glaciation; INHG = Intensification of Northern Hemisphere Glaciation; MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event. Figure is modified after Viaggi (2018a).

Figure 8. Time-series panel of the ODP Site 846 SST component-2 (transition obliquity to short eccentricity) and component-3-4 (obliquity) and FFS of the comp-2 for three time intervals, highlighted by red labelled boxes (a, b, c). The transition nature from obliquity (41-54-kyr) to short eccentricity (93-71-kyr) is less evident in this equatorial SST site, but still present. The increasing amplitude of the short eccentricity is clearer in subcomp-1-4 (blue lines). The obliquity damping during the MPT is very clear. X-axis scale variable (°C). The y-axis scale of the frequency spectrum increases upwards. ONHG = Onset of the Northern Hemisphere Glaciation; INHG = Intensification of Northern Hemisphere Glaciation; MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event. Figure is modified after Viaggi (2018a).

Figure 9. Time-series of obliquity components: (a) LR04 δ¹⁸O component-3-4; (b) ODP Site 846 SST component-3-4. The envelope of two or three low-amplitude 41-kyr obliquity peaks (glacial/interglacial) embedded in a weak ~74/92-kyr framework after the beginning of the MPT is noteworthy (red boxes and magnified time-series to the right). ONHG = Onset of the Northern Hemisphere Glaciation; INHG = Intensification of Northern Hemisphere Glaciation; MPT = Mid-Pleistocene Transition. Figure is modified after Viaggi (2018a).

Figure 10. Proxy records at Site U1308 (eastern North Atlantic): (a) wavelet analysis of benthic δ¹⁸O variation; (b) Ca/Sr and (c) Si/Sr records; (d) 41-kyr (blue) and 100-kyr (red) filtered components of the Si/Sr signal. Note the inverse amplitude coupling between obliquity and short eccentricity components. MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event. Figure is modified after Hodell et al. (2008).

Figure 11. Ca/Ti orbital signals at Site U1385 (Shackleton Site) on the SW Iberian Margin with bandpass filters for eccentricity (brown), obliquity (green), and precession (blue). Note the inverse amplitude coupling between obliquity (amplitude decrease) and short eccentricity/precession (amplitude increase) filtered components. MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event. Figure is modified after Hodell et al. (2015).

Figure 12. Global δ¹⁸O variability over the last 2 Ma: (a) benthic-planktic δ¹⁸O stack (thick line) of Huybers (2007) and benthic δ¹⁸O stack of Lisiecki and Raymo (2005); (b) evolutionary spectrum of benthic-planktic δ¹⁸O stack. Note the strong post-MPT increase in 100-kyr (and more moderate 22-kyr) spectral power coupled with a relevant decrease in the 41-kyr power; (c) Wavelet spectral power of detrended benthic LR04 is from Imbrie et al. (2011). Again, note the strong post-MPT increase in the 100-kyr spectral power coupled with decrease in the 41-kyr power. MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event. Figures are modified after (a, b) Huybers (2007), (c) Imbrie et al. (2011).

Figure 13. Wavelet power spectra of SST reconstructions from: (a) MD06-3018 (Coral Sea, subtropical SW Pacific); (b) MD97-2140 (Western Equatorial Pacific); (c) ODP 846 (Eastern Equatorial Pacific). The strong increase in 100-kyr spectral power is coupled to a decrease in the 41-kyr power. MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event. Figure is modified after Russon et al. (2011).

Figure 14. Benthic δ¹⁸O record from the eastern Mediterranean ODP Sites 967/968: (a) Morlet wavelet power spectra of the detrended δ¹⁸O signal; (b) Morlet wavelet power spectra of the δ¹⁸O 41-kyr component (filter centred on frequency 0.0241 ± 0.005 kyr^-1); and (c) δ¹⁸O 41-kyr time series. MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event. Data resampled at 1-kyr. Original data from Konijnendijk et al. (2015).

Figure 15. Morlet wavelet SST 41-kyr component (filter centred at frequency 0.0241±0.005 kyr^-1): (a) wavelet power spectra and time series from the North Atlantic ODP Site 982 (data from Lawrence et al., 2009); (b) wavelet power spectra and time series from the Arabian Sea ODP Site 722 (data from Herbert et al., 2010). MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event. Original data resampled at 1-kyr.

Figure 16. Morlet wavelet 41-kyr components (filter centred on frequency 0.0241±0.005 kyr^-1) at ODP Site 1143 (South China Sea). Wavelet power spectra and time series of (a) SST (data from Li et al., 2011); (b) δ¹⁸O planktic (data from Cheng X. et al., 2004); (c) δ¹⁸O benthic (data from Ao et al., 2011). MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event. Original data resampled at 1 kyr.

Figure 17. Variations in obliquity 40-kyr components of planktic δ¹⁸O and SST composite records from the sites RC11-120 and E49-18 (southern Indian Ocean) based on two different age models: (a) ELBOW; (b) TUNE-UP time scales. The reduction in amplitude of the obliquity cycles (red lines) over the last 400 kyr is highly evident. Figure is modified after Hays et al. (1976).

Figure 18. (a) Pliocene-Pleistocene time series of LR04 δ¹⁸O SSA-components: long-term trend (blue line, left y-axis, ~76% variance); short eccentricity (black line, ~6.5% variance) and obliquity (green line, 9.9% variance) (right y-axis); (b) short eccentricity and (c) obliquity Morlet wavelet spectra. Note the increase in amplitude/power of orbital responses (precession and half precession not shown) that is linked with long-term cooling (Viaggi, 2018a; this work, Fig.s 5 and 6) and the early onset (Piacenzian) of weak short eccentricity response. Blue arrows: subtrend I-II-III-IV, curvilinear step-wise isotopic enrichments of the long-term trend (Viaggi, 2018a). Red asterisks: multi-thresholds of the mean climate state that coincide with subtrends boundary (TH 1 to 4). The start of subtrend II and III is marked by ONHG and INHG, respectively. Subtrend IV includes MPT. Fitting of equations: full long-term trend (exponential); individual subtrends (linear). Maximum isotope enrichment rates for subtrends II and IV. Dashed black boxes are 532-kyr binned time intervals. Yellow arrows: transition patterns of positive to negative PCs spread suggesting competing interaction between obliquity and short eccentricity R_s under the influence of the long-term cooling (Figure 6); TRA-1 = Transition-1, including ONHG; TRA-2 = Transition-2, including INHG; TRA-3 = Transition-3, including MPT and MBE. The strengthening of the short eccentricity signal appears to be associated to the obliquity damped response whose maximum expression is post-MPT, and is preceded by an earlier and weaker expression since ONHG and INHG. Data from Viaggi (2018b).

Figure 19. Orbital nominal solutions and Morlet wavelet spectra over the last 2 Ma: (a) eccentricity (La2010 from Laskar et al., 2011); (b) obliquity; (c) precession (La2004 from Laskar et al., 2004). The amplitude trends (red lines) and the power spectra during the last 1.2 Ma exhibit signal pattern that is different from those of the climate proxies shown in this study, in particular concerning obliquity (nominal increase) and short eccentricity/precession (nominal decrease). MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event.

Figure 20. Theoretical secular changes of obliquity as a function of the ice-sheet response phase lag (degrees) for different positive/negative changes of oblateness ΔJ₂/J₂ (%). Theoretical constrains for ODH are highlighted in the coloured boxes (obliquity response lag <44° and positive/negative change of oblateness). Interglacial damping (high obliquity) requires positive change of oblateness, which leads to negative secular variation of obliquity; glacial damping (low obliquity) needs negative change of oblateness, which leads to positive secular variation of obliquity. Obliquity mean phase lags (*) and standard deviation intervals (vertical lines) are recalculated in the present study from different Pleistocene climate proxies (fuchsia: EPICA stack and SST; blue: benthic δ¹⁸O) for the last 800 kyr (Table 4). Figure is modified after Levrard and Laskar (2003).

Figure 21. Time series of the superposed SSA-filtered EPICA stack_s (red line, Viaggi, 2021a) and global benthic δ¹⁸O stack (blue line, Lisiecki and Raymo, 2005; Viaggi, 2018a) with cross-correlation statistics. Time series are standardised (0-mean, 1-SD) and δ¹⁸O inverted. The CCF exhibits the highest coefficient (0.91) at –5 lag number, equivalent to a δ¹⁸O ‘bulk’ lags of 2.5 kyr. The Pearson correlation at 0-lag is 0.87, the highest with respect to the single δD, CO₂ and CH₄ filtered records that are in the range of 0.80-0.83 (not shown). MPT = Mid-Pleistocene Transition; MBE = Mid-Brunhes Event.

Figure 22. Time series of standardised orbital SSA-components of the Red Sea RSL (present study, original data from Grant et al., 2014) and global δ¹⁸O benthic stack (Lisiecki and Raymo, 2005; Viaggi, 2018a) with nominal forcings (eccentricity La2010, Laskar et al., 2011; obliquity and precession La2004, Laskar et al., 2004) for the last 500 kyr: (a) short eccentricity; (b) obliquity; (c) precession. Nominal short eccentricity from Morlet wavelet filtering (frequency range 0.00678-0.01117 kyr^-1). δ¹⁸O data and nominal precession inverted to obtain the same paleoclimatic polarity. Records are plotted with their own time scale.

Figure 23. Plio-Pleistocene power spectral density (PSD) of (a) 65°N mean monthly insolation (21 june-20 july); (b) 85°N mean monthly insolation (21 june-20 july); (c) mean annual insolation. Colored lines (50%, 90%, 95%, 99%, and 99.9% significance levels) are first-order autoregressive AR(1) curves to check statistically significant peaks vs. red/white noise background. Nominal solutions from Laskar et al. (2004).

Figure 24. Panel of EPICA orbital stacks (Viaggi, 2021a) and related nominal solutions (eccentricity La2010, Laskar et al., 2011; obliquity, precession and mean annual insolation La2004, Laskar et al., 2004) for the last 800 kyr. Nominal short eccentricity components (La2010; MAI La2004) by filtering-out the 400-kyr band. Nominal precession inverted to align the climate polarity. EPICA stacks and nominal solutions are standardised. Records are with their own time scale.

Table 1. The maximum amplitude amplification (R_s max) and the maximum amplitude damping (R_s min) rate of change of the signals (EPICA δD, CO₂, CH₄ and LR04 δ¹⁸O) over the last 800 kyr, normalised in % for kyr by forcing (80-kyr bin). Note the marked asymmetry of the climate responses with extremely robust amplitude magnifications and moderate damping.

Rs statistic	Signal	Short eccentricity	Obliquity	Precession
Max	D	6.40	1.56	2.78
	CO2	6.18	4.21	2.58
	CH4	7.32	4.97	4.48
	18O	7.70	2.11	2.25
Min	D	–1.02	–0.45	–1.10
	CO2	–0.88	–0.97	–0.96
	CH4	–0.46	–1.08	–0.96
	18O	–0.04	–0.70	–0.97

Table 2. Results of quantitative estimation of the signal magnitude (percentage variance) and Fourier frequency data of the three Red Sea RSL reconstructed components isolated by singular spectrum analysis. The frequency peaks are the most significant peaks above a 99.9% critical limit. Original data are from Grant et al. (2014).

RSL component rank	RSL component variance	Frequency (kyr^-1)	TISA power (%)	Period (kyr)	Forcing
1-2;7-10	61.4 %	0.01001973	100.0	100	Short eccentricity
3-4;13-14	13.1 %	0.02460637	84.8	41	Obliquity
3-4;13-14	13.1 %	0.03244360	15.2	31	Obliquity
5-6;11-12	11.7 %	0.04393214	80.9	23	Precession
5-6;11-12	11.7 %	0.05342083	19.1	19	Precession
15-200	13.8 %				Suborbital + noise

Table 3. Comparison among the Red Sea RSL and EPICA Pliocene-Pleistocene rescaled variance vs. the LR04 δ¹⁸O. Δσ is the difference in variance between the RSL/EPICA and δ¹⁸O. Variance estimates are in terms of LR04 (Viaggi, 2018a) and EPICA (Viaggi, 2021a). Original RSL data are from Grant et al. (2014).

Forcing	Red Sea¹		Antarctica¹		LR04
Forcing	RSL (%)	Δσ (%)	EPICA (%)	Δσ (%)	δ¹⁸O (%)
Short eccentricity	14.5	8.0	12.2	5.7	6.5
Obliquity	3.1	–6.8	4.5	–5.4	9.9
Precession	2.8	0.8	2.0	0.0	2.0
¹ rescaled variance

Table 4. Cross-spectral analysis (Hann window) among obliquity components from different Pleistocene climate records (EPICA, SST, benthic δ¹⁸O) and La2004 obliquity nominal solution (Laskar et al., 2004) for the last 800 kyr. Analysis of uniform sampling interval (data resampled at 1-kyr, with the exception of EPICA 0.5-kyr). The data were detrended and standardised (0-mean, 1-SD) and, where necessary, inverted to obtain the same paleoclimatic polarity. Negative phase values indicate that the climate signal lags forcing. Records with their own time scale.

Obliquity component	Cross- spectrum freq. (kyr-1)	Cross- spectrum period (kyr)	Coherency	Phase shift (Deg, kyr)
EPICA D (0-800 kyr)1	0.02500	40.0	0.93	–37	–4.1
EPICA CO2 (0-800 kyr)1	0.02500	40.0	0.23a	–63	–7.0
EPICA CH4 (0-800 kyr)1	0.02500	40.0	0.76	–29	–3.2
EPICA stack (0-800 kyr)1	0.02500	40.0	0.79	–38	–4.2
SST East Equatorial Pacific ODP 846 (6-800 kyr)2	0.02500	40.0	0.59	–9	–1.0
SST North Atlantic DSDP 607 (250-2000 kyr)*3	0.02444	40.9	0.64	–34	–3.9
SST West Tropical Pacific IODP 1146 (6-800 kyr)*4	0.02500	40.0	0.86	–52	–5.7
SST Arabian Sea ODP 722 (8-800 kyr)*4	0.02500	40.0	0.83	–34	–3.8
18O benthic LR04 global stack (0-800 kyr)2	0.02500	40.0	0.88	–50	–5.5
18O benthic Atlantic stack (0-800)*5	0.02500	40.0	0.86	–52	–5.8
18O benthic Pacific stack (0-800)*5	0.02500	40.0	0.88	–51	–5.7
18O benthic East Mediterr. ODP 967/968 (0-800 kyr)*6	0.02500	40.0	0.69	–40	–4.4
Mean (EPICA stack + SST)	0.02489	40.2	0.74	–33 ± 15	–3.7 ± 1.7
Mean (benthic18O)	0.02500	40.0	0.83	–48 ± 6	–5.3 ± 0.6
* Morlet wavelet 41-kyr component (filter centred on main frequency 0.0241±0.005 kyr-1)

^a low coherency, ¹ obliquity SSA-component from Viaggi (2021b), ² obliquity SSA-component from Viaggi (2018b), ³ original data from Lawrence et al. (2010), ⁴ original data from Herbert et al. (2010), ⁵ original data from Lisiecki and Raymo (2009), ⁶ original data from Konijnendijk et al. (2015).

Table 5. Data of cross-spectral analysis (Hann window) among nominal forcings (eccentricity La2010, Laskar et al., 2011; obliquity and precession La2004, Laskar et al., 2004) and related SSA-signals from the Red Sea RSL (present study, original data from Grant et al., 2014) and global δ¹⁸O benthic stack (Lisiecki and Raymo, 2005; Viaggi, 2018b) for the last 500 kyr. δ¹⁸O data and nominal precession are inverted to obtain the same paleoclimatic polarity. Standardised records (0-mean, 1-SD). Uniform sampling interval at 1 kyr. Negative phase values indicate that the signal lags forcing. Records are with their own time scale.

Forcing	Signal	Cross-spectrum freq. (kyr^-1)	Cross-spectrum period (kyr)	Coherency	Phase shift (Deg, kyr)
Short eccentricity*	RSL	0.010	100.0	0.78	–12.6	–3.50	RSL lags short ecc.
	¹⁸O	0.010	100.0	0.75	–9.8	–2.72	¹⁸O lags short ecc.
	RSL vs.¹⁸O	0.010	100.0	0.94	2.8	0.78	¹⁸O leads RSL
Obliquity	RSL	0.024	41.7	0.88	–48.8	–5.65	RSL lags obl.
	¹⁸O	0.024	41.7	0.86	–66.1	–7.65	¹⁸O lags obl.
	RSL vs.¹⁸O	0.024	41.7	0.83	–17.3	–2.00	¹⁸O lags RSL
Precession	RSL	0.044	22.7	0.92	–46.3	–2.92	RSL lags prec.
	¹⁸O	0.044	22.7	0.88	–62.7	–3.96	¹⁸O lags prec.
	RSL vs.¹⁸O	0.044	22.7	0.92	–16.5	–1.04	¹⁸O lags RSL

* Morlet wavelet short eccentricity component (frequency range 0.00678-0.01117 kyr^-1).

Table 6. Data of cross-spectral analysis (Hann window) among nominal forcings (eccentricity La2010, Laskar et al., 2011; obliquity, precession and mean annual insolation La2004, Laskar et al., 2004) and related EPICA stacks (Viaggi, 2021a) for the last 800 kyr. Nominal precession inverted to align the climate polarity. Records rescaled in the range -1 to 1. Uniform sampling interval at 0.5 kyr. Negative phase values indicate that the signal lags forcing. Records are with their own time scale.

Forcing	EPICA signal	Cross-spectrum freq. (kyr^-1)	Cross-spectrum period (kyr)	Coherency	Phase shift (Deg, kyr)
Eccentricity	Short ecc. stack	0.01125	88.9	0.94	17.4	4.30	EPICA signal leads
Short eccentricity¹	Short ecc. stack	0.01125	88.9	0.94	18.3	4.52	EPICA signal leads
Mean annual insolation (short eccentricity¹)	Short ecc. stack	0.01	100.0	0.72	13.8	3.83	EPICA signal leads
Obliquity	Obliquity stack	0.025	40.0	0.98	-37.9	-4.21	EPICA signal lags
Precession	Precession stack	0.0425	23.5	0.67	-60.8	-3.97	EPICA signal lags
¹400-kyr band filtered out by wavelet analysis.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Global Evidence of Obliquity Damping in Climate Proxies and Sea-Level Record during the Last 1.2 Ma: A Missing Link for the Mid-Pleistocene Transition

Abstract

1. Introduction

1.1. Obliquity-oblateness feedback

1.2. Key role of Obliquity Forcing on the Earth’s Climate System

2. Materials and Methods

2.1. Materials

Global proxies

Antarctica proxies

Atlantic proxies

Pacific proxies

Mediterranean proxies

Indian proxies

2.2. Statistical methods

3. Results and Discussion

3.1. Evidence of post-MPT Obliquity Damping

3.1.1. EPICA record

3.1.2. Sea-level Record of the Red Sea

3.1.3. LR04 δ18O and equatorial ODP Site 846 SST

3.1.4. Atlantic, Pacific, Mediterranean, and Indian proxies

3.2. Long-Term Cooling sets Boundary Conditions for Glacial/Interglacial cycles

3.3. Amplitude relationships between Orbital Forcings and Proxies

3.4. Why does Obliquity’s Response Damping?

3.4.1. Remarks on Obliquity phase lag and Temperature/Ice-volume proxies

3.4.2. Observations of Orbital phase lags between δ18O and Red Sea RSL records

3.4.3. Changes in Earth’s Oblateness

3.5. Earth’s eccentricity and the 100,000-year issue

4. Summary and conclusions

4.1. Main results

4.2. Obliquity damping hypothesis

4.3. Outlook

Funding

Data Availability Statement

Acknowledgements

Conflicts of Interest

Abbreviations

Note

References

MDPI Initiatives

Important Links

Subscribe

3.1.3. LR04 δ¹⁸O and equatorial ODP Site 846 SST

3.4.2. Observations of Orbital phase lags between δ¹⁸O and Red Sea RSL records