Viaggi Progress in Earth and Planetary Science
https://doi.org/10.1186/s40645-018-0236-z(2018) 5:81

Progress in Earth andPlanetary Science

PAPER WITH FULL DATA ATTACHEDOpen Accessδ18O and SST signal decomposition anddynamic of the Pliocene-Pleistocene
climate system: new insights on orbital
nonlinear behavior vs. long-term trend
Paolo Viaggi

Abstract
The global LR04 δ18O, the tropical ODP Site 846 sea surface temperature (SST), and the global ΔSST stack records
were investigated using the advanced method for time-series decomposition singular spectrum analysis to outline
the quantitative role of orbital forcings and to investigate the nonlinear dynamics of the Pliocene and Pleistocene
climate system. For the first time, a detailed quantitative evaluation is provided of the δ18O and SST variance paced
by long-period orbital modulation, short eccentricity, obliquity, precession, and half-precession cycles. New insights
into the nonlinear dynamic of the orbital components suggest considering astronomical signals as composite
feedback lagged responses paced by orbitals and damped (Early Pliocene) or amplified (Mid-Late Pleistocene) in arange of − 100 to + 400% the forcing. The Early Pliocene asymptotic decay of the δ18O and SST response sensitivity
up to − 100% observed for the first time in all orbital responses is interpreted as damping effect of a wide global
forest cover along with a possible high ocean primary productivity, through the CO2-related negative feedbacks
during time of global greenhouse. An anomalous post-Mid-Pleistocene Transition (MPT) sharply declines to near-zeroin obliquity response sensitivity observed in both global δ18O and tropical SST, suggesting an attenuation mechanismof the obliquity driving force and a reduction of the related feedback amplification processes. It is hypothesized thepost-MPT obliquity damping has contributed to the strengthening of the short eccentricity response by mitigating theobliquity “ice killing”, favoring a long-life ice sheet sensitive to a synergistic ~ 100-kyr amplification of positive feedback
processes during the time of a global icy state. The global δ18O, the tropical SST, and the global ΔSST trend
components, all explaining ~ 76% of the Plio-Pleistocene variance and significantly modifying the mean climate state,appear to be related to the long-term pCO2 proxies, supposedly controlled by plate tectonics through the global
carbon cycle (CO2 outgassing, explosive volcanism, orography and erosion, paleogeography, oceanic paleocirculation,and ocean fertilization). Finally, singular spectrum analysis provides a valuable tool in cyclostratigraphy with theremarkable advantage of separating full-resolution time series by variance strength.
Keywords: Singular spectrum analysis, Nonlinear, Chaotic, Damping, Paleoclimate, Milankovitch, Pliocene-Pleistocene,
Greenhouse

Correspondence: paolov6363@yahoo.it
Museo Geologico Giovanni Capellini, SMA, Bologna University, Via Zamboni
63, I-40126 Bologna, Italy© The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0
International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to
the Creative Commons license, and indicate if changes were made.

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Introduction
The Milankovitch theory elegantly describes the collective
effects of changes in the Earth’s movements on its climate.
Such orbital changes (eccentricity, obliquity, and precession) alter the amount and the distribution of solar radiation reaching the Earth. However, despite the studies ofKominz and Pisias (1979) and Wunsch (2004), in some
textbooks or papers, there is still the misleading notionthat the insolation changes are a major element controlling climate. A lot of studies have analyzed a variety of records finding the Milankovitch frequencies, but few havequantitatively investigated the hypothesis that most of theenergy conveyed in paleoclimate records is not related toorbital forcing (Kominz and Pisias 1979; Raymo 1994;
Wunsch 2004; Mudelsee and Raymo 2005; Hansen et al.
2013). This study applies a novel nested analytical approach based on the advanced method for signal decomposition singular spectrum analysis (SSA) with the aim of
separating full-resolution time series by signal variance,
the latter a useful “proxy” of the energy conveyed inpaleoclimate records. SSA is a relatively new and powerful
technique to processing time-series data, applicable tobroad areas such as climatology, marine science, geophysics, engineering, image processing, and medicine (Vautard
and Ghil 1989; Elsner and Tsonis 1996; Ghil et al. 2002;
Hassani 2007). However, the SSA approach is little used in
the study of paleoclimate records compared to conventional Fourier analysis, despite its capability to estimate
frequencies with higher resolution than Fourier approachand to detect signal structural changes, as highlighted in
the present study. Plio-Pleistocene long-term climatic
cooling and the increasing glacial cycle amplitude towards
the Late Pleistocene are evidence of significant changes in
the climate system dynamic described in the literature,
but the cause of these features remains poorly understood.
The aim of this study is to investigate these topics by decomposing two selected δ18O and sea surface temperature(SST) records to outline the quantitative role of astronomical forcings and the nonlinear behavior in the orbital responses and to assess the Earth’s long-term controlling
factors as a “stand-alone” climate system and its interaction with the astronomical one. Through the SSA, it
provides for the first time a detailed quantitative evaluation of the δ18O and SST signal variance in a wide
spectrum of orbital components from long-period orbital
modulation cycles to the half-precession, both unidentified in the previous studies of the original records. The
nonlinear dynamic of the astronomical climate responses
is discussed regarding new insights obtained by SSA
time-series processing. A link among the atmospheric
pCO2, the sea surface temperature, and the δ18O
long-term trend components is also provided to support
the hypothesis of a dominant non-astronomical control
on the Earth’s climate.

Page 2 of 37

Methods/Experimental
This work contributes to the research issues by applying
SSA to the global LR04 δ18O stack (Lisiecki and Raymo
2005b) and the tropical Ocean Drilling Program (ODP)
Site 846 SST record (Herbert et al. 2010b) (referred toas HT846 in this study), both high-resolution time series
covering ~ 5 million years to capture long-term climate
dynamics. The global LR04 benthic stack was obtained
from 57 oceanic series globally distributed over a wide
latitudinal range, including the Atlantic, Pacific, andIndian Oceans (Lisiecki and Raymo 2005a). Benthic δ18O
measures changes in global ice volume and deep watertemperatures, which are controlled by high-latitude surface temperatures (Lisiecki and Raymo 2007). The icevolume component lags the deep water temperature
signal, and the net phase of the δ18O record most likely isweighted towards the slow ice volume component(Lisiecki and Raymo 2005a; Herbert et al. 2010a; Berger et
al. 2016). Because the LR04 stack incorporates information from so many sites, it has a high signal-to-noise ratio
and more accurately reflects changes in global climate(Lisiecki and Raymo 2005a). The original resolution of theLR04 record ranges from 1 to 5 kyr within four time intervals with constant data spacing (Lisiecki and Raymo
2005a). The average resolution is 2.5 ± 1.4 kyr. The
HT846 is an alkenone-based tropical SST record locatedin the Eastern Equatorial Pacific (EEP). The ODP Site 846
is considered not “typical” tropical SST signal because itcools more than the ODP 662 and 722 sites, probably as
the result of upwelling changes in the EEP (Herbert et al.
2010a; Lisiecki 2010b); nevertheless, its high-resolution
SST record extending ~ 5 Myr makes it suitable to capture
both long-term dynamics and high-frequency cycles. The
original data spacing of the HT846 record is uneventhroughout the 5070 kyr of the time series (Herbert et al.
2010b), with an average value of 2.3 ± 1.2 kyr. SSA requires evenly spaced records. Therefore, both records
were resampled at constant time intervals of 1.0 kyr. Further details on the methodology used in this study and a
discussion of the quality and reliability of the time series
can be found in Additional file 1. The valuable global
ΔSST stack of Martinez-Boti et al. (2015) extends 3.4 Myr
with a resolution lower than that of HT846 (which is oneof the ten sites included in the stack). Thus, the ΔSST
stack was used in this study as a comparison record to assess some common features between the tropical Site 846
and the global SST records. The LR04 stack is orbitally
tuned to an ice model driven by June 21 insolation at
65° N from La93(1,1) orbital solutions (Laskar et al.
1993), with correction of the sedimentation rate (Lisiecki
and Raymo 2005a). The HT846 record is synchronizedusing benthic foraminiferal δ18O measured in the same
sediments at the study site and aligned to the LR04 stack(Herbert et al. 2010a). Unfortunately, high-resolution

Viaggi Progress in Earth and Planetary Science(2018) 5:81

extended records of CO2 proxies to study the long-term
dynamics are not currently available. The valuable simulated time series of Van de Wal et al. (2011) and Stap et al.
(2016) are not suitable for the purposes of the present
study because the data are not δ18O independent. Nevertheless, the long-term links among climate proxies can also
be studied with low-resolution records, provided that thedata are temporally extended and in sufficient number andquality to ensure reliable equation fitting. In this study, thecomposite atmospheric boron-based pCO2 proxy of the
ODP Sites 668B-999A (Hönisch et al. 2009; Bartoli et al.
2011) and the atmospheric alkenone-based pCO2 and SST
proxies of the ODP Site 999A (Seki et al. 2010) were selected.
Seki et al. (2010) proxies are rather poor in resolution but
cover a long period of time. The high-resolution atmospheric
CO2 reconstructions of Martinez-Boti et al. (2015), based on
multi-site boron-isotope records, were considered despite
the short time series from 3.3 to 2.3 million years ago. The
composite European Project for Ice Coring in Antarctica(EPICA) pCO2 record (Luthi et al. 2008; Berger et al. 2016)
and the δ13C-based pCO2 proxy (Lisiecki 2010b) are also
considered. In order to evaluate the trend relationships
among δ18O, SST, and pCO2 proxies, the high-resolution
LR04 and HT846 SSA-component time series were downsampled at the approximate age of the low-resolution pCO2
proxies from the ODP Sites 668B and 999A (Hönisch et al.
2009; Seki et al. 2010; Bartoli et al. 2011). By applying the advanced method for signal decomposition singular spectrum
analysis (Vautard and Ghil 1989; Elsner and Tsonis 1996;
Ghil et al. 2002; Hassani 2007), the 1.0-kyr resampled LR04
and HT846 records are partitioned in a new time series of its
data components. The SSA is achieved by the computationof forward-backward covariance matrices with an order of
50 (embedding dimension) using the singular value decomposition (SVD) procedure. This method for time-series analysis can be used to isolate full-resolution independent signal
components based on signal strength (variance). In contrast
with Fourier analysis with fixed basis of sine and cosine functions, SSA uses an adaptive basis generated by the time series
itself. As a result, the underlying model in SSA is more general and can extract amplitude-modulated wave components
with higher resolution than Fourier analysis. Moreover, SSA
can be effectively used as a non-parametric method for
detection of structural changes in the time series, reflecting
the main physical phenomena underlying the data (Vautard
and Ghil 1989; Hassani 2007). In the present study, the SST
and δ18O SSA-components capture important aspects of thenonlinear dynamic response of the climate system in avariance-oriented way. However, some high-variance δ18O
and SST signal components can remain characterized by astill complex frequency spectrum with some bands dispersed
in different components. A novel nested-SSA approach has
been applied to these 50-order components to better separate the frequency bands and correctly estimate the variance

Page 3 of 37

fraction, by increasing the embedding dimension up to 300
and partitioning the component into subcomponents. Moredetails on the SSA methodology and the complete
SSA-components dataset are provided in Additional files 1 and2. The SSA-components are then investigated by a Fourier
frequency spectrum (FFS) with a data window to analyze thefrequency power and to reduce spectral leakage in the Fourier spectra. The data tapering used in this study is thecs2-Hann window and the best exact N fast Fourier transform (FFT) algorithm. The peak detection uses peak-based
critical limit significance levels. In this type of confidence
limit, one seeks to disprove the null hypothesis where one
postulates either a white noise signal or a red noise signal.
The frequency peaks shown in this study are the most significant peaks above a 99.9% critical limit, unless otherwise
indicated. For each peak, the time-integral squared amplitude(TISA) power is computed in relative percent. The summed
power reported in the frequency analysis is merely the sumof the component powers, not the power of the composite
signal that would result from the addition of the components. Because clear non-stationarity in some components
complicates the Fourier representation, these signals are investigated by FFS for time intervals to capture dynamics overtime. Finally, by using the La93 orbital solutions (Laskar etal. 1993) (see Additional file 1), the phase relationships to therelated orbitally driven LR04-SST components were examined by cross-spectral analysis. The data were standardizedon the entire Plio-Pleistocene record (mean = 0 and standard
deviation = 1), and some time series were inverted in orderto have the same paleoclimatic polarity (positive forcing or
warming for positive standardized values). The software usedin this research is AutoSignal™.
The chronostratigraphic reference charts used in thisstudy are the ICS International Chronostratigraphic Chart(2017-02) for the Pliocene and Pleistocene and the ICS
global chronostratigraphical correlation table for the last
2.7 million years (v. 2016a) for the Quaternary.

Results
Analysis of δ18O SSA-components

The SSA allowed the resampled LR04 signal to be decomposed into five δ18O time series, which cumulatively explain 99.5% of the original variance (LR04-filtered). The
spectral framework of these δ18O components (abbreviated as comp.) is diversified, with fluctuations ranging
from 1330-kyr to 11-kyr periodicity (Fig. 1, Table 1). The
features of the δ18O components are described as follows.
δ18O component-1

This component explains 81.5% of the total variance (eigenvalue 10.89) and constitutes the most significant fraction of
the LR04 signal associated with a strong trend of isotopic
enrichment (long-term LR04 component, LT-LR04). The
FFS shows a broad spectrum of nonsignificant peaks

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Page 4 of 37

Fig. 1 Time-series panel of the five δ18O components isolated from the LR04 signal by singular spectrum analysis, together with the LR04-filtered
record (left) obtained by combining these data components (99.5% total variance). The x-axis scale (‰) is variable in relation to the variancefraction of the data components. Original LR04 δ18O data from Lisiecki and Raymo (2005b) resampled at 1 kyr. Time scale in kyr(< 50% of critical significance level) with periods equalto ~ 1.33 Myr (TISA power 77.2%), 409 kyr (10.8%),294 kyr (6.9%), and 195 kyr (5.1%), overlapping thegeneral trend (Fig. 1, Table 1). To separate the oscillatory terms from component-1 and estimate the variance fractions, an SSA was performed by increasing
the embedding dimension to 300. Three main subcomponents (abbreviated as subcomp.) were identified
equal to a reconstruction R2 of 0.999 (Table 2).
The component-1 subcomp-1 incorporates 96.2% of thecomponent-1 variance (81.5%), equivalent to an absolute
fraction of 78.4%. It exhibits the non-Milankovitch featureof a strong isotopic trend with mild oscillations and,
together with component-1, allows the subdivision of themain drift into four subtrends (Fig. 2).
These subtrends are characterized by a mild curvilinear
shape caused by a terminal acceleration of the isotopic enrichment, separated by a late warming step that temporarily
interrupts the drift. The first subtrend in chronologicalorder (subtrend-I) is in the range of 5.33–3.30 Myr(Zanclean-Early Piacenzian) and culminates with a major
event of isotopic enrichment at approximately 3.30 Myr, related to the first Pliocene transition of climatic deterioration
known as the onset of the Northern Hemisphere Glaciation(ONHG) (Shackleton and Opdyke 1977; Blanc et al. 1983;
Lourens and Hilgen 1994; Lisiecki and Raymo 2007).
This subtrend is interrupted by a rapid warming at

approximately 3.2 Myr, when a new isotopic drift(subtrend-II) begins and culminates at approximately
2.5 Myr (Middle Piacenzian to Gelasian) in a second
major climatic event known as the intensification of theNorthern Hemisphere Glaciation (INHG) (Zagwijn 1974;
Lisiecki and Raymo 2007; Seki et al. 2010). Because of its
global significance, this episode of climatic deterioration at
2.5 Myr was recently used to ratify the new Plio-Pleistocene
limit to the base of the Gelasian (GSSP 2.58 Myr) which isformally assigned to the Pleistocene (Gibbard and Head2009). Subtrend-III begins with warming at approximately
2.35 Myr and reaches its peak at 1.54 Myr at the end of theSanternian (Middle Gelasian-Santernian). The most recent
segment of isotopic enrichment (subtrend-IV) starts withthe warming at 1.45 Myr (Early Emilian) and culminateswith a noticeable acceleration at 0.65 Myr, close to the beginning of the Middle Pleistocene, correlating to Marine
Isotope Stage (MIS) 16. The Mid-Pleistocene Transition(MPT)—recently dated to start at approximately 1.2 Myr(Head and Gibbard 2005; Clark et al. 2006; Head et al.
2008)—which marks the beginning of the classic “Pleistocene Ice Age”, falls within subtrend-IV but is not directly
visible in this component. The end of the subtrend-IV
marks the final stage of the MPT at about 0.7 Myr (Clark et
al. 2006). After MIS-16, the cooling trend is broken by awide swing of δ18O depletion that peaks at approximately
0.40–0.30 Myr (MIS-11 and MIS-9e), referred in this studyViaggi Progress in Earth and Planetary Science(2018) 5:81

Page 5 of 37

Table 1 Main descriptive parameters and Fourier frequency results of the five δ18O components isolated from the LR04 signal by
singular spectrum analysis. The spectral peaks of the trend component-1 are not significant (< 50% of critical significance level) andare shown only for the purpose of discussion. Original LR04 data from Lisiecki and Raymo (2005b) are resampled at 1 kyr
δ18O component
rank

δ18O comp.
variance (%)

1

81.5

2

7.3

3–4

8.1

5–7

2.0

8–11

0.6

> 11

Frequency (kyr−1)

TISA power(%)

Period (kyr)Weighted mean
period (kyr)Forcing

0.0007519

77.2

1330

0.0024474

10.8

409

Not significant peaks(< 50% of criticalsignificance level)

Long-term trend + long-period
orbital modulation andeccentricity cycles

0.0034051

6.9

294

0.0051187

5.1

195

0.0105636

20.9

95

0.0140889

16.8

71

0.0186578

22.3

54

0.0244330

40.1

41

0.0244339

83.1

41

0.0253170

12.2

39

0.0186593

4.7

540.0422565

25.7

24

0.0447325

57.6

22

0.0487779

16.7

21

0.0700410

27.0

14

0.0764743

46.5

13

0.0893457

26.5

11

Transition obliquity to shorteccentricity

41

Obliquity

22

Precession

13“Half-precession” (equatorial)

0.5

Noise

as the Mid-Brunhes Oscillation (MBO), containing the
Mid-Brunhes Event (MBE) at approximately 0.43 Myr(Berger and Yin 2012). This subtrend pattern shows four
steps of gradual buildup of the NHG, followed by a relative
stasis of temperature cooling and ice growth during the last
600 kyr (MBO).

The frequency spectrum of comp-1 subcomp-1 stillshows a nonsignificant residual peak at 1.33 Myr, very
close to the ~ 1.2 Myr obliquity amplitude modulation
cycle (Laskar et al. 2004, 2011; Boulila et al. 2011),
tightly enveloped to the long-term trend, and to which
the subtrend’s mild terminal oscillations may be related

Table 2 Descriptive parameters and Fourier frequency results of the δ18O subcomponents isolated from the LR04 component-1 bysingular spectrum analysis. The processing of the subcomponents was conducted to better separate and estimate the variancefraction of specific frequency bands from still complex data components. The italicized percentage values are the absolute varianceof the subcomponents (referred to the total LR04 variance)
δ18O comp.
rank
1

δ18O comp.
variance (%)
81.5

δ18O
subcomp.
rank

Frequency (kyr−1)

TISA power(%)

Period (kyr)Forcing

78.4

0.0007519

96.6

1330

Long-term trend +obliquity modulation
cycle

1.0

0.8

0.0006351

94.9

1575

Obliquity modulation
cycle

0.0024431

5.1

409

Eccentricity (long)

2.8

2.3

0.0034061

