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

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.