16.0

294

Earth’s orbit secularfrequency

0.0051192

24.5

195

0.0061178

16.6

163

Obliquity modulation
cycle

0.0075215

15.6

133

0.0093870

12.8

107

0.0105656

14.5

95

δ18O
subcomp.
var. (rel.) (%)δ18O
subcomp.
var. (abs) (%)

1

96.2

2

3–10

Not significant peak(< 50% of criticalsignificance level)

Eccentricity (short)Viaggi Progress in Earth and Planetary Science(2018) 5:81

Page 6 of 37

Fig. 2 Time-series panel of the δ18O component-1 and its subcomponents isolated by singular spectrum analysis. δ18O comp-1 exponential trend(‰) = 1.44194 + 2.92378 e(−t/6963.268). The vertical brown arrows show four subtrends of component-1 and their correlation to the NHG transitions
known in the literature (MBO this study). To the right, the Fourier frequency spectrum of δ18O detrended comp-1 subcomp-3-11 (2.23% variance)
for three time intervals are highlighted by red boxes. Note the increase in amplitude after the ONHG and MPT and the power increase of the ~ 93-kyr
eccentricity cycle. x-axis scale variable (‰). The y-axis scale of the frequency spectrum increases upwards. Time scale in kyr(Fig. 2, Table 2). The comp-1 subcomp-2 contains 0.8%
only of the δ18O variance dominated by a 1.57-Myr periodicity (power 94.9%), which could be related to the 1.33
obliquity modulation cycle, and an ancillary term of
409 kyr. A wide spectrum of long-term cycles characterized the comp-1 subcomp-3-10, which explains 2.3% of
the δ18O variance (294, 195, 163, 133, 107, and 95 kyr).
To estimate the variance fraction of the “pure” long-term
trend, component-1 was detrended by subtracting an exponential fitting with least squares minimization (R2 = 0.94)
and the oscillatory signal was processed by a 300-order
SSA (Fig. 2, Table 3). The LT-LR04 exponential fit contains
93.8% of the comp-1 variance, equivalent to an absolutecontribution of 76.4%, whereas the detrended oscillatory
terms explain an absolute 5.1% variance. The 2% variancedifference between the comp-1 subcomp-1 (78.4%) and thecomp-1 exponential fit (76.4%) may be attributed to the
1.33-Myr cycle still superimposed to the comp-1 subcomp-1.
Considering the variance fraction of the comp-1 subcomp-2(0.8%, Table 2), a ~ 2.8% contribution of the 1.57–1.33-Myr
cycle may be estimated. Three subcomponents were isolated
from the detrended comp-1 (Table 3). The detrended
comp-1 subcomp-1 explains 2.25% of the δ18O variance andexhibits a spectrum dominated by a 3.75-Myr periodicity(power 95.8%) with an ancillary term of 855 kyr. Oddly, thecycle at about 1.57–1.33 Myr disappears in the detrended
component-1, while a new 3.75-Myr cycle is ascribable tothe 3.7-Myr g9-s9 terms of the Earth’s orbit secular frequency (Laskar et al. 2011). Prevailing 409-kyr (80.6%) and
294-kyr (19.4%) cycles characterized the detrended comp-1
subcomp-2, explaining only 0.61% of the δ18O variance.
These periodicities may be related to the astronomicallystable 405-kyr-long eccentricity cycle and to the 304-kyr
Earth orbit secular frequency, respectively (Laskar et al. 2004,
2011; Boulila et al. 2011).
The detrended comp-1 subcomp-3-11 explains 2.23%of the δ18O variance and shows a wide spectrum of frequency bands. A 294-kyr cycle (power 15.2%) may be related to the 304-kyr Earth orbit secular frequency.
Oscillations of 163 kyr (17.2%) and 195 kyr (24.9%) may
be linked to the long-period obliquity modulation cycles(~ 160–200 kyr) (Boulila et al. 2011). Finally, the shorteccentricity cycles of 133 kyr (15.9%), 107 kyr (12.7%),
and 95 kyr (14.1%) are also recognizable. The FFS of thissubcomponent by time intervals (Fig. 2) shows two time
horizons of amplitude shifting at about 3.32 MyrViaggi Progress in Earth and Planetary Science(2018) 5:81

Page 7 of 37

Table 3 Main descriptive parameters and Fourier frequency results of the δ18O subcomponents isolated from the detrended
component-1 and component-2 by singular spectrum analysis. The processing of the subcomponents was conducted to better
separate and estimate the variance fraction of specific frequency bands from still complex data components. The italicized
percentage values are the absolute variance of the subcomponents (referred to the total LR04 variance)
δ18O
comp.
rank
1

δ18O comp.
variance (%)
81.5

δ18O subcomp.
rank and var. (abs.)

2

7.3

δ18O subcomp.
var. (abs) (%)

Frequency(kyr−1)TISA power(%)

Period (kyr)Forcing

93.8

76.4

–

–

–

Long-term trend

1

44.2

2.25

0.0002668

95.8

3748

Earth’s orbit secularfrequency

0.0011701

4.2

855

Eccentricity (long)

2

11.9

0.61

0.0024443

80.6

409

0.0034017

19.4

294

3–11

43.8

2.23

0.0034062

15.2

294

0.0051199

24.9

195

0.0061180

17.2

163

0.0075215

15.9

133

0.0093867

12.7

107

exp.fitdetrend 6.2%
(5.1%)

δ18O subcomp.
var. (rel.) (%)

0.0105655

14.1

95

1–4, 9–15

75.3

5.5

0.0105637

56.6

95

0.0140872

43.4

71

5–8

24.4

1.8

0.0244321

67.5

41

0.0186582

32.5

54

(ONHG) and 1.24–0.70 Myr (MPT) that delineate threeintervals, in which occur a progressive increase in signal
amplitude and a relative strengthening of the short eccentricity cycle (93 kyr) at the expense of the obliquity
modulation cycles (195–205 kyr), revealing a transient nonlinear behavior of this small fraction of isotopic response.
δ18O component-2

Component-2 explains 7.3% of the total variance (eigenvalue 0.97) and exhibits four significant oscillatory components in the 41-kyr (power 40.1%) and 54-kyr (22.3%)
obliquity bands and in the 95-kyr (20.9%) and 71-kyr(16.8%) short eccentricity bands (Fig. 3, Table 1).
Two shifting horizons of isotopic amplitude (3.32 and1.24 Myr) are distinctly observed, correlating to the
ONHG and MPT start, respectively, and defining threetime intervals of increasing amplitude and a significant
change in orbital periodicity. The oldest segment between 5.32 and 3.32 Myr (Zanclean-Early Piacenzian)exhibits the lowest amplitude and power and a relative
prevalence of the obliquity cycle (41-kyr and 54-kyr).
The eccentricity bands at 89 kyr and 102 kyr are particularly weak (Fig. 3a). From 3.32 to 1.24 Myr (Early
Piacenzian-Emilian), a relative increase of the short eccentricity power (71-kyr and 91-kyr) over the obliquity
cycle occurs (Fig. 3b). The most recent segment from1.24 Myr exhibits a further increase in the average amplitude and an overwhelming power prevalence of the

Earth’s orbit secularfrequency
Obliquity modulation
cycle
Eccentricity(short)Eccentricity (short)

Obliquity

eccentricity 93-kyr cycle over the 41-kyr response, which
maintains about the same power (Fig. 3c). To further
isolate the two orbital terms from component-2, the SSA
was performed by increasing the embedding dimension to300, producing two subcomponents equal to a reconstruction R2 of 0.998 (Fig. 3; Table 3). Subcomponent-1-4, 9-15
explains 75.3% of the component-2 variance, equivalent toan absolute variance fraction of 5.5%. It exhibits frequencies in the short eccentricity band (95-kyr, power 56.6%;
71-kyr, 43.4%). Subcomponent-5-8 explains 24.4% of thecomponent-2 variance (1.8% absolute variance) in the obliquity frequency band (41-kyr, 67.5%, 54-kyr, 32.5%).
Interestingly, the amplitude of this obliquity’s subcomponent is decreasing after the MPT (Fig. 3). The change in
the nonlinear dynamic of the component-2 response (transition from obliquity to short eccentricity) is configured asa sort of “switching” of the orbital power and amplitude
between the 41-kyr obliquity and the 93-kyr eccentricity cycles. From this perspective, due to the change in periodicity
depending on the amplitude, component-2 may be the expression of a nonlinear phase-locking mechanism between
obliquity and eccentricity astronomical forcing (Tziperman
et al. 2006; Ruggieri et al. 2009). Notably, by modeling theNorthern Hemisphere ice volume over the last 3 Myr, Berger et al. (1999) simulated similar power spectrum changes
using a combination of insolation and linearly decreasing
CO2 concentration. Imbrie et al. (2011) simulated the same
spectral behavior using the Pleistocene LR04 record with aViaggi Progress in Earth and Planetary Science(2018) 5:81

Page 8 of 37

Fig. 3 Time-series panel of δ18O component-2 and its subcomponents and Fourier frequency spectrum of comp-2 for three time intervals a, b
and c, highlighted by red-labeled boxes. Note the evident increase in comp-2 amplitude after the ONHG and MPT, and the transition inperiodicity from a obliquity to c eccentricity with a gradual upwards strengthening of the short eccentricity power. Interestingly, the amplitudesof the eccentricity in comp-2 subcomp-1-4, 9-15 are increasing, whereas the amplitude of the obliquity in comp-2 subcomp-5-8 is decreasing
after the MPT. x-axis scale variable (‰). The y-axis scale of the frequency spectrum increases upwards. Time scale in kyr

new phase-space model of orbital forcings but with a state
threshold function of the ice volume, showing that ice volume plays an important role in pacing the orbital responseof the climate system. These simulated data, together withthe observed features of the SSA comp-2 and detrended
comp-1 subcomp-3-11 (Fig. 2), strengthen the evidence ofa paleoclimatological transition where the ONHG marks
the beginning of a weak and gradual sensibility of the climate system to a 93-kyr eccentricity forcing, supporting the
hypothesis of an ice volume CO2-induced control on eccentricity sensibility (Raymo 1998; Berger et al. 1999; Hansen et al. 2007; Lisiecki and Raymo 2007; Imbrie et al. 2011;
Berger and Yin 2012).
δ18O component-3-4This signal explains 8.1% of the δ18O variance (eigenvalue 1.08) with a 95.3% power at an obliquity-driven
41–39-kyr period (Table 1). Component-3-4 is characterized by a clear increase in the average 41-kyr amplitude and power towards the beginning of the MPT(Fig. 4a, b), followed by an amplitude stabilization andpower reduction during the MPT and post-MPT (Fig. 4c).
Nonstationary behavior is observed at the beginning of
the MPT, where a change in the shape of the 41-kyr

cycle that appears enveloped in a very weak 93-kyr signal
can be noted (Fig. 4c), suggesting an apparent influence
between the eccentricity and the obliquity cycles after
the beginning of the MPT. The ONHG transition at
3.32 Myr does not seem to be associated with any particular features of component-3-4, except for a strong
power increase in the 41-kyr cycle after the transition.
δ18O component-5-7This component explains 2.0% of the δ18O variance(eigenvalue 0.27) driven by precession-related bands at~ 22 kyr (Table 1). This precession δ18O component also
shows the nonlinear shape of the upwardly increasing
isotopic amplitude (Fig. 1).
δ18O component-8-11

Component-8-11 exhibits 0.6% of the δ18O variance(eigenvalue 0.07) and contains a weight mean period at~ 13 kyr (Table 1). Furthermore, the feature of increasing mean amplitude in this weak component is still observed (Fig. 1). From an interpretative viewpoint,
component-8-11 is related to the “half-precessional” cycleof equatorial origin modulated at ~ 11 kyr (Berger et al.
2006). Hagelberg et al. (1994) found precession-relatedViaggi Progress in Earth and Planetary Science(2018) 5:81

Page 9 of 37

Fig. 4 Time-series panel of the δ18O component-3-4 (8.1% var.) and component-2 subcomp-5-8 (1.8% var.) related to obliquity. The Fourier frequency
spectrum of component-3-4 in three time intervals is highlighted by red-labeled boxes. The spectrum is always dominated by obliquity 41-kyr power.
Note a clear increase in the average 41-kyr amplitude and power towards the beginning of the MPT (box 4a and 4b), followed by an amplitude
stabilization and power reduction (box 4c) and an amplitude reduction (comp-2 subcomp-5-8) during the MPT and post-MPT. The envelope of two to
three obliquity cycles in a weak eccentricity 93-kyr framework after the beginning of the MPT is noteworthy (box 4c and magnified time series to theright). FFS scale power increases upwards from a to b and decreases from b to c. x-axis scale variable (‰). Time scale in kyr

climate variability at periods from 10 to 12 kyr in threesites for the Late Pleistocene (ODP 846, ODP 663, andDSDP 609). They concluded that this variability may bederived from the high sensitivity of the tropics to summertime insolation in both hemispheres, triggering an amplified response of tropical precipitation and temperaturethat may be transmitted to high latitudes via Atlantic
advective transport. The results of magnetic susceptibility and particle size analysis from the Northwestern
Chinese Loess Plateau reveal well-defined half-precessional
cycles during the last interglaciation, interpreted as a direct
response to low-latitude forcing through its modulation of
the East Asian summer monsoon (Sun and Huang 2006).
However, despite this interesting evidence of thesemi-precessional cycle of equatorial origin in the δ18O
deep-temperature and ice-volume signal, its contribution to
the total variance of the global δ18O record is less than 1%.
Analysis of the SST SSA-components

The resampled HT846 record was partitioned by SSA
into five SST time series, which cumulatively explain
99.6% of the original variance (HT846-filtered) (Fig. 5,
Table 4). The spectral frameworks of these SST components, ranging from 1266-kyr to 11-kyr periodicity, arevery similar to that of the LR04 record. The features of
the SST components are described and compared to theδ18O ones as follows.
SST component-1

This component explains 87.8% of the original variance(eigenvalue 128.1) and contains the most significant fraction of the SST signal associated with a strong cooling
trend (long-term HT846 component, LT-HT846). The
FFS shows a broad spectrum of significant peaks with periods of ~ 1.27 Myr (TISA power 83.3%), 406 kyr (8.6%),
298 kyr (5.0%), and 193 kyr (3.1%), overlapping the generaltrend (Fig. 5, Table 4). A high-order SSA was performed
to separate the oscillatory terms from the component-1
trend, from which three subcomponents were identified
equal to a reconstruction R2 of 0.999 (Table 5).
The component-1 subcomp-1 explains an absolute
variance of 81.7% and exhibits a strong cooling trendwith significant mild oscillations. The four LR04s’ subtrends are generally identifiable in the tropical SST
long-term components, with some differences especiallyin subtrend-III due to a strong cooling event at ~ 2.1 Myr(Herbert et al. 2010a), which is less evident in theδ18O signal (Fig. 6). In addition, the MBO feature isViaggi Progress in Earth and Planetary Science(2018) 5:81

Page 10 of 37

Fig. 5 Time-series panel of the five SST components isolated from the ODP Site 846 record by singular spectrum analysis with the SST-filtered
record (left) obtained by combining these data components (99.6% total variance). The x-axis scale (°C) is variable in relation to the component’s
variance fraction. Original SST data from Herbert et al. (2010b) resampled at 1 kyr. Time scale in kyr

less marked in the SST signal. These differences between
deep-temperature/ice-volume and sea surface temperature
signals could be due to temporal effects of upwelling in
the EEP “cool” region (Lisiecki 2010b).
The frequency spectrum of the SST comp-1 subcomp-1
still shows a significant peak of 1.27 Myr period (98.7%
power), related to the ~ 1.2-Myr obliquity modulation
cycle (Laskar et al. 2004, 2011; Boulila et al. 2011), towhich the mild subtrend oscillations may be related(Table 5). The comp-1 subcomp-2 explains 3.0% of theHT846 variance dominated by a 1.50-Myr periodicity(94.1%), probably related to the same obliquity modulation
cycle. A wide spectrum of long-term cycles characterized
the comp-1 subcomp-3-10, explaining 3.1% of the original
SST variance (399, 292, 213–193, 126, and 95 kyr).
To estimate the variance of the SST long-term trend, thecomponent-1 was detrended by subtracting an exponential
fitting with least squares minimization (R2 = 0.87), and theoscillatory signal was processed by a high-order SSA (Fig. 6,
Table 6). The LT-HT846 exponential fit explains 87.2% of
the comp-1 variance, equivalent to an absolute contributionof 76.6%. The detrended oscillatory signal explains an absolute 11.2% variance. The 5.1% variance difference between
the comp-1 subcomp-1 (81.7%) and the comp-1 exponential fit (76.6%) may be attributed to the 1.27-Myr cycle still
superimposed on the comp-1 subcomp-1. Three subcomponents were isolated from the SST detrended comp-1(reconstruction R2 = 0.999). The detrended comp-1

subcomp-1 explains 5.56% of the SST original variance andexhibits a spectrum dominated by 1.37-Myr (power 51.9%)
and 3.6-Myr (40.6%) cycles, with an ancillary term of
827 kyr. Again, the 3.7- to 3.6-Myr cycle appears in theLR04 and HT846 detrended comp-1 subcomp-1 only;
however, in the HT846 record, the 1.37-Myr cycle does notdisappear. Therefore, the 5.56% tropical SST variance(detrended comp-1 subcomp-1) related to the long-term
cycles (1.37-Myr obliquity modulation cycle and 3.6-Myr
Earth orbit secular frequency) is consistent with the 5.1%
variance estimate by difference from comp-1 subcomp-1(which includes the 1.27-Myr cycle only).
Prevailing 405-kyr (70.3%) and 300-kyr (29.7%) cycles
characterize the detrended comp-1 subcomp-2 explaining
2.81% of the SST variance. These frequencies may be related to the astronomically stable 405-kyr-long eccentricity
cycle and to the 304-kyr Earth orbit secular frequency(Laskar et al. 2004, 2011; Boulila et al. 2011). The detrended
comp-1 subcomp-3-11 explains 2.84% of the SST varianceand shows a wide spectrum of frequency bands. A 295-kyr
cycle (power 15.4%) can be attributed to the 304-kyr Earthorbit secular frequency. Oscillations of 193 kyr (45%)
and 148 kyr (6.7%) may be linked to the long-period obliquity modulation cycles (~ 160–200 kyr) (Boulila et al.
2011). Finally, weak short eccentricity cycles of 136 kyr(9.2%), 126 kyr (14%), and 95 kyr (9.7%) are also
recognizable. The FFS of this subcomponent by time intervals (Fig. 6) shows a relative strengthening of the short(2018) 5:81

Viaggi Progress in Earth and Planetary Science

Page 11 of 37

Table 4 Main descriptive parameters and Fourier frequency results of the five ODP Site 846 SST components isolated by singular
spectrum analysis. Original SST data from Herbert et al. (2010b) resampled at 1 kyr
SST componentrank

SST comp.
variance (%)

Frequency(kyr−1)TISA power(%)

Period (kyr)87.8

0.0007899

83.3

1266

0.0024651

8.6

406

0.0033527

5.0

298

0.0051684

3.1

193

0.0105090

32.6

95

0.0110927

26.3

90

1

2

5.2

3–4

4.5

5–7

1.6

8–11

0.5

> 11

Weighted mean period(kyr)ForcingLong-term trend + long-period
orbital modulation andeccentricity cycles

Transition obliquity to shorteccentricity

0.0144249

41.1

69

0.0244798–

41

Very low power peak
40

Obliquity

22

Precession

12“Half-precession” (equatorial)

0.0244852

69.3

41

0.0252735

18.4

40

0.0268223

12.3

37

0.0439726

31.6

23

0.0463271

43.1

22

0.0487111

25.3

21

0.0732060

36.7

14

0.0767562

30.5

13

0.0946970

32.8

11

0.4

Noise

eccentricity cycle (93-kyr) at the expense of the obliquity
modulation cycles (~ 195- to 158-kyr), whereas the increase in signal amplitude by time is not observable
post-MPT.
SST component-2

This component explains 5.2% of the total variance(eigenvalue 7.56) and exhibits a basically 95- to 90-kyr
short eccentricity periodicity (58.9%) with a 69-kyr term

(41.1%) and a very weak but significant 41-kyr obliquity
oscillation (Fig. 7; Table 4). Compared to the δ18O
comp-2, the two shifting horizons of increasing amplitude (ONHG and MPT start) and the change in orbital
periodicity from 41–54 to ~ 93 kyr by three time segments are less evident (Fig. 7a–c). The transition
from obliquity to short eccentricity is less evident but
still present as shown also by the subcomponent analysis (Table 6). Three subcomponents were isolated by

Table 5 Descriptive parameters and Fourier frequency results of the ODP Site 846 SST subcomponents isolated from the SST
component-1 by singular spectrum analysis. The processing of the subcomponents was done to better separate and estimate thevariance fraction of specific frequency bands from still complex data components. The italicized percentage values are the absolute
variance of the subcomponents (referred to the total SST variance)
SST comp.
rank
1

SST comp.
variance (%)
87.8

SST subcomp.
var. (rel.) (%)

SST subcomp.
var. (abs) (%)

Frequency(kyr−1)

TISA
power (%)

Period(kyr)

Forcing

1

93.0

81.7

0.0007899

98.7

1266

Long-term trend + obliquity
modulation cycle

2

3.4

3.0

0.0006678

94.1

1497

Obliquity modulation cycle

0.0024657

5.9

406

Eccentricity (long)

3–10

3.5

3.1

0.0025065

8.2

399

0.0034190

18.5

292

Earth’s orbit secularfrequency

0.0046896

15.3

213

Obliquity modulation cycle

0.0051832

36.9

193

0.0079507

12.7

126

0.0105043

8.4

95

SST subcomp.
rank

Eccentricity (short)Viaggi Progress in Earth and Planetary Science(2018) 5:81

Page 12 of 37

Fig. 6 Time-series panel of the ODP Site 846 SST component-1 and its subcomponents isolated by singular spectrum analysis. SST comp-1
exponential trend (°C) = 33.58041–11.84075 e(−t/8830.503). The vertical brown arrows show four subtrends of component-1 and their correlation to
the NHG transitions known in the literature (MBO this study). To the right, Fourier frequency spectrum of SST detrended comp-1 subcomp-3-11(2.84% variance) for three time intervals are highlighted by red boxes. Note the increase in amplitude after the ONHG (but not after the MPT) andthe power increase of the ~ 93-kyr eccentricity cycle. x-axis scale variable (°C). The y-axis scale of the frequency spectrum increases upwards. Timescale in kyr

a high-order SSA (R2 = 0.999). Subcomponent-1-4 explains 3.65% of the original SST variance and shows afrequency spectrum in the short eccentricity band (95-kyr,power 43.2% and 69-kyr, 56.8%). Subcomponent-5-8,
11-15 explains 1.27% variance in a still complex transition
spectrum from the obliquity (53-kyr, 60.2%) to eccentricity
band (126-kyr, 39.8%). Finally, subcomp-9-10 explains
only 0.27% of the SST variance in the obliquity frequency
band (45-kyr, 62.2%; 41-kyr, 37.8%).
SST component-3-4This component explains 4.5% of the SST variance(eigenvalue 6.56) at the 87.7% obliquity-driven 41–
40-kyr period (Table 4). The SST component-3-4 isstill characterized by a less evident increase in the
average amplitude by time approaching the beginningof the MPT (Fig. 8a, b), followed by a clear decreasein amplitude after the MPT start (Fig. 8c). A slightchange in the shape of the 41-kyr cycle that appears
enveloped in a very weak ~ 93-kyr signal can be observed after 1.2 Myr (Fig. 8c).

SST component-5-7This component explains 1.6% of the SST variance(eigenvalue 2.24) driven by precession frequency bandsat ~ 22 kyr (Table 4), similar to the precession δ18O
component.
SST component-8-11

Component-8-11 explains a total 0.5% of the SST variance (eigenvalue 0.74) and shows a weight mean periodat ~ 12 kyr (Table 4). An increasing mean amplitude in thisweak component is still observed (Fig. 5). This component
is similar to that of the δ18O identified in this study and canbe related to the “half-precessional” cycle of equatorial
origin modulated at ~ 11 kyr (Berger et al. 2006).
Summary of δ18O and SST forcing variance contribution

The SSA results allow the LR04 δ18O and ODP Site 846
SST forcing variance contribution to be estimated quantitatively in detail. A complex envelope of long-term orbital cycles accounts for a maximum of 5.1% of the δ18O variance(due to a small fraction of short eccentricity in detrendViaggi Progress in Earth and Planetary Science(2018) 5:81

Page 13 of 37

Table 6 Main descriptive parameters and Fourier frequency results of the ODP Site 846 SST subcomponents isolated from thedetrended comp-1 and comp-2 by singular spectrum analysis. The processing of the subcomponents was made in order to better
separate and estimate the variance fraction of specific frequency bands from still complex data components. The italicized
percentage numbers are the absolute variance of the subcomponents (referred to the total SST variance)
SST comp. SST comp.
SST subcomp.
rank
variance (%) rank and var. (abs.)
1

87.8

exp.fit

87.2

76.6

–

–

–

Long-term trend

detrend 12.8% 1(11.2%)

49.6

5.56

0.0002758

40.6

3626

Earth’s orbit secularfrequency

0.0007313

51.9

1367

Obliquity modulation
cycle

0.0012092

7.5

827

Eccentricity (long)

0.0024667

70.3

405

0.0033343

29.7

300

0.0033896

15.4

295

0.0051716

45.0

193

0.0067394

6.7

148

0.0073694

9.2

136

0.0079678

14.0

126

0.0104995

9.7

95

2

3–11

2

5.2

SST subcomp SST subcomp Frequency (kyr−1) TISA
Period (kyr) Forcing
var. (rel.) (%)
var. (abs) (%)
power (%)

1–4

25.1

25.4

70.2

2.81
2.84

3.65

0.0105084
0.0144299

5–8, 11–15

24.5

1.27

0.0079524
0.0187962

9–10

5.2

0.27

0.0220179
0.0244820

comp-1 subcomp-3-11; Table 3). Within the 2.23% varianceof the δ18O detrend comp-1 subcomp-3-11, a ~ 0.95% shorteccentricity contribution may be roughly estimated byusing the TISA power of the short eccentricity terms(42.7% of 2.23% variance). This approach is justified because of the low variance of the detrend comp-1
subcomp-3-11, causing a negligible error. Therefore, δ18O
long-period orbital modulation cycles (1.57–1.33 Myr and~ 160–200 kyr obliquity amplitude modulation; 3.75 Myrand ~ 300 kyr Earth orbit secular frequencies) and the ~405-kyr-long eccentricity cycles tightly overlapped to theδ18O trend, influencing the Earth’s climate with an estimated variance of ~ 4.2%. Regarding the tropical SST, asimilar complex envelope of long-term orbital cycles account for a maximum 11.2% of the HT846 variance(Table 6), and by subtracting the small variance of the shorteccentricity term (32.9% of 2.84%, variance = 0.93%), a contribution of ~ 10.3% can be estimated.
Long-term cycles (1.5–1.0 Myr) were previously observedin the planktonic δ18O record of ODP Site 625B from theGulf of Mexico and interpreted as a low-frequency component of the climate system, likely unrelated to orbital forcing(Joyce et al. 1990). Instead, according to Lourens and Hilgen(1994), the long-term 1.3-Myr period observed from some
Mediterranean and equatorial Atlantic and Pacific sites

43.2

95

56.8

69

39.8

126

60.2

53

62.2

45

37.8

41

Earth’s orbit secularfrequency
Obliquity modulation
cycle
Eccentricity (short)

Eccentricity (short)Transition obliquity
to short eccentricity
Obliquity

resulted from a long-term obliquity modulation cycle.
Through a correlation between Cenozoic icehouse sequences
and astronomical cycles, Boulila et al. (2011) demonstrate aclose correlation between ~ 1.2-Myr obliquity modulation
cycles and third-order eustatic sequences, suggesting glacioeustasy as the major driving process of third-order sequences in an icehouse world. A complex envelope of
Plio-Pleistocene long-term astronomical cycles and especiallythe 1.57–1.33-Myr obliquity amplitude modulation to which
the mild NHG’s subtrend terminal oscillations may be relatedby transient interaction on the long-term trend are further
geological evidences of the chaotic behavior of the
Solar System (Laskar et al. 2004, 2011; Boulila et al.
2011; Westerhold et al. 2017), and the SSA results
show that their influence on Earth’s climate, estimatedfor the first time in this study, have important differences in variance (global δ18O ~ 4.2%; tropical SST ~
10.3%). Interestingly, the tight overlap of these long
cycles could be the reflection of a close link between
these weak orbital forcings and the much more effective
forcing at the origin of the long-term trends through atemperature-ice volume-induced control on the long-cycle
sensibility. The lack of the ~ 1.2-Myr eustatic cycles in the
Mesozoic greenhouse world (ephemeral or no ice sheets)(Boulila et al. 2011) may support this hypothesis.

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Page 14 of 37

Fig. 7 Time-series panel of ODP Site 846 SST component-2 and its subcomponents and Fourier frequency spectrum of comp-2 for three timeintervals a, b and c, highlighted by red-labeled boxes. The transition nature from a obliquity to c short eccentricity is less evident in this tropical
SST site, but still present. The increasing amplitude of the short eccentricity is more clear in subcomp-1-4. The obliquity decrease after the MPT is
still shown by subcomp-9-10. x-axis scale variable (°C). The y-axis scale of the frequency spectrum increases upwards. Time scale in kyr

The short eccentricity does not manifest itself as a“pure” δ18O component; rather, it appears to move gradually to the Pleistocene climate signal as an interaction
subcomponent only, with a long-term response (detrend
comp-1 subcomp-3-11, ~ 0.95% variance) and theobliquity-eccentricity transitional component-2 (component-2 subcomp-1-4, 9-15, 5.5% variance). Therefore,
the short eccentricity cycle contributes ~ 6.5% to theglobal Plio-Pleistocene δ18O signal. Similarly, the SST
detrend comp-1 subcomp-3-11 explains ~ 0.93% of shorteccentricity variance, and a 3.65% fraction is recorded in
the SST comp-2 subcomp-1-4. A small fraction of shorteccentricity terms are still recorded in the SST comp-2
subcomp-5-8, 11-15 for which a 0.5% variance can be
estimated (39.8% of 1.27% variance), resulting in a SST
short eccentricity total variance of ~ 5.1%. These values
agree with the consideration that the amount of insolation perturbation at the eccentricity-driven ~ 100-kyr
period is too small to directly cause climate change ofan ice age amplitude (Imbrie et al. 1993; Berger et al.
1999, 2006; Lisiecki 2010a).
The 41-kyr obliquity cycle appears to transfer to theclimate system both as a “pure” δ18O component, withsignal strength of 8.1% (component-3-4), and as aninteraction subcomponent like the obliquity-eccentricity

transitional component-2 (comp-2 subcomp-5-8, 1.8%
var.). Therefore, the total δ18O signal variance in the
41-kyr obliquity band is estimated to be 9.9%. The
SST shows a similar behavior characterized by a“pure” component-3-4 (4.5% variance) plus transitional obliquity-eccentricity subcomponents (comp-2
subcomp-9-10, 0.27% var.; comp-2 subcomp-5-8, 11–15,~ 0.76%), from which a total obliquity contribution of
5.5% may be estimated.
The precession and semi-precession δ18O variance canbe evaluated at about 2.0% and 0.6%, close to the SST
contribution of 1.6% and 0.5%, respectively. In summary,the variance of the δ18O orbital components related tolong-term, short eccentricity, obliquity, precession, andsemi-precession cycles can be evaluated at approximately 4.2%, 6.5%, 9.9%, 2.0%, and 0.6%, respectively(total 23.2% variance), whereas the orbital SST components account for 10.3%, 5.1%, 5.5%, 1.6%, and 0.5%, respectively (total 23.0% variance). Within orbital forcings,
the impact of the obliquity cycle appears to be the mostconsistent in the LR04 stack. In the tropical HT846 record, the most consistent orbital signal is related to thelong-term cycles. These results agree with those of
Kominz and Pisias (1979), which state that the fractionof the record variance related to orbital changes is less

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Page 15 of 37

Fig. 8 Time-series panel of the ODP Site 846 SST obliquity’s component-3-4 (4.5% var.) and component-2 subcomp-9-10 (0.27% var.). Fourier frequency
spectrum of comp-3-4 by three time intervals highlighted by red labeled boxes. The spectrum is always dominated by obliquity 41-kyr power. Note a clear
increase in the average 41-kyr amplitude and power towards the beginning of the MPT (box 8a and 8b), followed by a clear decrease in amplitude after
the beginning of the MPT (box 8c and comp-2 subcomp-9-10). The envelope of two to three obliquity cycles in a weak eccentricity ~ 93-kyr framework
after the beginning of the MPT is less evident (box 8c and magnified time series to the right). FFS scale power increasing from a to b and was constantfrom b to c. x-axis scale variable (°C). Timescale in kyr

than 25%, and of Wunsch (2004), which states that thisfraction never exceeds 20%. Interestingly, the cumulative
influence on the Earth’s climate of orbital forcings
estimated from the global δ18O stack is equal to that of
the tropical SST single record, while the individual contribution of each astronomical term is quite different.
This result could be attributed to the latitudinal gradients of orbital forcings (Berger 1978; Liu and Herbert
2004), while the constancy of the long-term variancemight reflect a global latitude-independent nature of thepaleoclimatic trend’s control factors. The fact that theδ18O and SST orbital components explain ~ 23% of thetotal variance clearly shows that astronomical forcing
does not play a major role in guiding the response of theclimate system in the Plio-Pleistocene, in agreementwith Kominz and Pisias (1979), Raymo (1994), Wunsch(2004), and Mudelsee and Raymo (2005). However,~ 76% of the Plio-Pleistocene climate variance is related to
non-Milankovitch long-term controlling factors, drives
approximately three fourths of the energy conveyed inpaleoclimate records, and thus constitutes the main
Earth’s climate system. This robust and stable trend signal
needs to be explained.

Statistical analysis of the mean climate state through time

The SSA time series allow a quantitative analysis of themean climate state through time. The Plio-Pleistocene
δ18O and SST components and the solution of orbital
forcings (Laskar et al. 1993) were standardized (mean = 0,
standard deviation = 1) and binned at three arbitrary constant time intervals (532 kyr, 760 kyr, and 1064 kyr) to
avoid data density distortion in descriptive statistics andsignificance tests of mean differences among groups(Table 7). Some data were inverted in order to have thesame paleoclimatic polarity (Fig. 9).
As observed in the comparison of the 95% confidence
interval (CI) for the means (Fig. 9a, c) and in the significance tests of the mean differences in the analysis ofvariance (ANOVA) table (Table 7), the δ18O and SST
trend components (comp-1, comp-1 subcomp-1, andcomp-1 exponential fit) can very significantly change thesignal mean for different time intervals. Instead, theδ18O and SST orbital components (comp-2-11) and theastronomical forcings modified the mean with little significance or no significance. These results are similar
across all three binned time intervals showing that thisfeature is not a time frame artifact (Table 7). In moreViaggi Progress in Earth and Planetary Science(2018) 5:81

detail, the mean reaches a deviation from the standardized mean up to − 1.51–1.18 standard units through thetime segments in the δ18O component-1, up to − 1.50–
1.21 in the δ18O comp-1 subcomp-1, and in the range −1.72–1.40 in the δ18O comp-1 exponential fit. The SST
trend components exhibit similar results with deviations
up to − 1.45–1.27, − 1.49–1.33, and − 1.67–1.50 for thecomp-1, comp-1 subcomp-1, and comp-1 exponential
fit, respectively. The δ18O and SST orbital signals reachmaximum values of − 0.08–0.09 and − 0.09–0.04 standard
deviations, respectively. Eccentricity reaches a maximum
deviation from the standardized mean of − 0.30–0.28 units.
Obliquity ranges from − 0.12 to 0.12 and precession from− 0.02 to 0.01. The insolation mean (65° N at June 21) oscillates in the range of − 0.08–0.06. The ANOVA analysis was
also done to all single orbital components, obtaining similar
results. The histograms of the Plio-Pleistocene δ18O andSST astronomical components-2-11, compared to that of
the trend component-1 subcomp-1 (Fig. 9b, d), show aroughly normal distribution centered on the orbital response mean, while the trend component exhibits a multimodal distribution away from the standardized mean. This
analysis suggests Plio-Pleistocene orbital forcing could beconsidered a random walk at zero-sum (the standardizedmean): on long time scales, astronomical forcing oscillates
around a basically time-invariant normally distributedmean of the Plio-Pleistocene climate system response.
These results, consistent with those of the quantitative estimation of the variance, suggest that changes in the mean
climate state during the Plio-Pleistocene were primarily
caused by non-Milankovitch first-order controlling factors
recorded in the global LR04 and tropical HT846 trend
components (~ 76% variance). The orbital forcing acts at ahierarchically lower strength level (~ 23%), with cycles overlapping a primary trend forcing system.
Nonlinear response of δ18O and SST orbital vs. trend
components

Figure 10 exhibits the exponential relationships between theorbital δ18O and SST components variance vs. the average
δ18O-SST component-1 exponential fit trends by time segments binned at 532 kyr. This plot shows that the exponential increase in signal amplitude towards the Mid-Late
Pleistocene is observed in all δ18O and SST orbital components, from semi-precession to the eccentricity band. The
exponential growth of the orbital response variance isdependent on the long-term mean climate state (LT-MCS)
at a given time (mean ice volume-bottom watertemperature, average sea surface-atmospheric temperature),both at a global scale (LR04) and at a local tropical scale(ODP Site 846). The exponential R2 ranges from 0.7 to 0.9
for the δ18O components and is slightly lower in the SST
components (0.6 to 0.9). The minimum exponential fitting
does not fall below an R2 = 0.64 for SST component-2

Page 16 of 37

subcomp-1-4 and R2 = 0.72 for δ18O component-8-11. The
exponential coefficient ekx for the δ18O semi-precessional
component is 1.05 (comp-8-11) and increases with the orbital period; it is 1.78 for the precessional component(comp-5-7), 1.88 for the obliquity component (comp-3-4),
and 2.93 for the short eccentricity subcomponent (comp-2
subcomp-1-4; 9-15). The exponential coefficient is high(2.37) for the short eccentricity-obliquity modulation subcomponent (detrended comp-1 subcomp-3-11) (Fig. 10a).
These data indicate an increase in δ18O exponential growthfrom semi-precession to ~ 100-kyr eccentricity response towards the Mid-Late Pleistocene. The amplitude flattening of
the obliquity-driven δ18O component-3-4 during the MPT
and post-MPT time shown in Fig. 4 is clearly evident in
Fig. 10a also, a feature related to that observed by Lisiecki
and Raymo (2007). Regarding the tropical SST, the absolute
exponential coefficient ekx for the semi-precessional component is 0.31 (comp-8-11). It is 0.37 for the precessional component (comp-5-7), 0.48 for the obliquity component(comp-3-4), and 0.67 for the short eccentricity subcomponent (comp-2 subcomp-1-4) (Fig. 10b). These results intropical SST are similar to those of the global δ18O, althoughof lesser magnitude. The amplitude contraction of theobliquity-driven SST component-3-4 after the MPT start,
shown in Fig. 8, is very evident even in Fig. 10b.
Figure 11 shows a quantitative analysis of the δ18O
and SST response modulation relative to orbital forcing
as a function of the LT-MCS. The standardized response
sensitivity Rs is given by the following function
Rs ¼ σ 2 resp −σ 2 forc =σ2 forc  100
where σ2resp and σ2forc are the variance of the standardized (0—mean; 1—standard deviation) orbital responsecomponent and forcing calculated by time segments
binned at 532 kyr, respectively. For the purpose of thisexercise, the eccentricity and obliquity response are approximated because small variances in these terms cannot be merged in a composite signal, or a small fractionof different frequency terms are included. Since the variance is small and the aim of this analysis is an order ofmagnitude estimate, the errors should be negligible. The
δ18O short eccentricity response (~ 6.5% variance) was
considered both as a single main subcomponent(comp-2 subcomp-1-4, 9-15; 5.5% variance) containing~ 85% of the short eccentricity response and as theweighted mean by variance of the comp-2 subcomp-1-4,
9-15 (5.5% variance) and the detrend comp-1 subcomp-3-11(2.23% variance), the latter containing also a small fractionof long-term cycles. The δ18O obliquity response (9.9% variance) was calculated as the weighted mean of the comp-3-4(8.1%) and comp-2 subcomp-5-8 (1.8%). Similarly, the SST
short eccentricity response (~ 5.1% variance) wasViaggi Progress in Earth and Planetary Science(2018) 5:81

Page 17 of 37

Table 7 Significance tests of the mean differences among δ18O and SST components and the solution of orbital forcings (Laskar etal. 1993) binned at three arbitrary constant time intervals (532 kyr, 760 kyr, and 1064 kyr). Data standardized on the entire
Plio-Pleistocene record (mean = 0, standard deviation = 1) and, where necessary, inverted to have the same paleoclimatic polarity(positive values for positive forcing and warming)
Component/forcing

Time bin (kyr)Meana

Sum of squares

df

F

Sig.

LR04 comp-1

532− 1.47 ÷ 1.18

5059.52

911,504.34

0.00

760

− 1.51 ÷ 1.13

4996.27

6

13,708.49

0.00

1064

− 1.43 ÷ 1.14

4917.43

4

16,271.23

0.00

LR04 comp-1 subcomp-1

LR04 comp-1 exp.fit. trend

LR04 comp-2-11 (orbitals)

SST comp-1

SST comp-1 subcomp-1

SST comp-1 exp.fit. trend

SST comp-2-11 (orbitals)

Eccentricity

Obliquity

Precession

65°N Insolation at June 21532− 1.47 ÷ 1.21

5257.57

950,499.31

0.00

760

− 1.50 ÷ 1.17

5182.42

6

33,599.07

0.00

1064

− 1.44 ÷ 1.17

5116.42

4

33,559.24

0.00

532

− 1.72 ÷ 1.40

5262.71955,162.47

0.00

760

− 1.62 ÷ 1.35

5204.25

6

40,161.50

0.00

1064

− 1.49 ÷ 1.28

5094.35

4

30,132.42

0.00

532

− 0.08 ÷ 0.09

8.80

9

0.98

0.46

760

− 0.01 ÷ 0.02

0.74

6

0.12

0.99

1064

− 0.01 ÷ 0.01

0.33

4

0.08

0.99

532

− 1.45 ÷ 1.26

4644.7496235.990.00

760− 1.42 ÷ 1.27

4468.8766339.54

0.00

1064

− 1.38 ÷ 1.16

4389.5948244.210.00

532− 1.49 ÷ 1.26

4950.98

924,819.72

0.00

760

− 1.49 ÷ 1.33

4838.36

6

18,152.80

0.00

1064

− 1.42 ÷ 1.18

4722.96

4

17,566.50

0.00

532

− 1.67 ÷ 1.50

5006.76

949,988.56

0.00

760

− 1.57 ÷ 1.44

4949.41

6

36,725.39

0.00

1064

− 1.44 ÷ 1.36

4842.284

27,747.29

0.00

532

− 0.09 ÷ 0.04

6.25

9

0.69

0.71

760

− 0.01 ÷ 0.02

0.89

6

0.15

0.99

1064

− 0.02 ÷ 0.01

0.60

4

0.15

0.96

532

− 0.30 ÷ 0.25

169.60

9

19.43

0.00

760

− 0.28 ÷ 0.28

143.98

6

24.64

0.00

1064

− 0.22 ÷ 0.14

77.44

4

19.63

0.00

532

− 0.12 ÷ 0.11

35.42

9

3.96

0.00

760

− 0.12 ÷ 0.12

37.02

6

6.21

0.00

1064

− 0.11 ÷ 0.10

31.96

4

8.03

0.00

532

− 0.02 ÷ 0.01

0.30

9

0.03

1.00

760

0.00

0.02

6

0.00

1.00

1064

− 0.01 ÷ 0.01

0.12

4

0.03

1.00

532

− 0.08 ÷ 0.06

8.69

9

0.97

0.47

760

− 0.05 ÷ 0.06

8.89

6

1.48

0.18

1064

− 0.05 ÷ 0.04

6.37

4

1.59

0.17

df degrees of freedom, F Fisher test, Sig statistical significanceaMin-max of the mean among time segments

approximated as a single comp-2 subcomp-1-4 (3.65% variance) explaining ~ 72% of the short eccentricity response
and as weighted mean of the comp-2 subcomp-1-4 (3.65%
variance) and the detrended comp-1 subcomp-3-11 (2.84%

variance), the latter containing a small fraction of long-term
cycles. The weighted mean of the SST comp-3-4 (4.5% variance) and comp-2 subcomp-9-10 (0.27% variance) was calculated as the SST obliquity response, covering ~ 87% of theViaggi Progress in Earth and Planetary Science(2018) 5:81

Page 18 of 37

Fig. 9 Plots of the 95% confidence interval for the mean of a δ18O components and c SST components compared to orbital forcings (Laskar etal. 1993), binned at 532-kyr constant time intervals. Overlapped histograms of b Plio-Pleistocene δ18O and d SST orbital component-2-11 vs. trend
component-1 subcomp-1. The averages of the astronomical components oscillate around the standardized mean with little statisticallysignificance differences and a roughly normal distribution. Instead, the averages of the trend components differ very significantly, deviating
markedly from the mean with a general multimodal distribution. Data standardized on the entire Plio-Pleistocene record (mean = 0, standard
deviation = 1) and, where necessary, inverted to have the same paleoclimatic polarity (positive values for positive forcing and warming)

total signal (4.77% of 5.5%). A self-sustained climate system,
which is only paced by orbital cycles, would produce a response sensitivity near zero because the variance of the
response would match the variance in forcing, whereas
values much larger than zero would indicate a nonlinear
reinforcement of the signal response. Negative Rs values
show a response variance lower than orbital forcing, suggesting a damping mechanism for the climate system.
Figure 11a exhibits a δ18O nonlinear eccentricity Rs up to
400% towards the Mid-Late Pleistocene similar to that of aprecession response in shape but with a magnitude up to
200%. The two eccentricity curves are very similar, indicating a negligible error in the approximations. The δ18O obliquity Rs shows a nonlinear amplification up to 180%
towards the end of the MPT, followed by strong variancedepletion near to zero (balance line) post-MPT. This patternof low δ18O obliquity variance is similar to thatshown in Fig. 5b from Lisiecki and Raymo (2007). Remarkably, the asymptotic decay of the δ18O Rs towards the Early Pliocene becomes all negative up to− 100%, suggesting a damped orbital response during
the time of minimum average ice volume and higher
temperature. This Rs pattern is similar considering
the SST (Fig. 11b), although more scattered and witha smaller positive magnitude which does not reach
200%. The SST obliquity Rs also shows strong variance depletion up to − 20% post-MPT as well as the
new feature of asymptotic decay up to − 100% of all
the SST Rs towards the Early Pliocene.
Link between the pCO2 and the SST-δ18O long-term trend
components

Figure 12 shows the Plio-Pleistocene cross plots amongthe δ18O, SST, and pCO2 trend and orbital components.

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Page 19 of 37

Fig. 10 Exponential relationships of the orbital variance of a δ18O components and b SST components vs. the average of component-1 exponential
fit trends by time segments binned at 532 kyr. The right y-axis shows the variance of comp-2-11 (orbitals) with a different scale to not compress theleft y-axis values

These plots offer an overall view of the long-term trend
links among the proxies. Although it is possible to express the trend relationship as exponential regressions, it
was decided to use a simpler linear fitting since the correlation coefficients are similar, in order to quantify thetrend rates of change. The δ18O trend components exhibit agood linear interpolation vs. the boron-based pCO2 trend of
the Sites 668B-999A (Fig. 12a) and the alkenone-based
pCO2 trend of the Site 999A (Fig. 12b); differently, the δ18O
orbital components are unrelated to the same pCO2
long-term trends. Similar relationships are shown by theδ18O trend components and the alkenone-based SST trendsof the Site 846 (Fig. 12c) and Site 999A (Fig. 12d). The tropical SST trend components of the Site 846 are linearly related to the boron-based atmospheric pCO2 trend of theSites 668B-999A (Fig. 12e) and the alkenone-based pCO2
trend of the Site 999A (Fig. 12f); again, the SST orbital components are unrelated to the same pCO2 long-term trends.
However, in the plots with more data points, the amplitude increase of the orbital components is clearly visible as pCO2 and SST decrease, corroborating theanalysis of Figs. 10 and 11. These results show that theViaggi Progress in Earth and Planetary Science(2018) 5:81

Page 20 of 37

Fig. 11 Analysis of the order of magnitude in a δ18O and b SST response modulation relative to orbital forcing (standardized response sensitivity Rs,
see text) as a function of the LT-MCS. The balance line (0%) indicates that the variance of response matches the variance in forcing, suggesting a selfsustained climate system, which is only paced by orbital cycles. Values much larger than zero would indicate a nonlinear reinforcement of the signal
response. Negative Rs shows a response variance lower than orbital forcing, suggesting a damping mechanism for the climate system

climate sensitivity to orbital forcing is lower when CO2
and SST are higher (Early Pliocene) and that it increases
non-linearly when CO2 and SST are low (Mid-Late Pleistocene), in agreement to Martinez-Boti et al. (2015), afeature linked to the LT-MCS as shown in the present
study. Despite that Site 846 is considered not “typical”
tropical SST record, Fig. 13a shows a very close link

between its trend component and that of the global
ΔSST stack (Martinez-Boti et al. 2015) identified by SSA
and driving 76.2% of the variance (ΔSST comp-1), noteworthy the same magnitude of the Site 846 SST and global δ18O trend components (76.6% and 76.4%,
respectively). The inclusion of the Site 846 among theten of the ΔSST stack should not significantly affect theViaggi Progress in Earth and Planetary Science(2018) 5:81

global character of the signal. Moreover, Lisiecki (2010b)
calculates 0.28 °C/Myr the upwelling cooling effect at
Site 846 which roughly matches a 30% fraction of thetotal 1.0 °C/Myr trend (Herbert et al. 2010a), resulting
about 70% of the Site 846 cooling trend to be due tonon-regional effects, according to the results of the SSA.
Figure 13 shows similar links of Fig. 12 among the global ΔSST long-term trend and the pCO2 trends at Sites
668B and 999A (Fig. 13c) and 999A (Fig. 13d) and theglobal δ18O (Fig. 13b). Figure 13e, f exhibits the same relationships by converting the standardized ΔSST into
SST (°C) using the function of Fig. 13a. The SST vs.
boron-based CO2 trend slopes correspond to a tropicalrate of change of 25.6 μatm/°C (Fig. 12e), very close to
the global converted ΔSST stack of 25.5 μatm/°C(Fig. 13e). At 25 °C and atmospheric pressure, thesevalues are equivalent to a rate of change of 26.4 ppm/°C
and 26.3 ppm/°C, respectively. Coherently, as shown in
Fig. 14b, the two SST trends are subparallel in time, although the Site 846 is cool-shifted. Considering thealkenone-based CO2 trend slopes, these rates of changes
are of 42.6 ppm/°C (Fig. 12f ) and 39.2 ppm/°C (Fig. 13f ),respectively. Such variability in rate of changes (~ 26 ppm/°C
vs. ~ 41 ppm/°C) could be due to differences in the CO2proxy estimation (boron-based vs. alkenone-based). These
slope data offer an overall view of the Plio-Pleistocene
long-term climate sensitivity, and the oscillations around thelinear trend should indicate that it can change for discrete
time intervals. On the other hand, the climate sensitivity toorbital forcing change by time in response to the initial conditions set by LT-MCS (Tachikawa et al. 2014; Martinez-Boti
et al. 2015; Lo et al. 2017; this study).
These results authorize to be considered of global innature the non-Milankovitch trend links among pCO2,
SST, and δ18O. Considering that the main GHGs, suchas carbon dioxide and methane, operate physically as themost effective positive (warming) climate forcings (Hansen
et al. 2005; Hansen and Sato 2012; Van Der Meer et al.
2014), the hypothesis that the global long-term variation of
atmospheric GHG concentration plays a relevant role inexplaining the origin of the δ18O and SST trend components, driving ~ 76% of the Plio-Pleistocene signal variance,
appears plausible.
Figure 14 shows the correlation over time between theδ18O enrichment trend components (LR04 comp-1
subcomp-1 and LR04 comp-1 exponential fit) and theCO2 reduction trend of the ODP Sites 668B and 999A
(Fig. 14a) and the SST cooling trends of the global ΔSST
stack, ODP Sites 846 and 999A (Fig. 14b).
During the Zanclean, the pCO2 alkenone ranged from
450 to 330 ppm, the highest values of the entire
Plio-Pleistocene; the tropical SST was higher than the
Plio-Pleistocene average; and the global δ18O was depleted more than the average (Fig. 14). Conversely, thePage 21 of 37

mean orbital forcing was slightly negative (Fig. 9). Theseobservations indicate Zanclean as a period of a higher
than average positive climate state, supporting the hypothesis that the Zanclean warmer climate may havebeen driven by an important greenhouse effect due tohigher atmospheric CO2 (Bartoli et al. 2011). The lack ofsignificant Northern Hemisphere ice growth and 100-kyr
climate response prior to 3.3 Myr suggest that pCO2
alkenone above 320–330 ppm is associated with warming sufficient to prevent significant Northern Hemisphere ice sheet growth even during cold orbital phases(Seki et al. 2010). The pCO2 and SST show a significant
reduction during the Zanclean (subtrend-I), supporting
the hypothesis that the ONHG was triggered by the decreasing trend in atmospheric CO2 (Pagani et al. 2010;
Seki et al. 2010). Subtrend-II (Piacenzian) is also consistent with a reduction of pCO2 alkenone and SST, suggesting an INHG triggered by a further decrease inatmospheric CO2 according to Lunt et al. (2008), whostated that Greenland glaciation was primarily controlled
by a decrease in atmospheric CO2 during the Late
Pliocene. Nevertheless, the Piacenzian remains a periodof a δ18O-SST higher than average positive climate state,with slightly negative or neutral mean orbital forcing(Fig. 9). For this reason, the Piacenzian, excluding the
ONHG event, is still characterized by an average warming fingerprint due to a residual greenhouse effect.
According to this hypothesis, Lunt et al. (2012), using amodeling approach, showed that on a global scale, thelargest contributor to mid-Piacenzian warmth (3.26–
3.03 Myr) is high CO2. Subtrend-III (Gelasian-Calabrianp.p.) does not appear to comply with the pCO2 alkenone
data, probably due to the low-resolution proxy in thisrange, while it is more consistent with the decreasingboron-based pCO2 and the SST subtrends. For the firsttime during the Plio-Pleistocene, the long-term
pCO2and SST are below or close the mean and the δ18O
is enriched more than the average (Figs. 14 and 9). The
mean orbital forcing is neutral or slightly positive. The net
result of these two unbalanced forcings may be a significant cooling effect and the start of an average “icehouse”
world. Subtrend-IV (Calabrian-Middle Pleistocene p.p.) is
consistent with a further reduction of pCO2, reaching
peak concentrations of 240–250 ppm (MIS-16), and withthe continuation of the SST cooling trend. The long-term
pCO2, SST, and δ18O are consolidated below the mean,
while the mean orbital forcing appears slightly positive,
resulting in a lower than average negative climate state,
strengthening the icehouse effect. No CO2 variation corresponds to the MPT, although it falls at the beginning ofsubtrend-IV after a long time interval of surface/deep-water cooling and ice volume growth. Finally, the MBO(Mid-Late Pleistocene) marks a mild recovery of thelong-term pCO2 clearly visible in the EPICA trend comp-1

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Page 22 of 37

Fig. 12 Plio-Pleistocene cross plots between selected pCO2 proxies and the δ18O (a, b) and SST (e, f) SSA-components. Plio-Pleistocene cross plots
between selected SST proxies and the δ18O (c, d) SSA-components. A clear linear interpolation is shown between long-term trend components,
whereas orbital components-2-11 are unrelated to the long-term trends. Proxy data: composite atmospheric boron-based pCO2 proxy from ODP Sites
668B-999A (Hönisch et al. 2009; Bartoli et al. 2011), atmospheric alkenone-based pCO2 and SST proxies from ODP Site 999A (Seki et al. 2010), alkenonebased SST proxy from ODP Site 846 (Herbert et al. 2010b), and LR04 δ18O stack (Lisiecki and Raymo 2005b)

subcomp-3 and less distinctly in the CO2 proxy records (seeAdditional file 1) (Lisiecki 2010b; Seki et al. 2010) and in the
SST proxies (Fig. 14). The strong isotopic enrichment related to MIS-16 at the end of the subtrend-IV correlates wellwith the absolute CO2 minimum in the EPICA comp-1
subcomp-3 at ~ 0.65 Myr. Considering that the δ18O response lags CO2 (Shackleton 2000; this study), thecorrelation between the CO2 maximum in the EPICA
comp-1 subcomp-3 at ~ 0.4 Myr and the related but slightly
delayed isotopic depletion in the δ18O comp-1 subcomp-1 isnoteworthy (MBO). This observation supports the hypothesis that MIS-11, one of the most intense and prolonged
interglacial oscillations of the Middle-Late Pleistocene,
which is not proportional to eccentricity (Imbrie et al. 1993)

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Page 23 of 37

Fig. 13 Late Pliocene-Pleistocene cross plot between trend components of global ΔSST stack and tropical SST at Site 846 (a). Cross plots between the
ΔSST stack and b δ18O, c boron-based, and d alkenone-based pCO2 SSA-components. Cross plots between the converted ΔSST stack (°C by means of
the function in a) and e boron-based pCO2 and f alkenone-based pCO2 trends. The long-term trend components are clearly linked, whereas orbital
components are unrelated to the long-term trends. Proxy data: composite atmospheric boron-based pCO2 proxy from ODP Sites 668B-999A (Hönischet al. 2009; Bartoli et al. 2011), atmospheric alkenone-based pCO2 and SST proxies from ODP Site 999A (Seki et al. 2010), alkenone-based SST proxyfrom ODP Site 846 (Herbert et al. 2010b), LR04 δ18O stack (Lisiecki and Raymo 2005b), and global ΔSST stack (Martinez-Boti et al. 2015)

or cool insolation-induced interglacials (Berger et al. 2012),
may be an astronomical event boosted by GHG mid-term
recovery, which reached the highest value at ~ 400 kyr(Luthi et al. 2008; Berger and Yin 2012; Berger et al. 2016).
The origin of this GHG mid-term recovery could be partly
due to the long-term orbital cycles that reach a maximum
relative (warming) at ~ 300–200 kyr (Figs. 2 and 6); the400-kyr eccentricity cycle really reaches a maximum atabout 200 kyr (Laskar et al. 2011). However, the timing ofthese astronomical events does not seem consistent with theCO2 maximum of EPICA at ~ 400 kyr. In addition, we cannot exclude a contribution of the circum-Pacific explosive
activity (Jicha et al. 2009). The Late Pleistocene was a period
of quiescence of the circum-Pacific explosive activity,

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Page 24 of 37

Fig. 14 Comparison over time of the long-term trend components of δ18O and the a pCO2 and b SST proxies. Composite atmospheric boron-based
pCO2 proxy from ODP Sites 668B-999A (Hönisch et al. 2009; Bartoli et al. 2011), atmospheric alkenone-based pCO2 and SST proxies from ODP Site 999A
(Seki et al. 2010), alkenone-based SST proxy from ODP Site 846 (Herbert et al. 2010b), EPICA pCO2 record (Berger et al. 2016), δ13C-based pCO2 proxy(Lisiecki 2010b), LR04 δ18O stack (Lisiecki and Raymo 2005b), boron-based pCO2 proxy from ODP Site 999 and global ΔSST stack (Martinez-Boti et al. 2015).
The horizontal lines are the Plio-Pleistocene average of pCO2 alkenone (ODP Site 999A) and the standardized mean of the SST. All records are plotted on
their own age scale

which perhaps contributed to the mid-term CO2 increment by a reduction in volcanic aerosols and theconsequent minor algal ash-fertilization of the Pacific
Ocean waters (Jicha et al. 2009).
Lagged patterns of δ18O and SST orbital components

The phase relationships between the standardized La93
orbital forcings and the orbitally driven LR04-SST components are shown in Table 8 for three time intervals. The

results indicate that the δ18O and SST orbital components
are all delayed on average compared to the related forcingof eccentricity, obliquity, and precession (with the small
significant exception of the eccentricity in the pre-ONHG
time due to very low SST signal coherency), although thepatterns of coherency and phase lag are quite different
across the time intervals.
These cross-spectral results are in general agreementwith Imbrie et al. (1993), Shackleton (2000), Lisiecki andViaggi Progress in Earth and Planetary Science(2018) 5:81

Raymo (2005a), and Tziperman et al. (2006) and indicatethat the climate system takes a certain amount of timeto transfer the effect of orbital forcing to the global ice
volume/deep water temperature and to the sea surface
temperature. The coherency is generally lower in thetropical SST signal probably due to differences in theforcings controlled by high-latitude and low-latitude signals. The lagging of the δ18O signal is usually larger thanthat of the SST (Table 8) because the ice volume component lags the deep water temperature signal in the δ18O
record (Lisiecki and Raymo 2005a; Herbert et al. 2010a).
The coherency of the eccentricity signal increases frompre-ONHG time (incoherent) towards the Mid-Late

Page 25 of 37

Pleistocene, reflecting the progressive increase of eccentricity power both in global δ18O and local tropical SST
signals. The delta-lag in eccentricity between δ18O andSST reaches its maximum during the Late Pleistocene(4.31–0.44 kyr), suggesting an increased contribution ofthe δ18O ice volume component over the watertemperature signal. The δ18O and SST obliquity signals
have the largest coherency and exhibit an average increasing trend in the lags towards the Mid-Late Pleistocene, starting from 4.10–3.46 kyr to 5.78–5.28 kyr and,
finally, reaching values of 7.42–2.71 kyr (Table 8), in
general agreement with Lisiecki and Raymo (2005a). The
δ18O-SST delta-lag in obliquity is relatively small until

Table 8 Cross-spectral analysis results among the orbital forcings (Laskar et al. 1993) and the orbitally driven δ18O and SST
components for three time intervals. Data standardized (mean = 0, standard deviation = 1) and, where necessary, inverted to havethe same paleoclimatic polarity (positive values for positive forcing and warming)
Time interval

Orbital forcing

Response
component

Cross-spectrum
freq. (kyr−1)

Cross-spectrum
period (kyr)Coherency

Phase shift

Post-MPT start

Eccentricity

LR04 comp-2
subcomp-1-4, 9-15

0.01074293

0.72− 16.65°− 4.31 kyr

δ18O component
lags eccentricity

SST comp-2
subcomp-1-4

0.01074293

0.70− 01.70°− 0.44 kyr

SST component
lags eccentricity

LR04 comp-3-4

0.024414

410.93− 65.20°− 7.42 kyr

δ18O component
lags obliquity

SST comp-3-4

0.024414

410.68− 23.78°− 2.71 kyr

SST component
lags obliquity

LR04 comp-5-7

0.04199224

0.78− 74.50°− 4.93 kyr

δ18O component
lags precession

SST comp-5-7

0.04199224

0.72− 53.24°− 3.52 kyr

SST component
lags precession

LR04 comp-2
subcomp-1-4, 9-15

0.01074293

0.47− 27.57°− 7.13 kyr

δ18O component
lags eccentricity

SST comp-2
subcomp-1-4

0.01074293

0.45− 21.37°− 5.53 kyr

SST component
lags eccentricity

LR04 comp-3-4

0.024414

410.98− 50.79°− 5.78 kyr

δ18O component
lags obliquity

SST comp-3-4

0.02465841

0.95− 46.85°− 5.28 kyr

SST component
lags obliquity

LR04 comp-5-7

0.044678

22

0.89− 57.37°− 3.57 kyr

δ18O component
lags precession

SST comp-5-7

0.042236

24

0.28− 119.35− 7.85 kyr

SST component
lags precession

LR04 comp-2
subcomp-1-4, 9-15

0.009766

102

0.32

− 01.48°− 0.42 kyr

δ18O component
lags eccentricity

SST comp-2
subcomp-1-4

0.01074293

0.13+ 07.06°

1.83 kyr

SST component
leads eccentricity

LR04 comp-3-4

0.024414

410.95− 36.02°− 4.10 kyr

δ18O component
lags obliquity

SST comp-3-4

0.024414

410.81− 30.42°− 3.46 kyr

SST component
lags obliquity

LR04 comp-5-7

0.042480

24

0.42− 88.11°− 5.76 kyr

δ18O component
lags precession

SST comp-5-7

0.042480

24

0.59− 87.60°− 5.73 kyr

SST component
lags precession

Obliquity

Precession

ONHG-MPT start

Eccentricity

Obliquity

Precession

Pre-ONHG

Eccentricity

Obliquity

Precession

Viaggi Progress in Earth and Planetary Science(2018) 5:81

the beginning of the MPT and becomes large after thistransition, probably following global ice volume growth.
The cross-spectral results relative to precession are morevariable. The precession lag is basically constant (5.76–
5.73 kyr) during pre-ONHG time, then it becomes quite“unstable” (3.57–7.85 kyr) between ONHG and MPT
and reaches lags of 4.93–3.52 kyr post-MPT. The minimum δ18O phase lag relative to precession (3.57 kyr)during the ONHG-MPT time interval suggests the influence of a large global deep water temperature signal(Lisiecki and Raymo 2005a). During the same time interval, the precession’s SST lag reached its maximum(7.85 kyr) which could be due to the effect of a deep
water upwelling in the EEP (ODP Site 846). Table 9
exhibits, for the last 5 Myr, the cross-spectral analysis
results between SST and LR04 orbital components,
showing a slight δ18O lagged response vs. SST, in agreementwith Herbert et al. (2010a). The slightly earlier response ofthe δ18O eccentricity component (0.05 kyr) is due to the
pre-ONHG weak and unstable eccentricity signal. After the
ONHG, the eccentricity response is strengthened and becomes more stable resulting in a slight δ18O lag.

Discussion
Origin of the δ18O and SST trend components

The LR04 trend component is related to the Plio-Pleistocene
interval of the Cenozoic δ18O enrichment trend starting at
the end of the Early Eocene Climatic Optimum (EECO,
53–51 Myr), documenting the progressive developmentof the Antarctic and Boreal ice sheets (Zachos et al.
2001). The gradual climate cooling during the Cenozoic
is generally correlated to a CO2 concentration decreasein the atmosphere (Raymo and Ruddiman 1992; Zachos
et al. 2001; Royer 2006; Zachos et al. 2008; Van de Wal
et al. 2011; Zhang et al. 2013; Van Der Meer et al.
2014). Recent studies on Pliocene and Pleistocene alkenone δ13C and boron-based pCO2 proxies support thishypothesis (Seki et al. 2010; Pagani et al. 2010; Bartoli

Page 26 of 37

et al. 2011). Stap et al. (2014) find that the atmospherictemperature during the past 800 kyr has been controlled by a complex interaction of CO2 and insolation,and both variables serve as thresholds for the NHG.
Using CO2 proxies and geochemical carbon cycle
models, Royer (2006) shows that atmospheric CO2 exhibits a pervasive and tight positive correlation with
globally averaged surface temperatures, indicating that
CO2 played a strong role in driving the global climate
for much of the Phanerozoic. The Cenozoic climaticrecord provides an important case study to understandthe relationships between carbon cycling and climate(Zachos et al. 2001, 2008). In fact, much of the early
Cenozoic was characterized by noticeably higher concentrations of GHGs and a much warmer mean global
temperature. The most extreme event being between 53
and 51 Myr (EECO), when CO2 concentration was on
the order of 1000–1500 ppm, and the global temperatures reached a long-term Cenozoic maximum (Zachos
et al. 2001, 2008; Beerling and Royer 2011). Zhang et al.
(2013) estimate the new pCO2 data for the Middle
Miocene to be notably higher than published records
with average pCO2 concentrations in the range of 400–
500 ppm, suggesting that CO2 levels were highest
during a period of global warmth associated with the
Middle Miocene Climatic Optimum (MMCO, 17–
14 Myr). Instead, during the mid-Pleistocene ice ages,
CO2 fell to 170–190 ppm, and the global temperaturereached a long-term Cenozoic minimum (Beerling andRoyer 2011; Berger et al. 2016). Considering a CO2
variation from approximately 1000 ppm in the early
Cenozoic to 170 ppm in the icy Pleistocene, the resulting negative climate forcing computed by Hansen andSato (2012) for this CO2 range exceeded 10 W/m2,
which is consistent with the observed strong Cenozoic
cooling. These studies confirm the link among global
temperature, the cryosphere, and atmospheric CO2.
The present study suggests a link between the

Table 9 Cross-spectral analysis results between orbitally driven δ18O and SST components for the last 5 Myr. Data standardized(mean = 0, standard deviation = 1) and, where necessary, inverted to have the same paleoclimatic polarity (positive values for
positive forcing and warming)
SST component

δ18O component

Cross-spectrum
freq. (kyr−1)

Cross-spectrum
period (kyr)Coherency

Phase shift

SST comp-2
subcomp-1-4

LR04 comp-2
subcomp-1-4, 9-15

0.010498

95

0.89

0.20°

0.05 kyrδ18O leads SST eccentricity

SST comp-2
subcomp-1-4(post-ONHG)

LR04 comp-2
subcomp-1-4, 9-15
(post-ONHG)

0.01049895

0.84− 0.38°− 0.10 kyr

δ18O lags SST eccentricity

SST comp-3-4

LR04 comp-3-4

0.024414

410.92− 24.57°− 2.80 kyr

δ18O lags SST obliquitycomp.

SST comp-5-7

LR04 comp-5-7

0.044800

220.46− 16.40°− 1.02 kyr

δ18O lags SST precession
comp.

comp.

comp.

Viaggi Progress in Earth and Planetary Science(2018) 5:81

atmospheric pCO2 long-term variations and the
SST-δ18O trend components and highlights the primaryrole of the long-term change of the atmospheric GHG concentration in explaining the origin of the high-variance,
highly significant δ18O and SST trend components.
Factors controlling long-term atmospheric composition

Although it is generally accepted that the progressive
Cenozoic cooling of the global climate resulted from along-term decrease in atmospheric CO2, details on how
carbon cycle affects CO2 and a clearer distinction of theinvolved time scales have not been fully resolved. Because of the relatively small size of the atmospheric reservoir, minimal imbalances between CO2 sources andsinks cause rapid fluctuations in atmospheric CO2 levels
and global mean temperatures (Brady and Gislason
1997; Kent and Muttoni 2013). The results of thepresent study suggest ~ 76% of the Plio-Pleistocene climate variance could be driven by non-orbital carbon
cycle. On time scales longer than one million years,
atmospheric CO2 levels are primarily controlled by thebalance between the rate of volcanic and metamorphic
sources from the Earth’s interior and the rate of sinkthrough the chemical weathering of silicate minerals andthrough the burial of organic carbon (Walker and Hays
1981; Berner et al. 1983; Raymo and Ruddiman 1992;
Brady and Gislason 1997; Schrag et al. 2002; Hansen andSato 2012). Plate tectonics may affect the primary inputs
of atmospheric CO2 through mantle outgassing at seafloor spreading ridges, volcanoes, and the metamorphism of carbonate rocks along subduction zones (Van Der
Meer et al. 2014). Long-term sinks for CO2 are controlled by the burial of organic matter and by the weathering of silicate rocks to form carbonates (Raymo andRuddiman 1992; Berner and Kothavala 2001; Schrag etal. 2002; Van Der Meer et al. 2014). The weathering ofsilicate rocks is a chemical process whereby CO2 consumption is buffered by the temperature dependence on
the weathering rate, preventing the so-called “ice catastrophe” when temperatures fall (Walker and Hays 1981).
As the atmospheric CO2 concentration rises, thetemperature and precipitation increase. They thereby
enhance chemical weathering and atmospheric CO2
depletion, preventing “greenhouse” disasters on geological timescales. Chemical weathering is also affectedby changes in continental land area and paleoaltitude
distribution (Berner et al. 1983; Schrag et al. 2002). The
alteration of surficial marine basalts at low temperatures(< 40 °C) is another potentially important sink for atmospheric CO2 over geologic time (Brady and Gislason
1997). The numerical modeling of the carbon cycle byBerner et al. (1983) suggests that changes in volcanic
outgassing over the last 100 million years were theprimary cause of the transition between the greenhouse

Page 27 of 37

climates of the Cretaceous and the ice ages of the
Plio-Pleistocene. Their results indicate that the CO2 content in the atmosphere is highly sensitive to changes inthe seafloor spreading rate and continental land area. A
tropical concentration of continental area may lead to amore efficient burial of organic carbon through increased tropical river discharge and the chemical weathering of silicate rocks (Schrag et al. 2002). By affecting
CO2 volcanic degassing, metamorphic-magmatic decarbonation and changes in the continental land area
available for chemical weathering, plate tectonics is considered the main long-term controlling factor of the
global climate (Berner et al. 1983; Van Der Meer et al.
2014). However, it is a matter of debate as to whichmain mechanism can adjust the balance between theCO2 source through geotectonic volcanic degassing andthe CO2 sink by the chemical weathering of silicate
rocks. Some studies show conflicting results regarding
the role of the seafloor spreading rate, which is believed
to be the major control of CO2 degassing (Berner et al.
1983; Raymo and Ruddiman 1992). In fact, according toRowley (2002), the seafloor spreading rate remained
relatively constant over the past 180 Myr, while Cognéand Humler (2006) indicate a general increasing trend in
the last 50 Myr. These results are discordant with theobserved cooling trend of the Cenozoic climate. Recently, a study of carbon-cycle evolution indicated that
increasing physical erosion in response to widespread
mountain building is the most prominent driver of carbon cycling over the past 100 Myr (Li and Elderfield2013). Raymo and Ruddiman (1992) proposed that the
uplift of the Tibetan plateau was an important driving
force behind the Cenozoic long-term cooling by atmospheric CO2 consumption due to the increase of thechemical weathering of silicate rocks. However, theCenozoic long-term shift in Earth’s climatic state could
be really related to regional transient differences involcanic emissions by plate tectonic activity, which wereparticularly high during, e.g., the Late Paleocene and
Early Eocene epochs (North Atlantic rifting and India-Asian
northward migration) (Raymo and Ruddiman 1992; Zachos
et al. 2001; Mosar et al. 2002; Zachos et al. 2008; Le Breton
et al. 2012). From this viewpoint, a recent study on the
Himalayan-Tibetan orogenic system shows a multiple collision model in which a continent-continent collision beganaround the Eocene/Oligocene boundary (34 Myr) (Aitchison
et al. 2007), promoting the Early Oligocene glaciation(Zachos et al. 2001). Corroborating this hypothesis, theNorth Atlantic seafloor spreading rate (NASSR) exhibits aclear decreasing trend from high values in the Early Eocene(55–50 Myr) to a very low rate for most of the Oligocene(Mosar et al. 2002; Le Breton et al. 2012). On the other
hand, recent studies seem to converge towards a significant
250–200-Myr period variation of the seafloor spreading rate

Viaggi Progress in Earth and Planetary Science(2018) 5:81

(Becker et al. 2009; Coltice et al. 2013). Becker et al. (2009)
proposed a spreading rate long-term slowdown since
140 Myr with minor cycles superimposed over time, such asa relative increase at ~ 50 Myr and a minimum at ~ 32 Myr.
Even the timing of the global average spreading rate (GASR)
calculated by Conrad and Lithgow-Bertelloni (2007) appears
to support the hypothesis of a direct geotectonic “forcing” tothe long-term composition of atmospheric GHGs. On a
more detailed timescale, a study in the North Atlantic shows
the consistency of the NASSR relative maximum at~ 15 Myr and the MMCO (Merkouriev and DeMets
2014). The correlation of the NASSR minimum at approximately 1.2–1.0 Myr (the lowest value in the last
20 Myr) and the CO2 proxy minimum at 0.9–0.7 Myr(Fig. 14a), which correlate to the strong LT-LR04 coolingat the end of subtrend-IV, seems plausible. Indeed, thecontribution of the important negative forcing created byaerosols (Hansen et al. 2005) emitted into the atmosphere
during phases of explosive volcanic paroxysm in thecircum-Pacific (Jicha et al. 2009) must also be taken in account. In fact, hundreds of explosive eruptions occurred
during an extremely vigorous period of circum-Pacific volcanism beginning in the Late Eocene, which likely led tothe production of sulfur aerosols in the stratosphere andnet cooling. The feedback mechanism related to the enormous supply of ash that fertilized the Pacific Ocean waters
might have caused algal growth and CO2 uptake andthereby contributed to Oligocene glaciation at approximately 34 Myr (Jicha et al. 2009). It is interesting to observe how the Plio-Pleistocene in the circum-Pacific was aperiod of intense explosive volcanic activity characterized
by three activity peaks at approximately 4.8, 2.8, and
1.5 Myr (Jicha et al. 2009). We cannot exclude that these
events contributed in part to the three most pronounced
segments of δ18O enrichment highlighted in the LT-LR04
record (subtrend-I, subtrend-II, and subtrend-IV), culminating in the ONHG, INHG, and MIS-16 events, respectively. Moreover, the Late Pleistocene was a period ofquiescence of the circum-Pacific explosive activity, which
may have contributed to the mid-term CO2 incrementand δ18O impoverishment that occurred at approximately
400–300 kyr, a possible cause of the MBO. A further contribution to the global emission of CH4 and CO2 may alsocome from the thermogenic cracking of organiccarbon-rich sediments in large igneous provinces (LIPs),
which are associated with geotectonic spreading areas or
subduction zones (Svensen et al. 2004). Finally, anotherclimate-controlling factor is oceanic paleocirculation. The“Panama Hypothesis” states that the gradual closure of the
Panama seaway, between 13 and 2.6 Myr led to a decreased mixing of the Atlantic and Pacific water masses
and a strengthening of the Atlantic thermohaline circulation, increased temperatures and precipitation in Northern
Hemisphere high latitudes, and culminated with thePage 28 of 37

INHG during the Pliocene, at approximately 3.2–2.7 Myr(Bartoli et al. 2005; Lunt et al. 2008). However, Lunt et al.
(2008), using a more sophisticated general circulation
model (GCM), demonstrated that the ice volume difference between the “closed” and “open” seaway model configurations is small and concluded that it was not a major
forcing mechanism. This view accommodates δ18O andSST components-8-11 of inter-tropical semi-precessional
origin related to the oceanic heat transport towards highlatitudes (Hagelberg et al. 1994; Berger et al. 2006), which
contributes very little to the total variance of thepaleoclimatic response (~ 0.6–0.5%). In summary, thegeotectonic forcing appears a plausible hypothesis toexplain the long-term changes in atmosphere composition(GHGs and volcanic aerosols) through direct (CO2 outgassing and explosive volcanism) and indirect (orography anderosion, paleogeographic configuration, oceanic paleocirculation, and ocean fertilization) control by plate tectonics.
Nonlinear response and mean climate state

An interesting result of the present study shows that thenonlinear increase in signal amplitude from the Early
Pliocene to the Late Pleistocene is observed in all δ18O
and SST astronomical components and is exponentiallylinked to the LT-MCS (Figs. 10 and 11). These results
reinforce the notion that the climate sensitivity of orbital
forcings may be dependent on the mean climate state inagreement with Lisiecki and Raymo (2007) and Hansen
et al. (2013), with possible threshold effects typical of
nonlinear systems (Rial et al. 2004). Because the
Plio-Pleistocene orbital forcings do not exhibit a similar
increase in amplitude, it is believed that the process controlling the nonlinear behavior of the signal response
originates inside the climate system. The net balance of
positive and negative feedback mechanisms is the most
likely process behind most of the nonlinearities in the
climate (Berger and Loutre 1997; Archer et al. 2004;
Rial et al. 2004; Lisiecki and Raymo 2007; Brovkin etal. 2007; Hansen et al. 2007; Hansen et al. 2011;
Martinez-Boti et al. 2015).
The most important positive feedback mechanismsthat are able to amplify the climate response to primary
forcing are the CO2 exchange between the ocean andthe atmosphere, the water vapor feedback, the methane
release from permafrost peat bogs and from hydrates, allmechanisms related to greenhouse gases, and the surface
albedo feedback. Hansen et al. (2007) revealed a remarkable correspondence between the Vostok Antarctic
temperature and GHG climate forcing, confirming thatthe temperature causes gas changes over several hundred years and suggesting that increasing warming
temperature causes a net release of GHGs from oceans,soils, and the biosphere. In the EPICA ice cores, atmospheric greenhouse gases have strongly covaried with

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Antarctic temperature and δ18O benthic throughout
eight glacial cycles over the last 800 thousand years(Jouzel et al. 2007; Luthi et al. 2008; Loulergue et al.
2008; Berger et al. 2016). Even Schneider et al. (2013) at
EPICA Dome C, during the glacial inception at the endof MIS-5e, confirm a 5-kyr lag in the CO2 decline withrespect to temperatures. They explain this lag by terrestrial carbon release and carbonate compensation.
Brovkin et al. (2007), by simulating the effects ofchanges in atmospheric and oceanic circulations in
response to glacial changes in radiative forcing, find thatthe net effect of five feedback mechanisms related tooceanic circulation and marine biogeochemistry is able
to explain a 65 ppmv drop in glacial atmospheric CO2,
compared to 80–90 ppmv recorded in the ice cores(approximately 70–80% of the total drop). Cooler
oceanic water can absorb more CO2 than warmer water.
As the ocean warms, it releases CO2 to the atmosphere,
providing an amplifying climate feedback that causes
further warming. If the atmosphere is warmed, the saturation vapor pressure increases, and the amount of water
vapor in the atmosphere will tend to increase, leading tofurther warming, because water vapor is an important
greenhouse gas. Liptak and Strong (2016), in a modelingstudy of the Arctic Sea, confirm the existence of a positive thermodynamically driven sea ice-atmosphere feedback up to 35%. A recent study shows that permafrost isa very large carbon pool stored as peat and methane incircumpolar areas (Tarnocai et al. 2009). Lawrence et al.
(2008) suggest that a rapid melting of Arctic Sea ice may
start a feedback loop that rapidly melts Arctic permafrost, triggering further warming feedback by the releaseof methane and carbon dioxide. Methane hydrates are aform of water ice saturated by a large amount ofmethane within its crystal structure. Extremely large
deposits of methane hydrates have been found within
marine sediments along continental margins of the Earth(Dickens 2003). The release of this trapped methane has
the potential to cause a significant increase intemperature as a feedback mechanism triggered by bottom water warming. This hypothesis has been proposed
as a possible explanation of early Cenozoic hyperthermals, such as the Paleocene-Eocene Thermal Maximum(PETM) at approximately 55 Myr, a transient warmingthat shows characteristics of an extreme event ofshort-term positive feedback (Dickens 2003; Zachos etal. 2008). If the planet becomes warmer, the ice sheet
area tends to decrease, exposing a darker surface of landand open water that absorbs more sunlight, causing
additional heating. Instead, during times of global cooling, additional ice increases the reflectivity, which
reduces the absorption of solar radiation and results inmore cooling (albedo effect). According to Hansen et al.
(2013) and Hansen (2013), the net effect of fast positive

Page 29 of 37

feedback is to amplify the global temperature response
by a factor of two to three, and this reinforcement isstate dependent. Even Menviel et al. (2008) state that theatmospheric CO2 response is a delicately balanced sumof the terrestrial and marine inventory changes, emphasizing the importance of the initial climate state on theCO2 response to large-scale ocean circulation changes.
However, exponential growth cannot continue indefinitely because some negative feedback processes act asbalancing mechanisms. The new findings of Early
Pliocene asymptotic decays of Rs (Fig. 11) could represent the evidence of a damping effect of carbon cycle bynegative feedback processes. The continental weathering
is a chemical process of atmospheric CO2 consumptionthat requires tens of thousands of years to act (Walker
and Hays 1981; Schrag et al. 2002; Hansen 2013).
Instead, the iceberg cooling effect is a fast negative
feedback that comes into play as iceberg discharge
reaches a rate that cools the regional ocean surface(Hansen and Sato 2012). The interaction between
vegetation and climate is another interesting component
in climate system dynamics, but the question of theoverall role of the biogeophysical and biogeochemical
feedbacks is not yet solved because both feedbacks canbe positive or negative (Claussen 2009). Nevertheless,shifts in vegetation patterns have presumably provided anegative feedback in the climate system due to theirinfluence on carbon cycle and great areal changes.
During glacial inception, forests were reduced in theirareal extent leading to an increase in atmospheric CO2
concentration, while during deglaciation, atmospheric
CO2 concentration should have decreased due to thebuildup of forests (Claussen 2009). Palynological studies
attest to major changes in vegetation ecosystems both at
glacial and interglacial time scales (Ning et al. 1993;
Elenga et al. 2000; Bertini 2010). During a glacial
climate, tropical rainforests reduce in their extent totropical seasonal forests in tropical lowlands and toxerophytic taxa in tropical highlands (Elenga et al. 2000).
Boreal forests regress equatorward with a compression
and fragmentation of the forest zones, while cold steppe
vegetation expands (Ning et al. 1993; Bertini 2010).
However, after the Zanclean and towards the Late
Pleistocene, a progressive impoverishment of the subtropical forest ecosystem occurs with an increase inherbaceous taxa following the long-term cooling anddrying climate trend (Ning et al. 1993; Bertini 2010).
Another negative feedback mechanism is the increase inprecipitation rate in a greenhouse climate that may
intensify weathering rates as fresh water is flushed more
rapidly through soil and sedimentary rocks (Schrag et al.
2002). Removal of CO2 is accomplished through burialof organic carbon, and increasing river discharge maylead to a more efficient burial of organic carbon

Viaggi Progress in Earth and Planetary Science(2018) 5:81

(Schrag et al. 2002). The oceanic carbon pump by primary productivity can act as possible damping mechanism, but the complex feedback between climate and ocean
require a better understanding. In fact, on orbital timescale, export production by barite accumulation rates was
generally higher during cold periods (Ma et al. 2015). This
interpretation could imply a positive feedback on carbon
cycle that reinforces the cooler orbital response, but it isdifficult to predict the net changes in CO2 efflux (Ma et al.
2015). On the other hand, based on the long-term trends,
export production was on average higher during thewarmer Late Pliocene than the colder Pleistocene (Ma et al.
2015). On the contrary, the depositional conditions of the
Plio-Pleistocene Mediterranean sapropels exhibit warm andwet precession-related changes of phytoplankton productivity by river runoff increasing and bottom water anoxia(Emeis 1998; Gallego-Torres et al. 2007; Moller et al. 2012),
suggesting a negative feedback on carbon cycle (Schraget al. 2002) and a possible damping effect on warming
climate response.
On a geological time scale, almost all these mechanisms are sufficiently “fast” (less than 100 kyr) and canact over an orbital short-term hemicycle like precession
or “semi-precession.” Being geologically fast processes,
the feedback response should have in theory a similarfrequency modulation as the primary astronomical forcing, although the processes should exhibit a slight delay
compared to the latter and to the initial temperaturechange induced by insolation components (Berger andLoutre 1997; Ruggieri et al. 2009; Hansen et al. 2011).
The δ18O and SST lags shown in the present study
support this hypothesis. Moreover, the amplitude of the
feedback response should be linked to the orbital
hemicycle duration and to the scale effect of the surfaces/volumes involved in the processes (Schrag et al.
2002; Menviel et al. 2008; Hansen et al. 2013), a possible
explanation of the nonlinear features shown in Figs. 10
and 11. Therefore, the net cumulative effect of positive
and negative fast feedback mechanisms, added to theinitial temperature change, could lead to the exponential
response of each astronomically driven response as afunction of the LT-MCS at a given time (initial condition), the duration of the orbitally paced feedback processes, and the synergy of scale-dependent processes.
The different growth rates of the nonlinear orbital
responses (Fig. 10) and the asymmetry in Rs (Fig. 11)
suggest a high-efficiency paced system depending onwhether feedbacks amplify (Rs > 0) or suppress (Rs < 0)
response changes in the orbital forcings with astate-dependent magnitude related to the duration timeof the feedback processes. Interestingly, the results ofthe present study indicate an amplified orbital response
during times of a global icehouse state (Mid-Late Pleistocene, with the remarkable exception of obliquity),

Page 30 of 37

which can reach a magnitude of three to five times global δ18O or two to three times tropical SST forcing, anda maximum of two times the damped response for all
signals during the time of a global greenhouse state(Early Pliocene). In this study, it is proposed that orbital
response be considered basically as composite feedback
responses paced by orbital and damped or amplified in arange of − 100 to + 400% the forcing, as function of the
LT-MCS at a given time (initial conditions), the hemicycle duration, and the synergy of the feedback processes. Therefore, the sensitivity of the climate system toorbital forcings could depend on the long-term climate
factors in terms of average atmospheric GHG composition, the average concentration of volcanic aerosols, theglobal mean temperature, the cryosphere average volume, and, presumably, the global forest cover and ocean
primary productivity, factors affecting high-efficiency
feedback mechanisms related to carbon cycle andice-albedo. This interpretation may explain the energy
excess “paradox” of the astronomically paced signals
compared to the small energy of the orbital forcing, especially concerning the eccentricity bands (Imbrie et al.
1993; Wunsch 2004; Berger et al. 2005; Berger et al.
2012). The progressive increase in the strength of the
95-kyr δ18O response from the Zanclean (very weak) to
the Mid-Late Pleistocene may be interpreted as theeffect of increased sensitivity to eccentricity due to thegrowing trend of the global ice volume following thecooling trend. In this time interval, the LT-LR04 component undergoes an average enrichment from 2.99‰(Zanclean) to 4.18‰ (Middle Pleistocene-Holocene) due
to long-term changes in atmosphere composition, andthe climate system may become sensitive to the weakeccentricity insolation by overcoming the critical condition of global mean temperature and average ice volume,a hypothesis already formulated by Raymo (1998) andmany others (Berger et al. 1999; Raymo et al. 2006;
Berger and Yin 2012). The models of Imbrie et al. (2011)
and Abe-Ouchi et al. (2013) also support this hypothesis.
Their studies show that ice volume plays a very important role in pacing the orbital response of the climate
system. Considering orbital signals basically as feedback
responses paced by orbitals may explain why the ~ 100-kyr
eccentricity-driven response is so strong post-MPT, when
the obliquity “ice killing,” which prevents a long-life ice
sheet, was strongly mitigated by damping effect (Fig. 11).
Also, the hysteresis loop of the North American ice sheet
may have contributed to lead a long-life ice sheet.
After the inception of the ice sheet, its mass balanceremains mostly positive through several precession
cycles, whose amplitudes decrease towards an eccentricity
minimum (Abe-Ouchi et al. 2013). The larger the ice
sheet grows and extends towards lower latitudes, thesmaller is the insolation required to make the mass

Viaggi Progress in Earth and Planetary Science(2018) 5:81

balance negative. A long-life ice sheet sensitive to eccentricity during a time of a global icy state may have activated a synergistic response of positive feedback processes
modulated on ~ 100-kyr eccentricity response times. The
negligible response of δ18O in the eccentricity band during
the Zanclean could be due to the higher mean
temperature and the lower mean global ice volume, which
would fail to move permanently to the climate signal due
to obliquity “ice killing,” a powerful forcing during thisperiod. Considering the asymptotic decay of the δ18O andSST Rs up to a maximum of − 100% towards the Early
Pliocene (Fig. 11), a working hypothesis could be a
damping effect due to the action of negative feedbacks
during the time of higher global average temperature,such as, for instance, a higher global forest coverage, agreatest ocean primary productivity, and a more efficientburial of organic carbon. Zanclean was a period of global
greenhouse state, characterized by high atmospheric
GHGs and mean global temperature (Bertini 2010;
Seki et al. 2010; Bartoli et al. 2011). During Zanclean,
northern and central Italy were covered by vegetation
dominated by hygrophilous and thermophilous forest taxa
typical of a humid subtropical to warm-temperate climate. Precipitation was sufficiently high for the persistence of a “broad-leaved evergreen/warm temperate mixed
forest” up to 3.5 Myr (Bertini 2010). Therefore, a wide
forest cover, along with a possible high ocean primary
productivity and organic carbon burial, could increase thequantitative role of the CO2-related negative feedbacks(Emeis 1998; Schrag et al. 2002; Gallego-Torres et al.
2007; Claussen 2009; Moller et al. 2012; Ma et al. 2015)
causing a damping effect of the orbitally paced responses.
In the Mid-Late Pleistocene, an icy state leads to a significant reduction in global forest cover, making these processes quantitatively less important but emphasizing
the vegetation-snow albedo feedback in synergy withthe sea-ice albedo feedback which tends to amplify
Northern Hemisphere and global mean temperaturechanges (Claussen 2009).
Mid-Late Pleistocene obliquity damped responses

To attempt the interpretation of these nonlinear features(Figs. 10 and 11), some possible mechanisms are discussed
in this section as working hypotheses for future studies.
The global δ18O and tropical SST sharply decline tonear-zero in obliquity Rs during the Mid-Late Pleistocene
suggesting an attenuation mechanism of the obliquity
driving force and, consequently, a reduction of the feedback amplification processes linked to the mean icy state.
Understanding this phenomenon is important because thepost-MPT obliquity forcing attenuation could have contributed to the strengthening of the eccentricity response,
favoring a long-life ice growth. Lisiecki and Raymo (2007)
state that the change in δ18O glacial dynamics at ~ 1.4 Myr

Page 31 of 37

is associated with an abrupt decline in 41-kyr power and adecrease in modulation sensitivity to obliquity. Interestingly, the shape of the δ18O and SST components-3-4 after
the MPT (Figs. 4c and 8c), characterized by the envelope of two to three low-amplitude obliquity cycles ina weak ~ 90–93-kyr framework, suggests a connection
between eccentricity and obliquity, similar to that proposed by Huybers (2007) of obliquity-cycle skipping. The
obliquity-pacing results indicate that the 100-kyr variability could result from the skipping of one or two obliquity
beats, corresponding to 80-kyr or 120-kyr glacial cycles
that, on average, give 100-kyr periodicity, resulting fromlong-term changes in the climate system, such as ice sheet
growth (Huybers 2007). Why does obliquity-cycle skipping occur? Climate friction is a dissipative feedback between obliquity variations and climate which may cause asecular drift of the spin axis (Laskar et al. 2004). In response to obliquity forcing, glacial and interglacial conditions drive the redistribution of the ice/water mass andthe isostatic adjustment to the surface loading, affecting
the dynamical ellipticity of the Earth (obliquity–oblatenessfeedback) (Rubincam 1993, 1995; Levrard and Laskar
2003; Laskar et al. 2004). Delayed responses in bothclimatic and viscous relaxation processes may introduce asecular term in the obliquity evolution (Rubincam 1995).
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 theentire system (Laskar et al. 2004). Therefore, simplifying
assumptions have to be made to estimate the magnitude
and the direction of the secular drift (increase or decreasein the obliquity). Specifically, the sign of the secular obliquity change is particularly sensitive to the obliquity-lag
estimates and changes of oblateness (Williams et al. 1998;
Levrard and Laskar 2003). Levrard and Laskar (2003)
showed that climate friction impact is positive and likely
negligible over the last 3 Myr, based on the model assumptions that the ice sheet response time-lag could not
be longer than the obliquity period and explicitly assumed
it to be 8 kyr. However, the secular change of obliquity becomes negative for low obliquity lag response (< 4.2 kyr)and high positive change of oblateness (Levrard and Laskar
2003). Response lags are difficult to estimate. According to
Lisiecki and Raymo (2005a), the δ18O lag in obliquity
reaches values of ~ 6.8 kyr in the Pleistocene. The results ofthe present study indicate a δ18O obliquity lag of 7.4 kyrpost-MPT (Table 8), close to the value used by Levrard andLaskar (2003). In contrast, the obliquity lag estimation of
Hilgen et al. (1993) provides a lower value of 5.6 kyr withrespect to the Late Pleistocene if the phase relation withprecession is kept constant. Considering the margin oferror of these estimates, it is currently not possible to establish whether the obliquity lags may be lower post-MPT tojustify the observed obliquity damping effect. Furthermore,Viaggi Progress in Earth and Planetary Science(2018) 5:81

other simplified assumptions in the dissipative equations
model, together with age model bias, could affect the final
results. A recent study by Adhikari and Ivins (2016) using
space geodetic and satellite gravimetric present-day data
has shown that global-scale continent-ocean mass transportchanges (terrestrial water storage and global cryosphere)
explain nearly the entire amplitude (83 ± 23%) and mean
directional shift (within 5.9° ± 7.6°) of the observed polar
motion on a very short present-day time scale. The study ofAdhikari and Ivins (2016) suggests a high sensitivity of thepoles to the mass transport changes and to the isostatic
rebound. In conclusion, new studies are needed to clarify
the effective role of climate friction or to identify othermechanisms that might be responsible for this unusualtransition in obliquity response.
Tropical SST coupled with benthic δ18O cycles

Despite the differences in signal-to-noise ratio andoceanographic context between the LR04 δ18O and the
ODP Site 846 SST records (global benthic stack andupwelling tropical single record), the SSA results show
the tropical SST tightly coupled with benthic δ18O overthe past five million years, both at long and orbital timescales, in agreement with the original work of Herbert etal. (2010a). Some differences certainly exist (e.g., a strong
cooling event at ~ 2.1 kyr or MBO warming less marked
in the SST signal; the tropical trend cooling slightly
negative shifted compared to global ΔSST, nonethelesssubparallel in time) (Fig. 14b), but they may be explainedin terms of temporal and regional upwelling effects at
the ODP Site 846. These similarities suggest common
causes for both the δ18O long-term trend and cycles andtheir expression in the tropical SST. The most likely
candidate for some orbitally paced cycles lies inhigh-latitude oceanographic processes tied to the global
carbon cycle (Liu and Herbert 2004; Herbert et al.
2010a). These processes drive the oceanic circulation
through changes in deep water properties, affecting vertical density stability and upwelling in the tropicaloceans. Changes in the atmospheric GHG composition
controlled by insolation processes in the high latitudes
may also act to imprint obliquity on the tropical Pacificwhere the local direct effect of obliquity is weak (Liuand Herbert 2004; Medina-Elizalde and Lea 2005;
Masson-Delmotte et al. 2006). The δ18O-SST delta-lag
in obliquity became large after the beginning of the MPT,
suggesting a signal of global ice volume growth (Table 8).
Such observations, together with the high coherence
between SST and δ18O obliquity responses (Table 9) may
suggest remote connections from high latitudes towards
the tropics. These interpretations could explain why thetropical SST obliquity component contains a small variancein comparison to that of δ18O (5.5% vs. 9.9%). The smaller
coherency between tropical SST and δ18O precession

Page 32 of 37

components (Table 9) might indicate that the response in
the EEP was caused by low-latitude precession forcing andregional effects and was less coupled to the high-latitude
ice volume (Liu and Herbert 2004). Tachikawa et al. (2014)
found residual SST variability characterized by a 23-kyr
periodicity in the western Pacific warm pool (WPWP) region. They concluded a wind-driven northward surface
current advection towards the WPWP to be responsible forthe regional SST precessional periodicity, attributed to meridional shifts of the Intertropical Convergence Zone. Even
the “half-precession” cycle of equatorial origin (Berger et al.
2006), isolated in this study by SSA from both tropical SST
and global δ18O records, may be transmitted to highlatitudes via advective transport (Hagelberg et al. 1994).
The last interglacial half-precession cycle from theNorthwestern Chinese Loess Plateau, interpreted as a direct
response to low-latitude forcing through its modulation on
the East Asian summer monsoon (Sun and Huang 2006),
supports this interpretation. The coherency increase of theshort eccentricity’s δ18O and SST signals and the maximum
delta-lag in eccentricity between δ18O and SST (4.31 vs.
0.44 kyr) towards the Late Pleistocene (Table 8) may reflect
the progressive rise of eccentricity power and an increased
contribution of the high-latitude δ18O ice volume component over the water temperature signal. Berger et al. (2006)
show the presence of significant 100-kyr and 11-kyr cycles
in the amplitude of the seasonal cycle in the intertropicalregions, but the amplitude decreases rapidly when moving
away from the equator. The small differences in variance
between the short eccentricity’s SST and δ18O components(5.1% vs. 6.5%) might depend on these considerations.
Conversely, the big difference in variance of the closelyenveloped long-term cycles (Earth’s orbit secular frequencies, obliquity modulation cycles, and long eccentricity
cycles) between SST (10.3%) and δ18O (4.2%) components
could suggest a greater sensitivity of the tropics to theselong-term cycles. The tropical oceans should be shieldedfrom feedback processes that produce large temperaturesensitivity in the high latitudes, such as ice-albedo
feedbacks, sea ice-atmosphere feedback, and permafrost
feedback (Lawrence et al. 2008; Tarnocai et al. 2009;
Herbert et al. 2010a; Liptak and Strong 2016). With the
exception of obliquity, the marked difference in eccentricity
and precession Rs in the Late Pleistocene between
δ18O (200–400%) and SST (100–170%) (Fig. 11) may
be due to a lower synergy of the positive feedback
mechanisms in tropical regions. In summary, both
remote connections from high-latitude and direct
low-latitude orbital forcings could couple Plio-Pleistocene
global δ18O and tropical SST orbital cycles throughoceanic and atmospheric circulations. Finally, the remarkable constancy (~ 76%) of the Plio-Pleistocene trend
variance among global δ18O, local tropical SST (Site 846),
and global SST (ΔSST stack) might reflect theViaggi Progress in Earth and Planetary Science(2018) 5:81

latitude-independent nature of the long-term factors controlling the atmosphere composition and global mean
temperature.

Conclusions
The main conclusions of this study can be summarized
as follows:
1) A detailed evaluation is provided of the variancepaced by long-period orbital modulation, short eccentricity, obliquity, precession, and half-precession
cycles, estimated for the δ18O at approximately
4.2%, 6.5%, 9.9%, 2.0%, and 0.6%, respectively, andfor the tropical SST at about 10.3%, 5.1%, 5.5%,
1.6%, and 0.5%, respectively, both resulting in total
orbital contributions of ~ 23.0%. On long timescales, the Plio-Pleistocene astronomical forcing
could be considered a random walk at zero-sum inwhich the orbital responses oscillate around a basically time-invariant quasi-normal distributed mean,
supporting the hypothesis that orbital forcings did
not have a major role in guiding the PlioPleistocene climate system response.
2) The exponential growth of the orbital variancetowards the Late Pleistocene is function of the LTMCS at a given time and regards all δ18O and SST
components with an increasing rate related to theorbital period, from half-precession to eccentricity.
The asymmetry of the Rs suggests an amplified orbital response during times of a global icehousestate (Mid-Late Pleistocene, with the exception of
obliquity), which can reach a magnitude of + 200 to+ 400% global δ18O forcing, or + 100 to + 170% thetropical SST forcing, and a minimum up to − 100%
the damped response for all signals during time of a
global greenhouse state (Early Pliocene).
3) It is proposed that orbital responses be considered
as composite feedback lagged responses paced byorbitals and damped or amplified in a range of −100 to + 400% the forcing, depending on the initialconditions set by LT-MCS, the hemicycle duration,and the synergy of the feedback processes related tocarbon cycles and albedo effects. This interpretation
may explain the energy excess “paradox” of theastronomically paced signals compared to the smallenergy of the orbital forcing, especially concerning
the eccentricity bands.
4) The asymptotic decay of the δ18O and SST Rs up to− 100% towards the Early Pliocene could be due toa damping effect of negative feedbacks during atime of higher mean global temperature. Zanclean
was a period of global greenhouse state, and a wide
forest cover along with a possible high ocean
primary productivity and a more efficient burial ofPage 33 of 37

organic carbon could have increased thequantitative role of the CO2-related negative
feedbacks, causing a damping effect of the orbitally
paced responses.
5) The global δ18O and tropical SST anomalous
sharply decline to near-zero in obliquity Rs during
the Late Pleistocene suggesting an attenuation
mechanism of the obliquity driving force and a reduction of the related feedback amplification processes. The post-MPT obliquity damping may havecontributed to the strengthening of the eccentricity
response by mitigating the obliquity “ice killing,” favoring a 100-kyr long-life hysteresis ice growth volume. A long-life ice sheet sensitive to eccentricity
during the time of a global icy state may have activated a synergistic response of positive feedback
processes modulated on short eccentricity response
times. In the Mid-Late Pleistocene, an icy stateleads to a significant reduction in global forest cover
making its negative feedbacks quantitatively less important but emphasizing the positive feedback processes, such as the vegetation-snow albedo insynergy with the sea-ice albedo and permafrost
feedbacks. However, the origin of the post-MPT obliquity attenuation mechanism is an open question.
6) The tropical SST orbital cycles are tightly coupledwith the global benthic δ18O cycles over the past
five million years, suggesting both Plio-Pleistocene
oceanic and atmospheric remote connections fromhigh-latitude and direct low-latitude orbital forcings.
7) The tropical SST, the global SST, and the global
δ18O trend components, each explaining ~ 76% of
the Plio-Pleistocene climate variance, could be
related to the long-term pCO2 proxies, suggesting
non-orbital latitude-independent control factors.
These trend components significantly change theclimate mean using arbitrary time intervals, suggesting the Plio-Pleistocene mean climate state waschanged primarily by long-term carbon cycle
controlling the Earth’s atmospheric composition.
The geotectonic forcing appears to be a plausible
hypothesis to explain the long-term changes inatmospheric composition (GHGs and volcanic
aerosols) through direct (CO2 outgassing and explosive volcanism) and indirect (orography and erosion, paleogeography, oceanic paleocirculation, andocean fertilization) control by plate tectonics.
8) The NHG gradually developed in four phases(subtrends) of progressive lowering of the averagetemperatures and ice volume growth, triggered bylong-term changes in the atmospheric composition.
Nevertheless, transient interaction effects with very
long-term orbital cycles, expression of the chaotic
behavior of the solar system, may have contributed

Viaggi Progress in Earth and Planetary Science(2018) 5:81

to setting these transitions. After MIS-16, thecooling trend was broken by a wide swing of δ18O
depletion and SST recovery, a relative stasis of ice
growth and temperature cooling during the last
600 kyr (MBO), which could be related to midterm GHG recovery.
The idea of the geotectonic forcing is the most attractive
hypothesis by a long time, and the recent improvements inthe seafloor spreading rate evaluation lead in the right
direction. However, further studies are needed to better
estimate high-resolution seafloor spreading rate and toassess the role of the volcanic CO2 degassing and the explosive volcanic aerosols in the budget of the carbon cycle.
Even much higher resolution CO2 proxies, covering a longtime period, are necessary to better address carbon cycle.
Nevertheless, the relevance of the primary long-term climate system, controlling the Earth’s atmospheric composition and modifying the average global climate, must be
underlined. Finally, singular spectrum analysis provides avaluable tool in cyclostratigraphy with the remarkableadvantage of separating full-resolution time series by
variance strength.

Additional files
Additional file 1: Methodology. (DOC 1687 kb)
Additional file 2: Dataset. (XLS 3816 kb)
Additional file 3: High-resolution figures. (PDF 6502 kb)

Abbreviations
ANOVA: Analysis of variance; EECO: Early Eocene Climatic Optimum;
EEP: Eastern Equatorial Pacific; EPICA: European Project for Ice Coring in
Antarctica; FFS: Fourier frequency spectrum; FFT: Fast Fourier transform;
GASR: Global average spreading rate; GCM: General circulation model;
GHG: Greenhouse gas; GSSP: Global Stratotype Section and Point;
HT846: Herbert et al. (2010b) sea surface temperature of ODP Site 846;
INHG: Intensification of Northern Hemisphere Glaciation; LGM: Last Glacial
Maximum; LIP: Large igneous province; LT-HT846: Long-term HT846 SST
component; LT-LR04: Long-term LR04 δ18O component; LT-MCS: Long-termmean climate state; MBE: Mid-Brunhes Event; MBO: Mid-Brunhes Oscillation;
MIS: Marine Isotope Stage; MMCO: Middle Miocene Climatic Optimum;
MPT: Mid-Pleistocene Transition; NASSR: North Atlantic seafloor spreading
rate; NHG: Northern Hemisphere Glaciation; ODP: Ocean Drilling Program;
ONHG: Onset of the Northern Hemisphere Glaciation; PETM: PaleoceneEocene Thermal Maximum; Rs: Response sensitivity; SSA: Singular spectrum
analysis; SST: Sea surface temperature; TISA: Time-integral squared amplitude;
WPWP: Western Pacific warm pool
Acknowledgements
Special thanks to Gian Battista Vai for the useful suggestions and for theencouragement to publish the study. I am grateful to Francesco Mulargia forthe methodological review. Thanks also to Nicola Scafetta for the useful
discussion. I thank Lorraine Lisiecki for kindly providing me the δ13C-based
pCO2 proxy data. Special thanks to Alessandra Negri for improving theEnglish language and other constructive comments.
Funding
This research did not receive any specific grant from funding agencies in the
public, commercial, or not-for-profit sectors.

Page 34 of 37

Availability of data and materials
Details on the methodology used in this study and a discussion of thequality and reliability of the time series can be found in Additional file 1. The
complete SSA-components dataset and a complete set ofhigh-resolution figures are also provided:– Additional file 1 - Methodology.doc– Additional file 2 - Dataset.xls– Additional file 3 - High resolution figures.pdf
Authors’ contributions
The author read and approved the final manuscript.
Competing interests
The author declares that he has no competing interests.

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.
Received: 16 April 2018 Accepted: 29 October 2018

References
Abe-Ouchi A, Saito F, Kawamura K, Raymo ME, Okuno J, Takahashi K, Blatter H(2013) Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet
volume. Nature 500. https://doi.org/10.1038/nature12374
Adhikari S, Ivins ER (2016) Climate-driven polar motion: 2003–2015. Sci Adv 2.
https://doi.org/10.1126/sciadv.1501693
Aitchison JC, Ali JR, Davis AM (2007) When and where did India and Asia collide?
J Geophys Res 112:B05423. https://doi.org/10.1029/2006JB004706
Archer D, Martin P, Buffett B, Brovkin V, Rahmstorf S, Ganopolski A (2004) The
importance of ocean temperature to global biogeochemistry. Earth Planet
Sci Lett 222(2004):333–348. https://doi.org/10.1016/j.epsl.2004.03.011
Bartoli G, Hönisch B, Zeebe RE (2011) Atmospheric CO2 decline during the
Pliocene intensification of northern hemisphere glaciations.
Paleoceanography 26:PA4213. https://doi.org/10.1029/2010PA002055
Bartoli G, Sarnthein M, Weinelt M, Erlenkeuser H, Garbe-Schonberg D, Lea DW(2005) Final closure of Panama and the onset of northern hemisphere
glaciation. Earth Planet Sci Lett 237(2005):33–44. https://doi.org/10.1016/j.
epsl.2005.06.020
Becker TW, Conrad CP, Buffett B, Dietmar MR (2009) Past and present
seafloor age distributions and the temporal evolution of plate tectonic
heat transport. Earth Planet Sci Lett 278(2009):233–242. https://doi.org/10.
1016/j.epsl.2008.12.007
Beerling DJ, Royer DL (2011) Convergent Cenozoic CO2 history. Nat Geosci 4, pp.
418-420.
Berger A, Li XS, Loutre MF (1999) Modelling northern hemisphere ice volume
over the last 3 Myr. Quat Sci Rev 18(1999):1–11
Berger A, Loutre MF (1997) Long-term variations in insolation and their effects on
climate, the LLN experiments. Surv Geophys 18:147–161
Berger A, Loutre MF, Mélice JL (2006) Equatorial insolation: from precession
harmonics to eccentricity frequencies. Clim Past Discuss 2:519–533 https://doi.org/10.5194/cp-2-131-2006
Berger A, Mélice JL, Loutre MF (2005) On the origin of the 100-kyr cycles in theastronomical forcing, Paleoceanography, 20:PA4019, https://doi.org/10.1029/
2005PA001173.
Berger A, Yin Q (2012) Modeling the interglacials of the last 1 million years. In:
Berger et al (eds) Climate change: inferences from paleoclimate and regional
aspects. Springer-Verlag, Wien 2012.
Berger A, Yin QZ, Herold N (2012) MIS-11 and MIS-19, analogs of our Holocene
interglacial. In: 3rd International Conference on Earth System Modelling,
Hamburg, https://meetingorganizer.copernicus.org/3ICESM/3ICESM-11.pdf
Berger A et al (2016) Past interglacials working group of PAGES, interglacials of
the last 800,000 years. Rev Geophys 54:162–219. https://doi.org/10.1002/
2015RG000482
Berger AL (1978) Long-term variations of daily insolation and quaternary climaticchange. J Atmos Sci 35:2362–2367
Berner RA, Kothavala Z (2001) GEOCARB III: a revised model of atmospheric CO2
over phanerozoic time. Am J Sci 301(2):182–204

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Berner RA, Lasaga AC, Garrels RM (1983) The carbonate-silicate geochemical cycle
and its effect on atmospheric carbon dioxide over the past 100 million years.
Am J Sci 283:641–683
Bertini A (2010) Pliocene to Pleistocene palynoflora and vegetation in Italy - stateof the art. Quat Int 225(2010):5–24
Blanc PL, Fontugne MR, Duplessy JC (1983) The time-transgressive initiation of
boreal ice-cap: continental and oceanic evidence reconciled. Palaeogeogr
Palaeoclimatol Palaeoecol 42:211–224
Boulila S, Galbrun B, Miller KG, Pekar SF, Browning JV, Laskar J, Wright JD (2011)
On the origin of Cenozoic and Mesozoic “third-order” eustatic sequences.
Earth Sci Rev 109:94–112
Brady PV, Gislason SR (1997) Seafloor weathering controls on atmospheric CO2
and global climate. Geochimica Et Cosmochimica Acta 61(5):965–973
Brovkin V, Ganopolski A, Archer D, Rahmstorf S (2007) Lowering of glacial
atmospheric CO2 in response to changes in oceanic circulation and marine
biogeochemistry. Paleoceanography 22:PA4202. https://doi.org/10.1029/
2006PA001380
Clark PU, Archer D, Pollard D, Blum JD, Rial JA, Brovkin V, Mix AC, Pisias NG, Roy
M (2006) The middle Pleistocene transition: characteristics, mechanisms, andimplications for long-term changes in atmospheric pCO2. Quat Sci Rev.
https://doi.org/10.1016/j.quascirev.2006.07.008
Claussen M (2009) Late Quaternary vegetation-climate feedbacks. Clim Past 5:
203–216, https://doi.org/10.5194/cp-5-203-2009
Cogné JP, Humler E (2006) Trends and rhythms in global seafloor generation
rate. Geochem Geophys Geosyst 7:Q03011. https://doi.org/10.1029/
2005GC001148
Coltice N, Seton M, Rolf T, Müller RD, Tackley PJ (2013) Convergence of tectonic
reconstructions and mantle convection models for significant fluctuations inseafloor spreading. Earth Planet Sci Lett 383(2013):92–100. https://doi.org/10.
1016/j.epsl.2013.09.032
Conrad CP, Lithgow-Bertelloni C (2007) Faster seafloor spreading and lithosphere
production during the mid-Cenozoic. Geology 35(1):29–32. https://doi.org/10.
1130/G22759A
Dickens GR (2003) Rethinking the global carbon cycle with a large,
dynamic and microbially mediated gas hydrate capacitor. Earth Planet
Sci Lett 213(2003):169–183. https://doi.org/10.1016/S0012821X(03)00325-X
Elenga H et al (2000) Pollen-based biome reconstruction for southern Europe
and Africa 18,000 years ago. J Biogeogr 27:621–634
Elsner JB, Tsonis AA (1996) Singular spectrum analysis: a new tool in time series
analysis. Springer
Emeis K.-C., Sakamoto T. (1998) - The Sapropel Theme of Leg 160. In: Robertson
A.H.F., Emeis K.-C., Richter C. & Camerlenghi A. (Eds.). Proceedings of theocean drilling program, scientific results. 160
Gallego-Torres D, Martínez-Ruiz F, Paytan A, Jiménez-Espejo FJ, OrtegaHuertas M (2007) Pliocene–Holocene evolution of depositional
conditions in the eastern Mediterranean: role of anoxia vs. productivity
at time of sapropel deposition. Palaeogeogr Palaeoclimatol Palaeoecol
246(2007):424–439. https://doi.org/10.1016/j.palaeo.2006.10.008
Ghil M, Allen RM, Dettinger MD, Ide K, Kondrashov D, Mann ME, Robertson A,
Saunders A, Tian Y, Varadi F, Yiou P (2002) Advanced spectral methods forclimatic time series. Rev Geophys 40(1):3.1–3.41. https://doi.org/10.1029/
2001RG000092
Gibbard PL, Head MJ (2009) IUGS ratification of the Quaternary system/period
and the Pleistocene series/epoch with a base at 2.58 Myr. Quaternaire20(4):411–412
Hagelberg TK, Bond G, deMenocal P (1994) Milankovitch band forcing of subMilankovitch climate variability during the Pleistocene. Paleoceanography9(4):545–558. https://doi.org/10.1029/94PA00443
Hansen J. (2013) - An old story, but useful lessons. http://www.columbia.edu/~
jeh1/mailings/2013/20130926_PTRSpaperDiscussion.pdf
Hansen J, Sato M (2012) Paleoclimate implications for human-made climatechange. In: Berger A, Mesinger F, Sijacki D (eds) Climate change. Inferencesfrom paleoclimate and regional aspects. Springer-Verlag, Wien 2012, p 350,
https://doi.org/10.1007/978-3-7091-0973-1_2
Hansen J, Sato M, Kharecha P, Russell G, Lea DW, Siddall M (2007) Climatechange and trace gases. Phil Trans R Soc A 365:1925–1954. https://doi.org/10.
1098/rsta.2007.2052
Hansen J, Sato M, Kharecha P, von Schuckmann K (2011) Earth’s energy
imbalance and implications. Atmos Chem Phys 11:13421–13449, https://doi.
org/10.5194/acp-11-13421-2011

Page 35 of 37

Hansen J, Sato M, Russell G, Kharecha P (2013) Climate sensitivity, sea level andatmospheric carbon dioxide. Phil Trans R Soc A 371:20120294. https://doi.
org/10.1098/rsta.2012.0294
Hansen et al (2005) Efficacy of climate forcings. J Geophys Res 110:D18104.
https://doi.org/10.1029/2005JD005776
Hassani H (2007) Singular spectrum analysis: methodology and comparison. J
Data Sci 5(2007):239–257
Head M.J. & Gibbard P.L. (2005) - Early–Middle Pleistocene transitions. An
overview and recommendation for the defining boundary. From:
Head MJ & Gibbard PL (eds). Early–middle Pleistocene transitions:the land–ocean evidence. Geological Society, London, Special
Publications, 247, 1–18. 0305–8719/05/$15© The Geological Societyof London 2005
Head MJ, Pillans B, Farquhar SA (2008) The Early–Middle Pleistocene transition:
characterization and proposed guide for the defining boundary. Episodes
31(2):255–259
Herbert TD, Peterson LC, Lawrence KT, Liu Z (2010a) Tropical ocean temperatures
over the past 3.5 million years. Science 328:1530–1534. https://doi.org/10.
1126/science.1185435
Herbert T.D., Peterson L.C., Lawrence K.T. & Liu Z. (2010b) - Plio-Pleistocene
tropical alkenone SST reconstructions LAST UPDATE: 6/2010, NOAA/NCDC
Paleoclimatology Program. ftp://ftp.ncdc.noaa.gov/pub/data/paleo/
contributions_by_author/herbert2010/herbert2010.txt
Hilgen FJ, Lourens LJ, Berger A, Loutre MF (1993) Evaluation of the astronomicallycalibrated time scale for the Late Piocene and the earliest Pleistocene.
Paleoceanography 8:549–565
Hönisch B, Hemming NG, Archer D, Siddall M, McManus JF (2009) Atmosphericcarbon dioxide concentration across the mid-Pleistocene transition. Science
324:1551. https://doi.org/10.1126/science.1171477
Huybers P (2007) Glacial variability over the last two million years: an extended
depth-derived agemodel, continuous obliquity pacing, and the Pleistocene
progression. Quat Sci Rev 26(2007):37–55
Imbrie J, Berger A, Boyle EA, Clemens SC, Duffy A, Howard WR, Kukla G,
Kutzbach J, Martinson DG, Mcintyre A, Mix AC, Molfino B, Morley JJ,
Peterson LC, Pisias NG, Prell WL, Raymo ME, Shackleton NJ, Toggweile JR(1993) On the structure and origin of major glaciation cycles 2. The
100,000-year cycle. Paleoceanography 8(6):699–735
Imbrie JZ, Imbrie-Moore A, Lisiecki LE (2011) A phase-space model for
Pleistocene ice volume. Earth Planet Sci Lett 307:94–102
Jicha BR, Scholl DW, Rea DK (2009) Circum-Pacific arc flare-ups and global
cooling near the Eocene Oligocene boundary. Geology 37(4):303–306.
https://doi.org/10.1130/G25392A.1
Jouzel J et al (2007) Orbital and millennial Antarctic climate variability over the past
800,000 years. Science 317:793–796. https://doi.org/10.1126/science.1141038
Joyce JE, Tjalsma LRC, Prutzman JM (1990) High-resolution planktic stable isotope
record and spectral analysis for the last 5.35 M.Y.: Ocean Drilling Program site
625 northeast Gulf of Mexico. Paleoceanography 5:507–529
Kent DV, Muttoni G (2013) Modulation of Late Cretaceous and Cenozoic climate
by variable drawdown of atmospheric pCO2 from weathering of basaltic
provinces on continents drifting through the equatorial humid belt. Clim
Past 9:525–546. https://doi.org/10.5194/cp-9-525-2013
Kominz M, Pisias N (1979) Pleistocene climate: deterministic or stochastic?
Science 204(lssue 4389):171–173. https://doi.org/10.1126/science.204.4389.171
Laskar J, Fienga A, Gastineau M, Manche H (2011) La2010: a new orbital solution
for the long-term motion of the Earth. Astron Astrophys 532:A89. https://doi.
org/10.1051/0004-6361/201116836
Laskar J, Joutel F, Boudin F (1993) Orbital, precessional and insolation quantities
for the Earth from −20 Myr to +10 Myr. Astron Astrophys 270:522–533
Laskar J, Robutel P, Joutel F, Gastineau M, Correia ACM, Levrard B (2004) A longterm numerical solution for the insolation quantities of the Earth. AstronAstrophys 428:La2004
Lawrence DM, Slater AG, Tomas RA, Holland MM, Deser C (2008) Accelerated
Arctic land warming and permafrost degradation during rapid sea ice
loss. Geophys Res Lett 35:L11506. https://doi.org/10.1029/2008GL033985
Le Breton E, Cobbold PR, Dauteuil O, Lewis G (2012) Variations inamount and direction of seafloor spreading along the Northeast
Atlantic Ocean and resulting deformation of the continental margin ofNorthwest Europe. Tectonics 31:TC5006. https://doi.org/10.1029/
2011TC003087
Levrard B, Laskar J (2003) Climate friction and the Earth’s obliquity. Geophys J Int
154:970–990

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Li G, Elderfield H (2013) Evolution of carbon cycle over the past 100 million years.
Geochim Cosmochim Acta 103(2013):11–25. https://doi.org/10.1016/j.gca.
2012.10.014
Liptak J, Strong C (2016) A modeling investigation of the Arctic Sea iceatmosphere feedback. Clim Dyn 47(Issue 7):2471–2480
Lisiecki LE (2010a) Links between eccentricity forcing and the 100,000-year glacial
cycle. Nat Geosci 3. https://doi.org/10.1038/NGEO828
Lisiecki LE (2010b) A benthic δ13C-based proxy for atmospheric pCO2 over thelast 1.5 Myr. Geophys Res Lett 37:L21708. https://doi.org/10.1029/
2010GL045109
Lisiecki LE, Raymo ME (2005a) A Pliocene-Pleistocene stack of 57 globally
distributed benthic δ18O records. Paleoceanography 20:PA1003. https://doi.
org/10.1029/2004PA001071
Lisiecki L.E. & Raymo M.E. (2005b) Appendix 1 - Global Plio-Pleistocene stack ofbenthic oxygen isotope records. DOI:https://doi.org/10.1594/PANGAEA.
701576. DOI:https://doi.org/10.1029/2004PA001071. In supplement to: Lisiecki& Raymo (2005)
Lisiecki LE, Raymo ME (2007) Plio–Pleistocene climate evolution; trends andtransitions in glacial cycle dynamics. Quat Sci Rev 26(2007):56–69. https://doi.
org/10.1016/j.quascirev.2006.09.005
Liu Z, Herbert TD (2004) High-latitude influence on the eastern equatorial Pacificclimate in the early Pleistocene epoch. Nature 427:720–723
Lo L, Chang S, Wei K, Lee S, Ou T, Chen Y, Chuang C, Mii H, Burr GS, Chen M,
Tung Y, Tsai M, Hodell DA, Shen C (2017) Nonlinear climatic sensitivity togreenhouse gases over past 4 glacial/interglacial cycles. Nat Sci Rep 7:4626.
https://doi.org/10.1038/s41598-017-04031-x
Loulergue L, Schilt A, Spahni R, Masson-Delmotte V, Blunier T, Lemieux B, Barnola
JM, Raynaud D, Stocker TF, Chappellaz J (2008) Orbital and millennial-scalefeatures of atmospheric CH4 over the past 800,000 years. Nature 453. https://doi.org/10.1038/nature06950
Lourens LJ, Hilgen FJ (1994) Long-period orbital variations and their relation
to third-order eustatic cycles and the onset of major glaciations 3.0
million years ago. In: Lourens LJ (ed) Astronomical forcing ofMediterranean climate during the last 5.3 million years. Thesis Utrecht
University, ISBN 90-393-0754-7
Lunt DJ, Foster GL, Haywood AM, Stone EJ (2008) Late Pliocene Greenlandglaciation controlled by a decline in atmospheric CO2 levels. Nature 454:
1102–1106. https://doi.org/10.1038/nature07223
Lunt DJ, Haywood AM, Schmidt GA, Salzmann U, Valdes PJ, Dowsett HJ, Loptson CA(2012) On the causes of mid-Pliocene warmth and polar amplification. Earth
Planet Sci Lett 321-322(2012):128–138. https://doi.org/10.1016/j.epsl.2011.12.042
Luthi D, Le Floch M, Bereiter B, Blunier T, Barnola JM, Siegenthaler U, Raynaud D,
Jouzel J, Fischer H, Kawamura K, Stocker TF (2008) High-resolution carbon
dioxide concentration record 650,000–800,000 years before present. Nature
453. https://doi.org/10.1038/nature06949
Ma Z, Ravelo AC, Liu Z, Zhou L, Paytan A (2015) Export production fluctuations in
the eastern equatorial Pacific during the Pliocene-Pleistocene: reconstruction
using barite accumulation rates. Paleoceanography 30. https://doi.org/10.
1002/2015PA002860
Martinez-Boti MA, Foster GL, Chalk TB, Rohling EJ, Sexton PF, Lunt DJ, Pancost
RD, Badger MPS, Schmidt DN (2015) Plio-Pleistocene climate sensitivity
evaluated using high-resolution CO2 records. Nature 518. https://doi.org/10.
1038/nature14145
Masson-Delmotte et al (2006) Past temperature reconstructions from deep ice
cores. Relevance for future climate change. Clim Past 2:145–165 www.climpast.net/2/145/2006/
Medina-Elizalde M, Lea DW (2005) The mid-Pleistocene transition in the tropical
Pacific. Science 310:1009–1012
Menviel L, Timmermann A, Mouchet A, Timm O (2008) Meridional
reorganizations of marine and terrestrial productivity during Heinrich events.
Paleoceanography 23:PA1203. https://doi.org/10.1029/2007PA001445
Merkouriev S, DeMets C (2014) High-resolution Quaternary and Neogene
reconstructions of Eurasia-North America plate motion. Geophys J Int.
https://doi.org/10.1093/gji/ggu142
Moller T, Schulz H, Hamann Y, Dellwig O, Kucera M (2012) Sedimentology andgeochemistry of an exceptionally preserved last interglacial sapropel S5 in
the Levantine Basin (Mediterranean Sea). Mar Geol 291-294(2012):34–48.
https://doi.org/10.1016/j.margeo.2011.10.011
Mosar J, Lewis G, Torsvik TH (2002) North Atlantic Sea-floor spreading rates:
implications for the Tertiary development of inversion structures of theNorwegian–Greenland Sea. J Geol Soc London 159(2002):503–515

Page 36 of 37

Mudelsee M, Raymo ME (2005) Slow dynamics of the Northern Hemisphere
glaciation. Paleoceanography 20. https://doi.org/10.1029/2005PA001153
Ning S, Jia-Xin C, Königsson LK (1993) Late Cenozoic vegetational history and the
Pliocene-Pleistocene boundary in the Yushe Basin, SE Shanxi, China. Grana32(4-5):260–271. https://doi.org/10.1080/00173139309429990
Pagani M, Liu Z, LaRiviere J, Ravelo AC (2010) High Earth-system climate
sensitivity determined from Pliocene carbon dioxide concentrations. Nat
Geosci 3. https://doi.org/10.1038/NGEO724
Raymo ME (1994) The initiation of northern hemisphere glaciation. Annu Rev
Earth Planet Sci 22:353–383
Raymo ME (1998) Glacial puzzles. Science 281(5382):1467–1468
Raymo ME, Lisiecki LE, Nisancioglu KH (2006) Plio-Pleistocene ice volume,
Antarctic climate, and the global δ18O record. Science 313:492–495
Raymo ME, Ruddiman WF (1992) Tectonic forcing of late Cenozoic climate.
Nature 359:117–122
Rial JA, Pielke RA Sr, Beniston M, Claussen M, Canadell J, Cox P, Held H, deNoblet-Ducoudré N, Prinn R, Reynolds JF, Salas JD (2004) Nonlinearities,
feedbacks and critical thresholds within the Earth’s climate system. Clim
Change 65(1–2):11–38
Rowley DB (2002) Rate of plate creation and destruction 180 Myr to present.
Geol Soc Am Bull 114:927–933
Royer DL (2006) CO2-forced climate thresholds during the Phanerozoic. Geochim
Cosmochim Acta 70(2006):5665–5675
Rubincam DP (1993) The obliquity of Mars and climate friction. J Geophys Res 98:
10827–10832
Rubincam DP (1995) Has climate changed the Earth’s tilt ? Paleoceanography
10(3):365–372
Ruggieri E, Herbert T, Lawrence KT, Lawrence CE (2009) Change point methodfor detecting regime shifts in paleoclimatic time series. Application to δ18O
time series of the Plio-Pleistocene. Paleoceanography 24:PA1204. https://doi.
org/10.1029/2007PA001568
Schneider R, Schmitt J, Köhler P, Joos F, Fischer H (2013) A reconstruction of
atmospheric carbon dioxide and its stable carbon isotopic composition fromthe penultimate glacial maximum to the last glacial inception. Clim Past 9:
2507–2523. https://doi.org/10.5194/cp-9-2507-2013
Schrag DP, Berner RA, Hoffman PF, Halverson GP (2002) On the initiation of asnowball Earth. Geochem Geophys Geosyst 3(6). https://doi.org/10.1029/
2001GC000219
Seki O, Foster GL, Schmidt DN, Mackensen A, Kawamura K, Pancost RD (2010)
Alkenone and boron-based Pliocene pCO2 records. Earth Planet Sci Lett
292(2010):201–211. https://doi.org/10.1016/j.epsl.2010.01.037
Shackleton NJ (2000) The 100,000-year ice-age cycle identified and found to lag
temperature, carbon dioxide, and orbital eccentricity. Science 289:1897–1902
Shackleton NJ, Opdyke ND (1977) Oxygen isotope and paleomagnetic evidence
for early northern hemisphere glaciation. Nature 270:216–219
Stap LB, de Boer B, Ziegler M, Bintanja R, Lourens LJ, van de Wal RSW (2016) CO2
over the past 5 million years: continuous simulation and new δ11B-based
proxy data. Earth Planet Sci Lett 439(2016):1–10. https://doi.org/10.1016/j.epsl.
2016.01.022
Stap LB, van de Wal RSW, de Boer B, Bintanja R, Lourens LJ (2014) Interaction ofice sheets and climate during the past 800 000 years. Clim Past 10:2135–
2152. https://doi.org/10.5194/cp-10-2135-2014
Sun J, Huang X (2006) Half-precessional cycles recorded in Chinese loess:response to low-latitude insolation forcing during the last interglaciation.
Quat Sci Rev 25(9–10):1065–1072. https://doi.org/10.1016/j.quascirev.2005.
08.004
Svensen H, Planke S, Malthe-Sørenssen A, Jamtveit B, Myklebust R, Eidem TR, Rey
SS (2004) Release of methane from a volcanic basin as a mechanism for
initial Eocene global warming. Nature 429:542–545
Tachikawa K, Timmermann A, Vidal L, Sonzogni C, Timm OE (2014) CO2 radiative
forcing and intertropical convergence zone influences on western Pacificwarm pool climate over the past 400 ka. Quat Sci Rev 86(2014):24–34.
https://doi.org/10.1016/j.quascirev.2013.12.018
Tarnocai C, Canadell JG, Schuur EAG, Kuhry P, Mazhitova G, Zimov S (2009)
Soil organic carbon pools in the northern circumpolar permafrost region.
Global Biogeochem Cycles 23:GB2023. https://doi.org/10.1029/
2008GB003327
Tziperman E, Raymo ME, Huybers PJ, Wunsch C (2006) Consequences of pacing
the Pleistocene 100 kyr ice ages by nonlinear phase locking to Milankovitch
forcing. Paleoceanography 21(PA4206):1–11. https://doi.org/10.1029/
2005PA001241

Viaggi Progress in Earth and Planetary Science(2018) 5:81

Van de Wal RSW, de Boer B, Lourens LJ, Kohler P, Bintanja R (2011)
Reconstruction of a continuous high-resolution CO2 record over the past 20
million years. Clim Past 7:1459–1469. www.clim-past.net/7/1459/2011/.
https://doi.org/10.5194/cp-7-1459-2011
Van Der Meer D.G., Zeebe R.E., van Hinsbergen D.J.J., Sluijs A., Spakman W. &Torsvik T.H. (2014) - Plate tectonic controls on atmospheric CO2 levels since
the Triassic. PNAS Early Edition, www.pnas.org/cgi/doi/10.1073/pnas.
1315657111
Vautard R, Ghil M (1989) Singular spectrum analysis in nonlinear dynamics, withapplications to paleoclimatic time series. Physica D 35(1989):395–424 NorthHolland, Amsterdam
Walker JCG, Hays PB (1981) A negative feedback mechanism for the long-term
stabilization of Earth’s surface temperature. J Geophys Res 86(C10):9776–9782
Westerhold T, Röhl U, Frederichs T, Agnini C, Raffi I, Zachos JC, Roy H, Wilkens RH(2017) Astronomical calibration of the Ypresian timescale: implications forseafloor spreading rates and the chaotic behavior of the solar system? Clim
Past 13:1129–1152 https://doi.org/10.5194/cp-13-1129-2017
Williams DM, Kasting JF, Frakes LA (1998) Low-latitude glaciations and rapidchanges in the Earth’s obliquity explained by obliquity–oblateness feedback.
Nature 396:453–455
Wunsch C (2004) Quantitative estimate of the Milankovitch-forced contribution
to observed Quaternary climate change. Quat Sci Rev 23(2004):1001–1012.
https://doi.org/10.1016/j.quascirev.2004.02.014
Zachos JC, Dickens GR, Zeebe RE (2008) An early Cenozoic perspective on
greenhouse warming and carbon-cycle dynamics. NATURE, Vol. 451,17
January 2008. https://doi.org/10.1038/nature06588
Zachos JC, Pagani M, Sloan L, Thomas E, Billups K (2001) Trends, rhythms, andaberrations in global climate 65 Myr to present. Science 292, pp.686-693.
Zagwijn WH (1974) The Pliocene--Pleistocene boundary in western and southern
Europe. BOREAS. Vol. 3:75–97
Zhang Y.G., Pagani M., Liu Z., Bohaty S.M. & DeConto R. (2013) - A 40-million-year
history of atmospheric CO2. PHIL. TRANS. R. SOC. A, 371: 20130096. https://doi.org/10.1098/rsta.2013.0096

Page 37 of 37