the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Bispectra of climate cycles show how ice ages are fuelled
Diederik Liebrand
Anouk T. M. de Bakker
The increasingly nonlinear response of the climate–cryosphere system to insolation forcing during the Pliocene and Pleistocene, as recorded in benthic foraminiferal stable oxygen isotope ratios (δ^{18}O), is marked by a distinct evolution in iceage cycle frequency, amplitude, phase, and geometry. To date, very few studies have thoroughly investigated the nonsinusoidal shape of these climate cycles, leaving precious information unused to further unravel the complex dynamics of the Earth's system. Here, we present higherorder spectral analyses of the LR04 δ^{18}O stack that describe coupling and energy exchanges among astronomically paced climate cycles. These advanced bispectral computations show how energy is passed from precessionpaced to obliquitypaced climate cycles during the Early Pleistocene (from ∼2500 to ∼750 ka) and ultimately to eccentricitypaced climate cycles during the Middle and Late Pleistocene (from ∼750 ka onward). They also show how energy is transferred among many periodicities that have no primary astronomical origin. We hypothesise that the change of obliquitypaced climate cycles during the midPleistocene transition (from ∼1200 to ∼600 ka), from being a net sink into a net source of energy, is indicative of the passing of a landice mass loading threshold in the Northern Hemisphere (NH), after which cycles of crustal depression and rebound started to resonate with the ∼110 kyr eccentricity modulation of precession. However, precessionpaced climate cycles remain persistent energy providers throughout the Late Pliocene and Pleistocene, which is supportive of a dominant and continuous fuelling of the NH ice ages by insolation in the (sub)tropical zones, and the control it exerts on meridional heat and moisture transport through atmospheric and oceanic circulation.
 Article
(9120 KB) 
Supplement
(1370 KB)  BibTeX
 EndNote
The recurrent ice ages of the Pliocene and Pleistocene, as captured in benthic foraminiferal δ^{18}O records, are characterised by longterm trends in glacial–interglacial cycle duration, amplitude, response time, and geometry (Figs. 1 and 2). These trends mainly reflect the increasingly nonlinear response of the (northern) highlatitude cryosphere and global deepsea temperatures to radiative forcing, i.e. the combined greenhouse effect of incoming solar radiation (insolation) and the partial pressure of atmospheric carbon dioxide (${P}_{{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}}$), when global climatic conditions deteriorated (Lisiecki and Raymo, 2005a, 2007; MartinezBoti et al., 2015; Chalk et al., 2017). Spectral analysis, in combination with age control independent from astronomical tuning, formed the decisive evidence in support of Milankovitch's theory of astronomical climate forcing (Broecker, 1966; Hays et al., 1976). These statistical methods, such as the (timeevolutive) fast Fourier transform or wavelet analysis using a Morlet transform, perform well in defining sinusoidal cycle properties, such as frequency, amplitude, and/or (cross)phase relationships with respect to an astronomical reference curve (Lourens and Hilgen, 1997; Lisiecki, 2010; Meyers and Hinnov, 2010). However, spectral analysis methods are not well suited to describe highly nonsinusoidal signals, such as the skewed and asymmetric cycles that characterise time series of Middle and Late Pleistocene climate (Fig. 2) (Hagelberg et al., 1991; King, 1996; Lisiecki and Raymo, 2007). To mitigate the shortcomings of statistical tests that implicitly assume sinusoidality, higherorder spectral analysis techniques were introduced into the research field of palaeoceanography/climatology in the 1990s (Hagelberg et al., 1991, 1994; King Hagelberg and Cole, 1995; King, 1996).
Bispectral analysis was conceived in the 1960s within the research field that studies ocean waves (Hasselmann et al., 1963). It is an accepted method of quantifying nonlinear energy transfers among nearshore waves that reach breaking point when the seafloor shallows (Elgar and Guza, 1985; Doering and Bowen, 1995; Herbers et al., 2000; de Bakker et al., 2015). Since the pioneering interdisciplinary studies by palaeoceanographer/climatologist Teresa (Terri) King Hagelberg and colleagues (notably the physical oceanographer/nearshore ocean wave researcher Steve Elgar) on Pleistocene and Holocene records (Hagelberg et al., 1991, 1994; King Hagelberg and Cole, 1995; King, 1996), relatively few studies have applied higherorder spectral analysis methods to the palaeoclimate archive (e.g. Muller and MacDonald, 1997a, b; von Dobeneck and Schmieder, 1999; Rial and Anaclerio, 2000; Rutherford and D'Hondt, 2000; Huybers and Wunsch, 2004; Huybers and Curry, 2006; Liebrand et al., 2017; Da Silva et al., 2018). However, two more recent developments make a reappreciation for the potential of bispectra for palaeoclimate science timely. First, advancements have been made in the understanding and interpretation of the bispectrum (Herbers et al., 2000; de Bakker et al., 2015). These constitute (i) a shift of focus from bicoherence, which quantifies the strength of the couplings between the frequencies that are present in both the real and the imaginary parts of the complexvalued bispectrum, to (just) the imaginary part of the bispectrum, which can be used to compute nonconservative relative energy exchanges, if time series are dominated by asymmetric wave forms/cycle shapes, and (ii) integration over the imaginary part of the bispectrum to quantify conservative net energy transfers, and potentially absolute energy transfers, if net transfers can be scaled to the power spectrum (de Bakker et al., 2015, 2016). Second, increasingly noisefree benthic foraminiferal δ^{18}O stacks, characterised by precise and accurate age models, have become available to the palaeoclimate community (Lisiecki and Raymo, 2005a; Ahn et al., 2017). Such data are a prerequisite for successful application of advanced higherorder spectral analysis methods, because they describe the distribution of nonsinusoidality (i.e. relatively small amounts of signal compared to the sinusoidal cycle properties) over time and frequency (King, 1996).
Both the 40 kyr problem of the Pliocene and Early Pleistocene, and the ∼110 kyr problem (a.k.a. the 100 kyr problem) of the Middle and Late Pleistocene, are defined by considerable mismatches in spectral power, at their designated periodicities, between benthic foraminiferal δ^{18}O records (dominant) and summer insolation at any particular latitude (weaker or absent, respectively) (Raymo and Nisancioglu, 2003; Raymo et al., 2006; Lisiecki, 2010). The midPleistocene transition (MPT) (from ∼1200 to ∼600 ka) constitutes the temporal link between the 40 and ∼110 kyr “worlds”, when climate cycles in benthic foraminiferal δ^{18}O become longer in duration, of higher amplitude, and asymmetric in shape (Hagelberg et al., 1991; King, 1996; Lisiecki and Raymo, 2007). In association with this increasingly nonlinear response of Earth's climate–cryosphere system to radiative forcing, several climate cycles without a straightforward astronomical origin have been identified, such as those with periodicities of semiprecession, ∼28, ∼56, and ∼80 kyr (Rutherford and D'Hondt, 2000; Berger et al., 2006; Lourens et al., 2010). Many, often contrasting, hypotheses have been proposed to address this complex climatic evolution of astronomically forced Pliocene and Pleistocene climate cycles, which range from the merging of the Cordilleran and the Laurentide ice sheets (Bintanja and van de Wal, 2008), a delayed isostatic rebound of the lithosphere–asthenosphere after landice mass loading (AbeOuchi et al., 2013), to regolith erosion (Clark and Pollard, 1998), and/or ${P}_{{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}}$ thresholds (Chalk et al., 2017), to name just a few (Roe and Allen, 1999; Huybers, 2011). However, more information needs to be extracted from geological archives to distinguish between these mechanisms. We examine the nonlinear response of the Pliocene–Pleistocene climate–cryosphere system to radiative forcing by applying stateoftheart bispectral analysis techniques (Herbers et al., 2000; de Bakker et al., 2015, 2016) to one of the most noisefree palaeoclimate records, namely the LR04 benthic foraminiferal δ^{18}O stack (Lisiecki and Raymo, 2005a), with the aim to shed new light on (i) the nonlinear origins of climate cycles with and without a direct astronomical connection, (ii) the 40 kyr problem of the Pliocene and Early Pleistocene, (iii) the midPleistocene transition, and (iv) the ∼110 kyr problem of the Middle and Late Pleistocene. We show in great detail through which frequency interactions energy is transferred to the 40 kyr, ∼80 kyr, and ultimately ∼110 kyr cycles. These new insights can be used to better link climate cycle geometries to nonlinear response mechanisms and to obtain a better understanding of Earth's energy balance on astronomical timescales.
2.1 Benthic foraminiferal δ^{18}O record of the Pliocene and Pleistocene
To quantify nonlinear energy transfers among astronomically paced cycles of Earth's climate–cryosphere system, we use the LR04 compilation of globally distributed records of stable oxygen isotope ratios (δ^{18}O) measured on the calcite of benthic foraminifera (Lisiecki and Raymo, 2005a). The δ^{18}O values are given in per mille (‰) against Vienna Pee Dee belemnite (VPDB). This stack spans the Pliocene–Pleistocene time interval, from 5333 to 0 ka (Lisiecki and Raymo, 2005a). We recompiled the original LR04 data set (39 473 data points, http://lorrainelisiecki.com/stack.html, last access: 3 January 2019) and resampled it at 1 kyr resolution using SiZer, a method that extracts structures in time series by assessing the statistically significant zero crossing of its derivatives (Chaudhuri and Marron, 1999), because bispectral analysis requires a constant sampling resolution (Figs. 1 and 2). For the purposes of this study, we selected the smooth (out of 41) that preserves most of the structure in the data, given the 1 kyr resampling resolution. This resampling resolution resulted in an oversampling of the original data for ∼0.83 % of the record (N=328), most of which falls in the earliest Early Pliocene part of the record, which is excluded from the bispectral analysis for ages > 5250 ka. The resultant resampled LR04 record is very similar in structure to the original LR04 stack (Lisiecki and Raymo, 2005a), and, because of this similarity, we hereafter refer to the SiZer resampled LR04 record more simply as the “LR04 stack”. Prior to bispectral analysis, we detrended the LR04 stack using a Gaussian notch filter (frequency of 0.0 Myr^{−1}, bandwidth of 2.0 Myr^{−1}) to remove periodicities ≥500 kyr (Paillard et al., 1996). We also multiplied the age scale of the LR04 stack by −1, because the direction of time, which determines the sign of asymmetry values and the direction of energy transfers in bispectral analysis, increases toward the present. However, for the figures, we follow the convention and plot ages that increase with geological age.
2.2 Quantifying geometries using central moments
Skewness and asymmetry are quantified using both the third central moment and bispectral methods (Figs. 2, S2, and S3 in the Supplement) (see Sect. 2.3.3 for bispectral method) (Elgar, 1987). Excess kurtosis is quantified using the fourth central moment only, because no trispectra were calculated (Fig. 2). Using thirdmoment quantities, skewness is determined by Eq. (1):
where the overbar indicates the mean value and where 〈〉 is the timeaveraging operator (i.e. window length over which the computation is performed) (Doering and Bowen, 1995), and asymmetry is determined following Eq. (2):
where H is the Hilbert transform (Kennedy et al., 2000). Skewness is a measure for the dissymmetry of a cycle relatively to a horizontal axis (Fig. 2b), whereas asymmetry is a measure for the dissymmetry of a cycle through time (Fig. 2c). Using fourthmoment quantities, excess kurtosis is defined by Eq. (3):
Excess kurtosis is a measure for the flatness or peakedness of the extrema of a cycle. Cycles with flat tops and bottoms (i.e. platykurtic cycles), such as a square wave, have negative excess kurtosis values, whereas sharptopped and sharpbottomed cycles (i.e. leptokurtic cycles), such as a cardiogram, have positive kurtosis value (Fig. 2d).
2.3 Bispectral analysis
2.3.1 The bispectrum
In contrast to spectral analysis (e.g. using the fast Fourier transform), which gives the distribution of variance of a sinusoidal signal with frequency, bispectral analysis describes the distribution of nonsinusoidality with frequency (King, 1996). The skewness of a cycle is described by the real part of the bispectrum, whereas the asymmetry of a cycle is deconvolved in the imaginary part of the bispectrum. The bispectrum shows nonconservative relative energy exchanges among frequencies of a single time series. Energy transfers in the bispectrum computed on the LR04 stack are expressed in ‰^{3} kyr^{2}, which constitutes a statistical measure of energy. (We refer to Sect. 2.4.1 for an explanation of conservative net energy exchanges and Sect. 4.2 for a discussion of the scaling of conservative to absolute energy exchanges). Energy transfers can only occur, and be analysed by the bispectrum, if both the frequencies (f) and phases (Phi, φ, in radians) of these nonsinusoidal cycles are coupled: i.e. ${f}_{\mathrm{1}}+{f}_{\mathrm{2}}={f}_{\mathrm{3}}$ and ${\mathit{\phi}}_{\mathrm{1}}+{\mathit{\phi}}_{\mathrm{2}}={\mathit{\phi}}_{\mathrm{3}}$, respectively (Hagelberg et al., 1991). Thus, for any possible combination of three frequencies, the bispectrum assesses whether there is a coupling, and if so, whether energy exchanges occur. The transfer of energy in these socalled triad interactions is nonlinear, because the changes in cycle amplitudes (A) do not sum and have a currently unknown, and probably variable, coupling through time (i.e. ${A}_{\mathrm{1}}+{A}_{\mathrm{2}}\ne {A}_{\mathrm{3}}$; see Sect. 4.2). We note that the meaning of “cycle”, “frequency” and/or “periodicity” is different when referring to either the time or frequency domain. For example, a “single” cycle in a time series is often composed of multiple frequency components, which can be deconvolved by (bi)spectral analysis. However, for textual purposes, we use these words interchangeably, regardless of their reference to a specific domain.
The bispectrum is defined by Eq. (4):
where E[…] is the ensemble average of the triple product of complex Fourier coefficients C at the difference frequencies f_{1}, f_{2}, and their sumf_{3} (i.e. $={f}_{\mathrm{1}}+{f}_{\mathrm{2}}$), and the asterisk indicates complex conjugation (Hasselmann et al., 1963). In this study, we focus on the imaginary part of the bispectrum (hereafter often referred to as “bispectrum”, unless indicated otherwise), following studies on ocean waves (e.g. Herbers et al., 2000; de Bakker et al., 2015), because according to bispectral theory most energy transfers are associated with asymmetric cycle shapes (i.e. wave forms), and because the climate cycles of the Middle and Late Pleistocene have high amplitudes (Fig. 1) and are highly asymmetric (Fig. 2). Furthermore, conservative net energy exchanges between the coupled frequencies that are located in the imaginary part of the bispectrum can be computed if we assume a simple coupling coefficient between frequencies. Prior to bispectral analysis, we detrend the data and subsequently apply a tapering function using a Hann (a.k.a. Hanning) window that also mitigates spectral leakage (Hagelberg et al., 1991) and then multiply the data by an energy correction factor to adjust for the change in amplitude that results from the windowing.
2.3.2 Interpreting the bispectrum
The x and y axes in the bispectrum, which correspond to f_{1} and f_{2}, respectively, are mirror images of each other and share a symmetry axis at x=y (Fig. 3) (Hasselmann et al., 1963). Therefore, only the positive, oneeighth part of the bispectrum is depicted, which is subdivided into 15 zones (see Sect. 2.4.2, Table A1). For triad interactions, only two outcomes exist: either two frequencies (by definition, f_{1} and f_{2}) gain energy (though not necessarily in equal measures), and one frequency (i.e. f_{3}) loses energy, or f_{1} and f_{2} lose energy (ditto), and f_{3} gains energy. Frequencies f_{1} and f_{2} are the socalled difference frequencies, while frequency f_{3} is referred to as a sum frequency. Frequencies often participate in multiple interactions, and depending on the interaction, they act as either a difference frequency or a sum frequency. A particular frequency can thus receive energy through one interaction and simultaneously lose energy in another interaction.
To help the interpretation of the bispectrum, we depict the main astronomical frequencies in the bispectrum with coloured lines (Fig. 3). Vertical and horizontal lines correspond to the difference frequencies (f_{1} and f_{2}, respectively), and diagonal lines correspond to the sum frequency (f_{3}). Junctions between lines highlight the locations in the bispectrum where interactions occur between three periodicities (i.e. triad interactions) of which at least one is an astronomically paced climate cycle, which can also interact with itself (Fig. 3; see also Table A2). No frequency axis is associated with sum frequency f_{3}, but its value can be read off at the crossing points with the x and y axes of the bispectrum or by summing the x and y coordinates of the difference frequency axes at any point along the diagonals (Fig. 3).
We follow the convention and define a triad interaction as negative or positive if f_{3} either loses (blue colours) or gains (red colours) energy, respectively (Fig. 4). In written form, energy gains at a particular frequency are marked by an upwardpointing arrow (↑) and losses by a downwardpointing arrow (↓). Bispectral notation of triad interactions are given as B_{f}(f_{1}, f_{2}, f_{3}) in the frequency domain (Myr^{−1}) and as B_{p}(p_{1}, p_{2}, p_{3}) in the periodicity domain (kyr). For palaeoclimatological purposes, we will mainly focus on the periodicities for which the corresponding frequencies can easily be looked up in Table A2. We refer to the relevant part of the bispectrum as “Re” (i.e. real) and “Im” (i.e. imaginary), which are given in superscript (B^{Re} and B^{Im}, respectively). Thus, if we consider a negative triad interaction, located in the imaginary part of the bispectrum, between p_{1}, a 40 kyr obliquitypaced cycle (i.e. f_{1}=25.0 Myr^{−1}), and p_{2}, a ∼95 kyr eccentricitypaced cycle (i.e. f_{2}=10.5 Myr^{−1}), that results in a loss of energy at p_{3}, the ∼28 kyr periodicity (i.e. f_{3}=35.5 Myr^{−1}), this is given as ${B}_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 95↑, 28↓) (i.e. ${B}_{\mathrm{f}}^{\mathrm{Im}}$(25.0↑, 10.5↑, 35.5↓)) (Fig. 4a).
We note that within this study, we almost exclusively focus on the imaginary part of the bispectrum. This approach contrasts with previous palaeoclimate investigations that used the bicoherence spectrum (see, e.g. Hagelberg et al., 1991; von Dobeneck and Schmieder, 1999; Wara et al., 2000; Huybers and Curry, 2006; Da Silva et al., 2018). In short, bicoherence is a measure of the coherency between the real and imaginary parts of the bispectrum (i.e. where both parts have strong interactions) and is marked by positive values only (see, e.g. Elgar and Sebert, 1989; de Bakker et al., 2014). Conversely, the real and imaginary parts of the bispectrum transform skewed and asymmetric cycle geometries into their composite frequency components, respectively, and both parts are characterised by positive as well as negative interactions.
2.3.3 Quantifying geometries using the bispectrum
Adopting the bispectral method to extract cycle geometries (Figs. 2, S2, and S3) (Elgar, 1987), skewness and asymmetry are computed from the biphase. The biphase reflects the ratio between the real and imaginary parts of the bispectrum. The biphase is defined by Eq. (5):
where n and l range from 1 to the Nyquist frequency N, with n > l and $n+l\le N$ (Elgar, 1987), and f_{p} refers to the primary frequency (a.k.a. natural or resonance frequency, or first harmonic). We note that positive asymmetry values correspond to a dominance of negative interactions (i.e. a loss of energy at f_{3}) in the imaginary part of the bispectrum, and vice versa, negative values for this geometry are indicative of an overall dominance of positive interactions (i.e. a gain of energy at f_{3}).
2.4 Power spectral density and net energy transfers
2.4.1 Integration of the spectrum and bispectrum
To quantify power spectral density, we apply a standard integration of the power spectrum. To obtain conservative net energy transfers per frequency (defined by the nonlinear source term S_{nl}) and hence to extract the total gain or loss in energy per climate cycle over a specific time interval (i.e. window), we integrate over the imaginary part of the bispectrum and multiply it by a coupling coefficient. For palaeoclimate time series, no coupling coefficient has been determined previously. Therefore, we make minimum assumptions and use a coupling coefficient that only corrects for frequency ${W}_{({f}_{\mathrm{1}},{f}_{\mathrm{2}})}=({f}_{\mathrm{1}}+\phantom{\rule{0.125em}{0ex}}{f}_{\mathrm{2}})$. A comparable correction for frequency is also part of the Boussinesq scaling that is used for ocean waves (e.g. Herbers and Burton, 1997; Herbers et al., 2000). The correction for frequency (${W}_{({f}_{\mathrm{1}},{f}_{\mathrm{2}})}$) insures energy conservation during triad interactions, because more energy is exchanged among the higher frequencies. Furthermore, this correction allows for qualitative interpretations of conservative net energy exchanges across the spectrum. If a more appropriate coupling coefficient for palaeoclimatic purposes were to be deduced, it could also scale the conservative net energy transfers of the bispectrum (already corrected using S_{nl} during integration) to absolute energy transfers that can be compared to (changes in) power spectral density (see Sect. 4.2). We express the integral of S_{nl} in terms of an integration, or summation, over the positive quadrant of the bispectrum alone or equivalently in sum and difference interactions following Eq. (6):
where the sum contributions are expressed by Eq. (7):
and the difference contributions are expressed by Eq. (8):
The sum contributions are obtained by integrating diagonally over the bispectrum, and the difference contributions are obtained by integrating along vertical and horizontal lines, perpendicular to the x and y axes, respectively (Fig. 3). We track the evolution of net energy transfers through time by integrating over bispectra that are determined for consecutive and partially overlapping windows of the time series (i.e. a moving window) (Figs. 5, 6, and S3). Integration can be performed including all frequencies (i.e. over the entire imaginary part of the bispectrum) or over specific zones only, in the event that subsets of interactions need to be further examined (Fig. 3, Table A1) (de Bakker et al., 2015, 2016). After integration over the bispectrum, either totally or zonally, energy transfers computed on the LR04 stack are expressed in ‰^{3}. Blue colours indicate a loss of energy at a specific frequency; red colours indicate a gain.
2.4.2 Spectral and bispectral zones
To obtain more insight on the energy that is exchanged among eccentricity, obliquity, and precession cycles, we integrate the spectra and bispectra over separate zones (Fig. 3, Table A1). For bispectra, such a zonation approach was first applied in research concerned with nearshore ocean waves to distinguish between infragravity wave and seaswell wave frequencies (de Bakker et al., 2015), and is adapted here to investigate the most important climate cycle bandwidths of the Pliocene and Pleistocene. The boundaries between the climate cycle zones are (arbitrarily) defined at frequencies of 17.8, 35.6, and 80.0 Myr^{−1} (i.e. periodicities of ∼56.2, ∼28.1, and 12.5 kyr). In total, 15 bispectral zones are defined (Fig. 3, Table A1), each of which represents a unique combination of three main groups of astronomically paced climate cycles and suborbital periodicities (i.e. supraorbital frequencies) that can interact and exchange energy with one other. We note that many periodicities without primary astronomical origin are included in these “eccentricity”, “obliquity”, and “precession” zones. We focus here on the results of the eight most important zones, namely those involving astronomically paced climate cycles (Figs. 7, 8, 9, S4, and S5). To better highlight the roles of all nonlinear interactions involving eccentricity, obliquity, or precession, we recombine (i.e. sum) the individual zones (Fig. 10).
3.1 Geometries of Pliocene and Pleistocene climate cycles
The geometry computations confirm the qualitative visual inspection of the LR04 globally averaged, benthic foraminiferal δ^{18}O stack and show that strong nonsinusoidal cycle shapes are present in the data, especially during the Middle and Late Pleistocene parts of the record (Fig. 2) (Lisiecki and Raymo, 2007). We focus the geometry interpretations on the time interval with the greater amplitude variability in the δ^{18}O record, which starts to increase (gradually) from ∼3000 ka onward. Between ∼3000 and ∼350 ka, the skewness of climate cycles varies between −0.5 and 1.0. Peak values of ∼1.0 are reached between ∼2500 and ∼1500 ka, after which skewness rapidly decreases at ∼1500 ka to negative values of −0.5 that slowly increase to 0.0 during the MPT (Fig. 2). Glacial–interglacial cycle asymmetry varies between about −0.5 and 1.5. The most prominent steps in asymmetry occur at ∼2500 and ∼1000 ka (the latter during the MPT) when values increase from about −0.5 to 0.5 and from approximately −0.2 to 1.2, respectively. These results compare well to those previously obtained on (stacks of) benthic foraminiferal δ^{18}O records (Hagelberg et al., 1991; King, 1996; Lisiecki and Raymo, 2007). Excess kurtosis computations identify leptokurtic (thinpeaked) cycle shapes that vary between 0.5 and 3.0 from ∼3000 ka onward. Unlike skewness and asymmetry, which are marked by multiMyr decreasing and increasing trends, respectively, no longterm trend can be discerned in Pleistocene cycle kurtosis (Fig. 2).
3.2 Bispectra of Pliocene and Pleistocene climate cycles
Nonsinusoidal cycle shapes may indicate the potential for obtaining more information about nonlinear interactions among frequencies using higherorder spectral analysis (Hagelberg et al., 1991; King, 1996). We present three examples of Pliocene and Pleistocene bispectra (i.e. their imaginary parts), based on the LR04 stack, which are characterised by clear triad interactions (Fig. 4; see also Sect. 2.3.2.). The first example is a bispectrum across the Middle to Late Pleistocene (i.e. the “∼110 kyr world”), which shows clear and punctuated triad interactions concentrated where p_{2} equals 95 kyr (Fig. 4a). The main interactions are located at the following triads or their connecting lines: at $\sim {B}_{\mathrm{p}}^{\mathrm{Im}}$(95↑, 405↑, 77↓), weakly at ${B}_{\mathrm{p}}^{\mathrm{Im}}$(95↑, 125↑, 54↓), from ${B}_{\mathrm{p}}^{\mathrm{Im}}$(95↑, 95↑, 48↓) via ${B}_{\mathrm{p}}^{\mathrm{Im}}$(69↑, 95↑, 40↓) to ${B}_{\mathrm{p}}^{\mathrm{Im}}$(80↑, 80↑, 40↓), broadly around ${B}_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 95↑, 28↓), at ${B}_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 51↑, 22↓), from ${B}_{\mathrm{p}}^{\mathrm{Im}}$(31↑, 95↑, 24↓) to ${B}_{\mathrm{p}}^{\mathrm{Im}}$(29↑, 95↑, 22↓), and weakly (but still visible) at the triple junction ${B}_{\mathrm{p}}^{\mathrm{Im}}$(24↑, 95↑, 19↓). Interactions at these triads (Table A2) are all negative, indicating that energy is mainly transferred to the 95 kyr (and 125 kyr) cycle from the three precession periodicities, the ∼28 kyr periodicity, and the 40 kyr obliquity periodicity. A very weak positive interaction is located at ${B}_{\mathrm{p}}^{\mathrm{Im}}$(125↓, 405↓, 95↑). Strong energy gains at the ∼110 kyr periodicity of eccentricitymodulated precession thus characterise the bispectrum of the Middle to Late Pleistocene interval.
The second example is a bispectrum across the midPleistocene transition (i.e. the “∼80 kyr world”) and shows several negative interactions where either p_{1} or p_{2} is equal to ∼80 kyr (Fig. 4b). Particularly strong negative interactions are located at $\sim {B}_{\mathrm{p}}^{\mathrm{Im}}$(80↑, 200↑, 57↓), ${B}_{\mathrm{p}}^{\mathrm{Im}}$(80↑, 80↑, 40↓), and ${B}_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 80↑, 27↓). Also, ${B}_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 58↑, 24↓) is strongly negative. We document one major positive interaction at ${B}_{\mathrm{p}}^{\mathrm{Im}}$(69↓, 95↓, 40↑), which connects to ${B}_{\mathrm{p}}^{\mathrm{Im}}$(59↓, 125↓, 40↑). Due to the relatively loosely distributed triad interactions in this bispectrum (Fig. 4b), compared to, for example, the bispectrum of the Late Pliocene to Early Pleistocene (Fig. 4c), it is more difficult to interpret dominant energy gains and losses directly from the bispectrum of the MPT. However, we observe an overall gain of energy at the 70 to 80 kyr periodicities, and a dual role for the 40 kyr periodicity, which operates as both the energy receiver and supplier simultaneously. We note that the ∼80 kyr periodicity is not an obvious triad of the primary astronomical cycles (Table A2), apart from the junction at ${B}_{\mathrm{p}}^{\mathrm{Im}}$(95↑, 405↑, 77↓), where no strong interaction is documented in this particular bispectrum. The 80 kyr periodicity thus starts as a “lower” harmonic or “subharmonic” (i.e. where “lower” and “sub” refer to the frequency domain) of the 40 kyr periodicity, which in turn is a first harmonic of obliquity forcing, as well as a “combination tone” of precession and nonastronomical periodicities (Table A2) (Pestiaux et al., 1988). This subharmonic origin for the 80 kyr periodicity is supported by the strong negative interaction at ${B}_{\mathrm{p}}^{\mathrm{Im}}$(80↑, 80↑, 40↓).
The third example is a bispectrum across the Late Pliocene to Early Pleistocene interval (i.e. the “40 kyr world”), which depicts negative and positive interactions, where both difference frequencies f_{1} and f_{2}, and sum frequency f_{3} are equal to 25 Myr^{−1} (i.e. 40 kyr), respectively (Fig. 4c). More specifically, negative interactions are concentrated at and between triads along the lines from ${B}_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 405↑, 36↓) to ${B}_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 40↑, 20↓), and from ${B}_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 40↑, 20↓) to ${B}_{\mathrm{p}}^{\mathrm{Im}}$(36↑, 40↑, 19↓), indicating a gain at the 40 kyr periodicity. Most notable is the strong negative interaction at $\sim {B}_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 93↑, 28↓) to $\sim {B}_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 83↑, 27↓). Other bispectral coordinates marked by negative interactions are those between ${B}_{\mathrm{p}}^{\mathrm{Im}}$(95↑, 95↑, 48↓) and $\sim {B}_{\mathrm{p}}^{\mathrm{Im}}$(70↑, 150↑, 48↓). Positive interactions are concentrated on the line between the bispectral coordinates ${B}_{\mathrm{p}}^{\mathrm{Im}}$(80↓, 80↓, 40↑) and ${B}_{\mathrm{p}}^{\mathrm{Im}}$(44↓, 405↓, 40↑), which also indicate a gain at the 40 kyr periodicity by transfers from numerous other frequencies. Overall, the Late Pliocene to Early Pleistocene bispectrum depicts strong energy gains at the obliquitypaced 40 kyr periodicity through many triad interactions, some of which are with other astronomically paced climate cycles.
3.3 Qualitative comparison of summerhalf insolation to the climate–cryosphere proxy record
Figure 5 shows three timeevolutive spectral and bispectral analyses for the Pliocene–Pleistocene time interval. It compares the total power spectral density of summerhalf insolation at 65^{∘} N (Fig. 5a) (Laskar et al., 2004) to the total energy transfers and total power spectral density of the LR04 benthic foraminiferal δ^{18}O stack (Fig. 5b and c, respectively). This insolation curve constitutes a classical shorthand for the complex external forcing of the ice ages (not shown, but see Figs. 1e and 2e for individual records of precession, obliquity and eccentricity) (Köppen and Wegener, 1924; Milankovitch, 1941). Figure 6 highlights five equally spaced (500 kyr apart) time slices of the three panels of Fig. 5. Figure 7 depicts three zonally averaged frequency slices of the three panels of Fig. 5, namely the “eccentricity”, “obliquity”, and “precession” zones (see Sect. 2.4.2.), as well as the total of all frequencies combined.
3.3.1 Power spectral density of insolation
The timeevolutive power spectra of insolation are marked by high spectral power at the precession and obliquity periodicities, and the near absence of spectral power at the eccentricity periodicities (Figs. 5a, 6, and 7a) (Hays et al., 1976; Berger, 1977). Furthermore, it is characterised by amplitude variability in spectral power at these periodicities, which is governed by the longterm (eccentricity) modulations of obliquity and precession. These amplitude modulations have durations of ∼180 kyr (not detected here due to window length) and 1.2, 2.4, and 4.8 Myr for obliquity, and 405 kyr and 2.4 Myr for eccentricitymodulated precession (Figs. 5a and 7a). Concurrent with the MPT, we document the sharpest decrease in total spectral power of the entire Pliocene and Pleistocene, which is largely due to a decrease in precession power from maximum values (Fig. 7a) (Laskar et al., 2004). This is indicative of a reduction in the perturbation of the Earth's system that results from a smaller difference in insolation forcing between the hemispheres (i.e. reduced seasonality).
3.3.2 Energy transfers among climate cycles
Figures 5b and S3 depict the conservative net energy transfers among (astronomically forced) climate cycles as present in the LR04 stack. The net energy transfers show how energy is redistributed over the power spectrum among astronomically paced and nonastronomically paced climate cycles. They reveal strong increases in both energy gains (in red) and losses (in blue) toward the present, and from ∼1000 ka onward several consecutive shifts in energy transfers from lower to higher periodicities occur (Fig. 5b). Energy gains are highly concentrated at the 40 kyr periodicity from ∼2500 ka onward and are marked by a modest 1.2 Myr beat in amplitude (Figs. 5b and 7b). A second location of strong energy gains is the ∼80 kyr periodicity from ∼1000 ka onward (Fig. 6b), which at ∼700 ka shifts (Fig. 5b) to a dominant gain for the 95 kyr component (Fig. 6a). Further energy gains are located at a ∼200 kyr periodicity from ∼750 ka onward, a 50 to 60 kyr periodicity between ∼1500 and 500 ka, and (very weak) at the 24 kyr periodicity from ∼650 ka onward (Fig. 5b). Energy losses are also observed in the total integration over the bispectra (Figs. 5b, 6, and 7b). However, these losses are generally of smaller amplitude per cycle and are distributed across a greater number of periodicities, in comparison to the highly concentrated energy gains (Figs. 5b and 6). Despite this more dispersed pattern of energy losses, we highlight the more conspicuous decreases: the 40 kyr periodicity from ∼550 ka onward, the 30 to 36 kyr periodicities between 1500 and 1000 ka, the 30 to 26 kyr periodicities between ∼2500 and ∼2000 ka, the ∼28 kyr periodicity from ∼800 ka onward, the 24 kyr periodicity between ∼1000 and ∼800 ka, the 22 kyr periodicity from ∼700 onward, the 19 kyr periodicity from 2000 ka onward, and the 15 kyr periodicity intermittently from ∼1000 ka onward (Fig. 5b).
3.3.3 Power spectral density of climate cycles
Timeevolutive power spectral analysis of the LR04 record is characterised by highamplitude variability at the 40 kyr obliquitypaced cycle, from approximately 4000 ka onward (Figs. 5c, 6, and 7c). In addition, distinct ∼80 and ∼110 kyr cycles of rapidly increasing amplitude are documented from ∼1200 ka onward, marking the MPT (Fig. 6b) and postMPT (Fig. 6a) intervals, respectively (Lisiecki and Raymo, 2005a, 2007). We note that the 125 kyr component of the ∼110 kyr signal in the power spectra (Fig. 5c) is not present in the integration over the bispectrum (Fig. 5b) A relatively weak ∼200 kyr cycle is present from ∼750 ka onward (Fig. 5c). No strong precession response is present (Raymo and Nisancioglu, 2003; Raymo et al., 2006). The 40 kyr cycle in δ^{18}O partially mimics the 1.2 Myr amplitude modulation of the obliquity cycle, especially from ∼3000 ka onward (Fig. 5c) (Lourens and Hilgen, 1997; Lisiecki and Raymo, 2007). We draw attention to a subtle “arc” in spectral power (absent in bispectral power), which starts with a strong 405 kyr cycle, at ∼3500 ka, that becomes gradually shorter in duration until it resonates with the 95 kyr eccentricity modulation of precession, at ∼2650 ka, and then culminates in several weaker 70 to 80 kyr cycles at ∼2000 ka (Fig. 5c; see also Fig. 1b) (Lisiecki, 2010; Meyers and Hinnov, 2010; Viaggi, 2018).
3.4 Zooming in on energy transfers
To better understand the different roles of specific (bandwidths of) climate cycles in nonlinear interactions, we determine conservative net energy transfers per bispectral zone (Figs. 8, S4, and S5). These zones reflect the “layers” that make up total energy exchanges as depicted in Fig. 5b. We assess if triad interactions are conservative of energy for each of the bispectral zones and for the entire bispectrum (Fig. 9). Lastly, we recombine the bispectral zones again for frequency bandwidths involving precession, obliquity, or eccentricitypaced climate cycles to highlight their influences in nonlinear interactions (Fig. 10).
3.4.1 Zonal energy transfers
In Zone 1, from 1000 ka onward, we document a modest exchange of energy between the (eccentricityrelated) periodicities in the range from ∼200 to 50 kyr (Fig. 8a). Notably, from ∼500 ka onward, periodicities of ∼80 to 50 kyr predominantly lose power, whereas those between ∼200 and ∼80 kyr mainly gain energy. Stronger exchanges are documented in Zone 2 (Fig. 8b). Especially the 40 kyr obliquitypaced cycle gains energy between 3500 and 900 ka in this zone, after which it loses energy from 900 ka onward. Also, the 50 to 60 kyr periodicity loses energy, especially from ∼600 ka onward. Between ∼1100 and ∼800 ka, a single periodicity of ∼87 kyr gains power. At ∼750 ka, this periodicity splits into an ∼80 kyr periodicity and a 95 kyr eccentricitypaced cycle, which both gain energy. In general, interactions in Zone 2 change from (weakly) positive (B^{Im}(E↓, E↓, O↑)) to negative (B^{Im}(E↑, E↑, O↓)) during the MPT (∼1000 ka), indicative of a change from escalating energy to higher frequencies, to a cascading from higher frequencies. We note that the energy exchanged through B^{Im}(E, E, O) interactions (Zone 2) is almost of equal amplitude to that of all other interactions combined. Zones 3 to 5 are mainly marked by increases in energy at the 40 kyr cycle, roughly from 2500 ka onward (Fig. 8c to e). A smaller gain is documented in these zones (especially in Zone 3) at periodicities between ∼70 and 50 kyr from ∼1300 ka onward, and in Zone 3 at the 95 kyr cycle from ∼650 ka onward. Furthermore, we document decreases in energy at the ∼30 to 27 kyr periodicities, the ∼28 to 22 kyr periodicities, and the main precession periodicities of 24, 22, and 19 kyr, for Zones 3, 4, and 5, respectively. The switchover during the MPT, at about 750 ka, of energy loss at the 24 kyr cycle to a loss at the ∼28 kyr periodicity is associated with B^{Im}(O, E, P) interactions (Zone 4) only. Zones 6 and 7 show weaker energy exchanges among the main astronomical periodicities. In Zone 6, the ∼27 and 24 kyr cycles show a small gain between ∼1000 and 700 ka, and from ∼700 ka onward, respectively. We document only very minimal direct fuelling of eccentricitypaced climate cycles by precessionpaced climate cycles in this zone. Energy is lost in Zone 6 from 1500 ka onward at the 19 kyr cycle. In Zone 7, relatively weak gains are present at the main obliquity and precession cycles from 1500 ka onward. Energy losses of Zone 7 are located mainly at subprecession periodicities of ∼15 to < 13 kyr (Fig. 8f and g). No energy exchanges of similar magnitude are documented in Zone 8 (Fig. 8h).
Overall, the zonal integrations are marked by negative triad interactions and thus reveal a cascade of energy from the (sub)precession bandwidth to the obliquity and (ultimately) eccentricity bandwidths through several successive interactions (Fig. 8). These analyses further show that many periodicities fulfil a dual role that often remains stable through time, namely they serve simultaneously and persistently as energy provider and receiver. For example, the 24 kyr precessionpaced cycle gains energy through B^{Im}(P, E, P) and B^{Im}(P, O, P) interactions, as documented in Zones 6 and 7, meanwhile losing energy through B^{Im}(O, O, P) and B^{Im}(O, E, P) interactions, as can be observed in Zones 4 and 5. As an exception to this persistent behaviour of many lower periodicities in transferring energy to higher periodicities, the role of the 40 kyr obliquitypaced cycle evolves but only in B^{Im}(E, E, O) interactions (Zone 2). In Zones 3, 4, 5, and 7, the 40 kyr obliquitypaced climate cycle (predominantly) partakes in triad interactions as a difference frequency, and it continuously receives energy, mainly from precession. However, during the MPT (from ∼900 ka onward) in Zone 2, when the 40 kyr periodicity partakes in interactions solely as a sum frequency (f_{3}) of two eccentricity components (f_{1}+f_{2}), it starts providing energy to the eccentricity bandwidth. The net (i.e. total) result of these opposing roles for the 40 kyr obliquitypaced cycle after the MPT is that energy losses overtake the gains (Figs. 5b and S3). However, the initial decrease of the integrated “obliquity” zone, between ∼1000 and 600 ka, is linked to the ∼28 kyr periodicity, not the 40 kyr periodicity (Fig. 5b).
3.4.2 Energy conservation
Computations of energy conservation, hereafter also referred to as “conservativity”, indicate that both the total as well as the zonal energy gains in triad interactions are commensurate with losses, given our assumed coupling coefficient (Fig. 9). The dominance of negative interactions (i.e. positive values for f_{1}+f_{2} and negative values for f_{3}) is indicative of a transfer within the Earth's climate–cryosphere system of insolation energy to lower frequencies. Energy is not well conserved in triad interactions documented by Zone 8. This is probably the result of the low number of interactions located in Zone 8 (Fig. 4) and the small area of this zone (Fig. 3).
3.4.3 Energy exchanges involving precession, obliquity, or eccentricitypaced climate cycles
To highlight the separate roles of the precession, obliquity, and eccentricity bandwidths during nonlinear interactions within the climate–cryosphere system, we recombine (i.e. sum) their respective bispectral zones of Fig. 8. To do so, we include those zones that contain either one, two, or three frequency components of a particular bandwidth (Figs. 3 and 10). The general pattern of an energy cascade is robust, but it becomes clearer that precession cycles do not contribute much to the fuelling of eccentricitypaced climate variability directly (Fig. 10a–c). The 40 kyr periodicity during the Middle and Late Pleistocene gains energy from precessiondominated interactions but loses energy in obliquity and eccentricitydominated interactions (Fig. 10a–c). A comparison of conservativity of the recombined zones that include at least one precession, obliquity, or eccentricity component indicates that approximately similar amounts of energy are exchanged in (triad) interactions involving obliquity as in those involving eccentricity (Fig. 10d).
4.1 Limitations
4.1.1 Age uncertainty and signaltonoise ratio
By combining more than 50 individual benthic foraminiferal δ^{18}O records, the LR04 stack simultaneously decreases the age uncertainty and increases the ratio of signaltonoise considerably in comparison to any record individually (Lisiecki and Raymo, 2005a). Reduced age uncertainties and improved signaltonoise ratios of stacked records may be obtained through the computation of averaged sedimentation rates in combination with statistical “binning” (Lisiecki and Raymo, 2005a), a probabilistic approach using a hidden Markov model (Ahn et al., 2017) or, as is the case in this reanalysis of the (1 kyr resampled/binned) LR04 stack, by additional statistical “smoothing” (Chaudhuri and Marron, 1999). For successful application of advanced bispectral analysis techniques, good age control and (very) high signaltonoise ratios are needed, because the distributions of only small amounts of signal over time and frequency are considered (Hasselmann et al., 1963; King, 1996). For consecutive computational steps that use bispectra (i.e. the bispectrum and the integration over the bispectrum), increasingly better constrained and more noisefree data are required, because each step zooms in further on these small amounts of signal that are associated with nonsinusoidal geometries. The bispectra of climate cycles presented here show clear interactions (Fig. 4), and the integrations of the bispectrum reveal qualitatively consistent patterns of energy transfers for analysis with varying window lengths (Figs. 5b and S3). Therefore, we argue that the LR04 stack has sufficiently accurate age calibrations and a high enough signaltonoise ratio for the successful application of advanced bispectral analysis.
4.1.2 Resolvability in time and frequency domains
Similar to the power spectrum, the bispectrum is marked by a tradeoff between resolution in the time versus frequency domains (Fig. S3). Bispectral results become more significant, and can yield greater degrees of freedom, if a larger number of similarly shaped cycles, or waveforms, can be included for analysis (i.e. a lengthening of the window). For example, studies on natural nearshore waves or those on flumegenerated waves (periodicities of seconds to minutes) use hourlong stable time series with high sampling resolutions and wave numbers (de Bakker et al., 2014). Bispectral analyses of such data sets are well resolved in the frequency domain. Furthermore, they permit the computation of robust confidence levels on the results and the selection of statistical settings that yield high degrees of freedom (e.g. Elgar and Guza, 1985; de Bakker et al., 2015). However, for the Pliocene and Pleistocene record, similar statistical objectives are unattainable, because no two climate cycles are alike and because of the low number of ∼110 kyr cycles in the Pleistocene (King, 1996; Lisiecki, 2010). As a result of this absence of “ergodicity” (i.e. a stable waveform through space and time), which characterises the unique Pliocene–Pleistocene climatic evolution as captured by the LR04 stack, relatively short window lengths are preferred to obtain wellresolved results in the time domain. The Hann windowing function visually clarifies the interactions in both the bispectra as well as the integrations over the bispectrum. Similarly to previous bispectral studies on the palaeoclimate archive (Hagelberg et al., 1991; King, 1996), we cannot define confidence levels for the bispectral analysis of the LR04 stack, because the number of cycles per window is too small. However, the bispectral couplings among climate cycles documented here are robust for varying bispectral settings (Figs. 5, 8, and S3 to S5).
4.1.3 Quantifying cycle geometries of nonergodic records
In contrast to the relatively more robust bispectral results (see previous section), absolute geometry values are more strongly affected by the nonergodic nature of the Pliocene–Pleistocene record. Choices of (i) window length, (ii) method of detrending the time series, and (iii) the application of a windowing function or not can have profound effects on the robustness of the results (Tim E. van Peer, personal communication, 2018). These somewhat arbitrary choices affect both ways of quantifying geometries similarly (i.e. using the central moments or higherorder spectra; Figs. 2, S2, and S3) (Tim E. van Peer, personal communication, 2018). The sensitivity of geometry computations to data processing techniques was not previously appreciated in full and may have resulted in overoptimistic confidence levels for cycle geometry values computed on an Oligocene and early Miocene climate record (Liebrand et al., 2017). We note that the uncertainty in skewness and asymmetry becomes larger for cycles with values near the mean (=0). The relatively high values of about −1 or +1, similar to those computed on the LR04 stack of the Middle and Late Pleistocene, are much less sensitive to data processing methods and have been independently computed on at least three occasions (Fig. 2) (Hagelberg et al., 1991; King, 1996; Lisiecki and Raymo, 2007). Remarkably, excess kurtosis values between ∼3000 and ∼350 ka even change sign between computations that do apply a Hann window (resulting in leptokurtic cycles; Fig. 2) and those that do not (yielding platykurtic cycles that reach values of −1.0 during the Pleistocene; not shown). By comparison, skewness and asymmetry values appear much more robust to various methods of data processing than values of excess kurtosis. This greater sensitivity of excess kurtosis to windowing is not surprising considering that fourthmoment statistics (and the trispectrum) zoom in on even smaller amounts of signal in comparison to the third central moments and the bispectrum (Collis et al., 1998).
4.2 Coupling coefficient
For ocean waves, the relationship between the integration of the bispectrum and absolute nonlinear energy transfers is determined by a coupling coefficient, which is based on the secondorder Boussinesq theory (Herbers and Burton, 1997; Herbers et al., 2000). This Boussinesq approximation describes how energy is transferred from the primary ocean waves to higher and lowerfrequency components by nearresonant nonlinear triad interactions. A coupling coefficient is needed (i) to ensure the conservation of energy within a single triad interaction (i.e. all energy lost at a certain frequency is compensated by an energy gain at the other two frequencies, and vice versa), and (ii) to obtain absolute energy transfers, which are directly comparable with the changes in the power spectrum when the waves propagate toward the coast (Herbers et al., 2000; de Bakker et al., 2014). Here, we apply a similar approach to palaeoclimate time series and demonstrate that energy conservation can be ensured assuming a simple coupling coefficient that multiplies the energy transfers at each specific sum frequency after integration, with that specific frequency (f_{3}) (see Sect. 2.4.1). Comparison of energy gains and losses, for both the integration over the entire imaginary part of the bispectrum, as well as for each bispectral zone separately, shows that sum (${S}_{{\mathrm{nl}}^{+}}$; see Eq. 7) and difference (${S}_{{\mathrm{nl}}^{}}$; see Eq. 8) interactions balance each other, which is indicative of energy conservation, and qualitatively trustworthy signals in energy transfers (Fig. 9). However, the second step toward absolute energy transfers is not developed here (see qualitative comparison of spectral and bispectral results in Fig. 6), because a more fundamental understanding is needed of the many climatologic and oceanographic, biologic, sedimentologic, and lithologic processes that affect the (globally integrated) benthic foraminiferal δ^{18}O record. Ultimately, these poorly constrained processes underpinning the LR04 stack determine how to scale bispectral energy transfers (in ‰^{3}) to changes in power spectral density (in ‰^{2} kyr) in absolute terms.
5.1 Revisiting the problem of the ice ages
The common denominator of (i) the 40 kyr problem of the Pliocene and Early Pleistocene, (ii) the ∼110 kyr problem of the Middle and Late Pleistocene, and (iii) the (nonlinear) origins of many climate cycles without primary astronomical pacemaker (e.g. the ∼28 and 80 kyr periodicities) is the mismatch in the distribution of power spectral density between benthic foraminiferal δ^{18}O and records and insolation (Raymo and Nisancioglu, 2003; Lisiecki, 2010; Lourens et al., 2010). Furthermore, (iv) during the midPleistocene transition, this discrepancy increases significantly (Lisiecki and Raymo, 2005a; Clark et al., 2006; Huybers, 2007) by the shifting of spectral power of climate records from the 40 kyr, via the ∼56 and the ∼80 kyr, to the ∼110 kyr periodicity (Rial and Anaclerio, 2000; Lourens et al., 2010; Viaggi, 2018). In light of the bispectral results obtained in this study, we revisit the problem of the ice ages (Agassiz, 1840; Milankovitch, 1941; Hays et al., 1976) in search of the most parsimonious explanation(s) for the origins and evolution of climate cycles with and without astronomical periodicities. The answer that we distil from the bispectra of Pliocene and Pleistocene climate cycles, in greater clarity than before (Hagelberg et al., 1991; Rial and Anaclerio, 2000; Lourens et al., 2010), is that these problems (i.e. i to iv) are intrinsically linked by the (many) increasingly nonlinear responses of Earth's climate–cryosphere system to insolation forcing (Pisias et al., 1990; Yiou et al., 1994). Namely, the bispectra confirm (Hays et al., 1976; Hagelberg et al., 1991) that energy is transferred from the lower to the higher periodicities from at least 2500 ka onward. More importantly, they show (i) how this happens, i.e. through which nonlinear triad interactions (Fig. 4), and (ii) how their involvement in transferring energy evolves through time (Figs. 5, 8, and 10), especially during MPT.
5.2 The origins of climate cycles: toward a better mechanistic understanding
To advance the understanding of the origins of Pliocene and Pleistocene climate cycles, we explore the potential climatic and cryospheric processes that cause the documented nonlinear triad interactions and energy transfers. Similar to the difficulties in linking spectral peaks of benthic foraminiferal δ^{18}O records, which describe the integrated effect of deep water temperatures and cryosphere (Hays et al., 1976; Elderfield et al., 2012; Rohling et al., 2014; Crowhurst et al., 2018), to particular mechanisms, it is not straightforward to identify (singular) causes for bispectral peaks (Fig. 4). Should we attempt to link every individual triad interaction to a specific climatic or cryospheric process? Or, alternatively, should we search for one, or a perhaps a couple of, mechanism(s) that can explain multiple, or even all, nonlinear couplings at once? A comparison to nearshore ocean waves, and the interactions that characterise them as they break in the surf zone, may be informative (Herbers et al., 2000; de Bakker et al., 2014, 2016). The strength and number of nonlinear interactions of shoaling and breaking waves is strongly determined by the shallowing of the water depth when the wave approaches the coast and starts to interact with the sea bed (de Bakker et al., 2015, 2016). Thus, a single mechanism, i.e. the interaction between a wave and the beach, causes multiple nonlinear triad interactions.
In contrast to the single cause for multiple nonlinear triad interactions among ocean waves, we tentatively link two mechanisms to energy transfers among climate cycles, because we observe two key features in the total and zonal integrations over the bispectrum: (i) a persistent fuelling of obliquitypaced climate cycles by precessionpaced climate cycles, which we link to astronomical forcing of atmospheric and oceanic circulation (“the climatic precession motor”) (Sect. 5.2.1), and (ii) a change of obliquitypaced climate cycles from net energy sink into net energy source, which we link to a resonance of cycles of crustal sinking and rebounds with the eccentricity modulation of precession, after Northern Hemisphere (NH) landice mass loading passed a critical threshold during the MPT (“the cryospheric obliquity motor”) (Sect. 5.2.2.).
5.2.1 The climatic precession motor
First, we speculate that the fuelling of the ice ages by (mainly) the precession periodicities is largely linked to the low to midlatitude monsoons (Rutherford and D'Hondt, 2000; Berger et al., 2006; Bosmans et al., 2015a). The (sub)tropical zones, where most insolation is received, constitute Earth's heat engine, and their monsoons function as both a source of moisture and as a meridional teleconnection to build up large polar (land)ice volumes. The interactions between NH landice volumes and the African and Asian monsoons are well documented for large parts of the Pliocene and Pleistocene (deMenocal, 1995; Sun et al., 2006; Cheng et al., 2016). The monsoonal (i.e. atmospheric) redistribution of heat and moisture was most probably aided by orbitalscale variability in the strength of Atlantic (i.e. oceanic) meridional overturning circulation (AMOC) (Lisiecki et al., 2008; Ziegler et al., 2010; Elderfield et al., 2012). Based on the bispectral results, we infer that during the Pliocene and Early Pleistocene this predominantly monsoonally driven precession motor fuels the 40 kyr obliquitypaced iceage cycles, aided by more linear climatic–cryospheric responses resulting from variability in insolation at this periodicity, either zonally/crossequatorially integrated (Huybers and Tziperman, 2008; Bosmans et al., 2015b) or locally at higher latitudes. The precession motor is especially strongly expressed in B^{Im}(O, O, P) interactions (Zone 5, Fig. 8e). However, uncertainty remains over the exact moisture sources for the Pleistocene landice volume fluctuations (e.g. Werner et al., 2001; de Boer et al., 2012).
5.2.2 The cryospheric obliquity motor
Second, we hypothesise that the role reversal of obliquitypaced climate cycles (from net energy receiver to net energy provider) across the MPT, in both B^{Im}(E, E, O) interactions (Fig. 4) and total net energy transfers (Figs. 5b and 8b), is linked to a resonance of (auto)cycles of crustal sinking and a delayed rebound with ∼110 kyr long (allo)cycles of eccentricitymodulated precession. In this hypothesis, crustal sinking is triggered by the passing of a landice mass loading threshold, and the delayed response to both eccentricitymodulated precession and obliquity is caused by the inertia of bulky, continentalsized ice sheets, which also results in a nonlinear positive icealbedo feedback and hence a hysteresis pathway for glacial–interglacial cycles (Calov and Ganopolski, 2005; AbeOuchi et al., 2013). This obliquity motor evolves during and after the MPT from a lower harmonic of two to three obliquitypaced cycles (Fig. 6b) (Ridgwell et al., 1999; Huybers and Wunsch, 2005) into a resonance (i.e. phase coupling/locking or synchronisation) with the ∼110 kyr eccentricitypaced cycles, mainly the ∼95 kyr component (Fig. 6a) (Hagelberg et al., 1991; Tziperman et al., 2006; Crucifix, 2012; AbeOuchi et al., 2013). We speculate that the resonance of crustal sinking and a delayed isostatic rebound of the lithosphere–asthenosphere also shifted the sensitivity of the climate–cryosphere system to insolation forcing at higher latitudes, where the relative influence of obliquity compared to precession is greater than at lower latitudes. We conjecture that this shift in sensitivity to obliquity forcing was especially enhanced during glacial terminations (Huybers and Wunsch, 2005; Huybers, 2011; Konijnendijk et al., 2015).
5.3 Climatic and tectonic boundary conditions
If our interpretations of the bispectra are indeed supportive of an astronomically paced monsoonal/AMOC “push” of (heat and) moisture, followed by a resonant lithospheric–asthenospheric “pull” through a cooling of the NH polar region after an initial phase of ice growth, for sufficiently large glacial landice masses to form, then they would highlight the importance of delayed isostatic rebound and the shape (i.e. extent, volume, and especially height) of the NH ice sheets in determining the pattern and geometry of climate–cryosphere variability observed for the Middle and Late Pleistocene (Roe, 2006; Bintanja and van de Wal, 2008; AbeOuchi et al., 2013). However, these mechanisms do not explain why largescale NH glaciations only developed during the Late Pliocene and Pleistocene, and not during earlier time periods with comparable astronomical configurations. The Pliocene–Pleistocene climatic–cryospheric evolution can therefore not be seen outside a context of longterm, multiMyr changes in climatic and tectonic boundary conditions, not captured by the bispectral analysis presented here. For example, the “arc” in spectral power (from ∼3500 to ∼2000 ka) (Fig. 5c), which is largely absent in the bispectral results (Fig. 5b), and which precedes the strong glaciations of the MPT and Late Pleistocene, may be indicative of these gradually changing boundary conditions. We speculate that this arc reflects a change in the response time of the cryosphere, potentially through glacial landscaping by earlier (smaller) glacial cycles (Clark and Pollard, 1998; Tabor and Poulsen, 2016) or through changing tectonic/geographic boundary conditions and associated oceanic circulation patterns (Haug and Tiedemann, 1998; Kender et al., 2018). Such developments could have preconditioned the Earth's system for the fullblown bipolar glacial cycles of the Middle and Late Pleistocene by moving the NH cryosphere further into the “Milankovitch window”. These glacial–interglacial cycles were amplified by variability in radiative forcing on astronomical timescales and set against a background of an allCenozoictime low ${P}_{{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{atm}}}$ (Shackleton, 2000; Beerling and Royer, 2011; MartinezBoti et al., 2015; Chalk et al., 2017).
We present considerable advances in the understanding of the origins of many Pliocene and Pleistocene climate cycles by applying advanced bispectral analysis to the palaeoclimate archive. However, several lines of research remain to be explored further. Here, we list those that we think need more attention in the future:

For nearshore ocean waves, energy transfers are (almost) fully described by the imaginary part of the bispectrum (Norheim et al., 1998), because (i) nearshore ocean waves are dominated by asymmetric wave forms, and (ii) the imaginary part of bispectra in bispectral theory in general describes the majority of energy that is transferred. The importance of the imaginary part in describing most/all energy exchanges can easily be observed when computing power spectra of synthetically skewed and asymmetric time series (not shown here). Such spectra indicate that the higher harmonic peaks of skewed cycles are about an order of magnitude smaller than those of asymmetric cycles, especially for the higher harmonic peaks (3f, 4f, etc.) in comparison to the first higher harmonic (2f). Despite these profound implications of bispectral theory, and the high asymmetry values of Middle and Late Pleistocene climate cycles, a more comprehensive understanding of the real part of the bispectrum may still be important for palaeoclimate studies, because it may enable an even more complete description of nonlinear energy exchanges (Figs. S1 and S2; e.g. see the ∼28 kyr periodicity). However, triad interactions documented in the real part of the bispectrum are not conservative of energy.

Although we obtain conservative energy exchanges with our simplified, assumed coupling coefficient, additional studies are needed to scale these exchanges to the power spectrum (Fig. 6). Coupled climate–cryosphere models may be used to derive or approximate coupling coefficients for palaeoclimates that are based on “first principles” (i.e. a physicochemical understanding of the Earth system). Alternative coupling coefficients may also change how the bispectrum redistributes energy from climate cycles with low to those with high periodicities. However, we note that energy conservation should be respected for every triad interaction in the bispectrum.

The proposed relationships between nonlinear triad interactions as documented in bispectra of climate cycles and specific processes of the climate–cryosphere system (i.e. monsoons, isostatic rebound) remain speculative at best. Therefore, we cannot currently separate these mechanisms from, for example, deep ocean carbon storage (e.g. Willeit et al., 2015; Ganopolski et al., 2016; Farmer et al., 2019), AMOC shutdown and subsequent overshoot (Liu et al., 2009; Wu et al., 2011), and/or changing sediment cover (Clark and Pollard, 1998; Ganopolski and Calov, 2011; Willeit et al., 2019) and the amplifying effects these processes may have had on the strength of the 40 and ∼110 kyr cycles. More data (analysis) and/or (isotopeenabled) modelling studies are needed to support our preferred hypotheses. For example, global climate models could potentially help identify which mechanisms correspond to the specific frequency interactions we observe in the LR04 stack, and energy balance models may provide further insights on the amounts of energy that are transferred.

There is a limited amount of information that can be extracted from studying a single climate record. Bispectral analysis of other (Pliocene–Pleistocene) records may shed more light on the dynamics of Earth's palaeoclimate system not only in elucidating the couplings among astronomical frequencies but also across the spectral continuum, in relation to 1∕f “red” noise, and/or the powerlaw scaling of energy exchanges between frequencies (see Sects. 2.4.1. and 4.2) (Hasselmann, 1976; Bak et al., 1987, 1988; Bak, 1996; Pelletier, 1998; Huybers and Curry, 2006). Examples of interactions that could be studied in this way are those between astronomical and subastronomical (e.g. millennial) climate cycles (Hagelberg et al., 1994; von Dobeneck and Schmieder, 1999; Wara et al., 2000; Da Silva et al., 2018), and among cycles with centennial, decadal, and perhaps even annual durations (King Hagelberg and Cole, 1995). We emphasise that other applications of bispectral analysis techniques to the palaeoclimate archive should be limited to records with good age control and high signaltonoise ratios (see Sect. 4.1.1).
In this interdisciplinary study, we present a new, higherorder spectral analysis and interpretation of Pliocene and Pleistocene climate cycles as present in the wellestablished LR04 globally averaged benthic foraminiferal δ^{18}O stack. We use the nonsinusoidal properties of these climate cycles as a proxy for the strength of nonlinear couplings and interactions among different components of the climate system that operate on distinctly separate timescales. These advanced bispectral analysis techniques were initially developed by the nearshore ocean wave community and are adapted here to make them more applicable for palaeoclimatic research purposes. They show, more thoroughly than before, the complexities of the Pliocene–Pleistocene climatic–cryospheric evolution, namely that (i) nonlinear energy exchanges result in the distinct asymmetric cycle geometries of the Middle and Late Pleistocene; (ii) energy is transferred across the power spectrum, among astronomical and nonastronomical climate cycles; (iii) precessionpaced climate cycles fuel both the 40 kyr and (largely indirectly) the ∼110 kyr iceage cycles of the Pliocene and Pleistocene; (iv) the role of obliquitypaced climate cycles changes from being a net sink into a net source of energy during the MPT (i.e. in B^{Im}(E, E, O) interactions); and (v) after the MPT, these obliquitypaced climate cycles help fuel the ∼110 kyr eccentricitypaced climate cycles of the Middle and Late Pleistocene.
We postulate two distinct fuelling mechanisms for the ice ages to explain the energy cascade of negative interactions from climate cycles with high frequencies to those with low frequencies: (i) a continuous Pliocene–Pleistocene climatic “precession motor” and (ii) a MiddletoLate Pleistocene cryospheric “obliquity motor”. We interpret the dominant precession fuelling of obliquitypaced climate cycles (i.e. the precession motor) as evidence of a lowlatitude, (sub)tropically driven climate system, in which the monsoons (and oceanic meridional overturning circulation) transport heat and moisture to higher latitudes. We argue that the role reversal of obliquitypaced climate cycles during the MPT is indicative of the passing of a landice mass loading threshold, after which autocycles of crustal sinking and (delayed) rebound start to resonate, i.e. become frequency and phase coupled, with the ∼110 kyr long allocycles of eccentricitymodulated precession. This resonance may have caused a shift of sensitivity to (obliquitypaced) insolation changes at increasingly higher latitudes, especially during terminations.
Despite the fact that the evidence for energy transfers agrees well with previously published mechanisms, their unprecedentedly detailed description here is an important step forward toward a more comprehensive solution of the problem of the ice ages, because they largely reconcile the mismatch in power spectral density between records that capture the combined effects of deepsea temperatures and landice volumes and those of insolation. Furthermore, we can now state with greater certainty than before that the asymmetry of the Pliocene and Pleistocene ice ages is in very close agreement with NeoMilankovitchism: i.e. the classical Milankovitch theory of astronomical climate forcing and its NeoMilankovitchistic derivatives/superlatives that give greater weight to Earth's internal nonlinear responses. If the bispectral energy transfers documented here indeed track the entire redistribution of power over the spectrum, then the 40 kyr problem of the Pliocene and Early Pleistocene and the ∼110 kyr problem of the Middle and Late Pleistocene are reduced to finding appropriate coupling coefficients that describe the relevant climatic and cryospheric mechanisms – a promising outlook.
The raw data that went into the LR04 stack, which was recompiled for the purposes of this study, can be downloaded at http://lorrainelisiecki.com/LR04cores.zip (last access: 3 January 2019). The original LR04 stack can be downloaded at http://lorrainelisiecki.com/LR04stack.txt and at https://doi.org/10.1594/PANGAEA.704257 (Lisiecki and Raymo, 2005).
The supplement related to this article is available online at: https://doi.org/10.5194/cp1519592019supplement.
DL and ATMdB designed the study, made the figures and wrote the manuscript. ATMdB adjusted the MATLAB scripts for the purposes of this study.
The authors declare that they have no conflict of interest.
We thank Xavier Bertin for inviting Diederik Liebrand to present this work at an early stage at the Université de La Rochelle. We are grateful to Tim E. van Peer for making us aware of the difficulties in quantifying robust geometries of nonergodic records. We thank Lorraine E. Lisiecki and Maureen E. Raymo, and Seonmin Ahn and colleagues for making their benthic foraminiferal δ^{18}O stacks available online. We thank Thomas H. C. Herbers and Gerben Ruessink for providing bispectral analysis MATLAB scripts, which Anouk T. M. de Bakker developed further for the purposes of this study. We thank Qiuzhen Yin for the editorial handling of the manuscript. We thank Mathieu Martinez and the anonymous referee for their valuable comments during the review process.
The article processing charges for this openaccess publication were covered by the University of Bremen.
This paper was edited by Qiuzhen Yin and reviewed by Mathieu Martinez and one anonymous referee.
AbeOuchi, A., Saito, F., Kawamura, K., Raymo, M. E., Okuno, J., Takahashi, K., and Blatter, H.: Insolationdriven 100 000year glacial cycles and hysteresis of icesheet volume, Nature, 500, 190–194, https://doi.org/10.1038/nature12374, 2013.
Agassiz, L.: Études sur les glaciers, Neuchâtel, Jent et Gassmann, available at: https://archive.org/details/etudessurlesgla00agasgoog/page/n10 (last access: 30 March 2019), 1840.
Ahn, S., Khider, D., Lisiecki, L. E., and Lawrence, C. E.: A probabilistic Pliocene–Pleistocene stack of benthic δ^{18}O using a profile hidden Markov model, Dynamics and Statistics of the Climate System, 2, 1–16, https://doi.org/10.1093/climsys/dzx002, 2017.
Bak, P.: How nature works: the science of selforganized criticality, Copernicus, New York, 1996.
Bak, P., Tang, C., and Wiesenfeld, K.: Selforganized criticality – an explanation of 1/f noise, Phys. Rev. Lett., 59, 381–384, https://doi.org/10.1103/PhysRevLett.59.381, 1987.
Bak, P., Tang, C., and Wiesenfeld, K.: Selforganized criticality, Phys. Rev. A, 38, 364–374, https://doi.org/10.1103/PhysRevA.38.364, 1988.
Beerling, D. J. and Royer, D. L.: Convergent Cenozoic CO_{2} history, Nat. Geosci., 4, 418–420, https://doi.org/10.1038/ngeo1186, 2011.
Berger, A., Loutre, M. F., and Mélice, J. L.: Equatorial insolation: from precession harmonics to eccentricity frequencies, Clim. Past, 2, 131–136, https://doi.org/10.5194/cp21312006, 2006.
Berger, A. L.: Support for astronomical theory of climatic change, Nature, 269, 44–45, https://doi.org/10.1038/269044a0, 1977.
Bintanja, R. and van de Wal, R. S. W.: North American icesheet dynamics and the onset of 100 000year glacial cycles, Nature, 454, 869–872, https://doi.org/10.1038/nature07158, 2008.
Bosmans, J. H. C., Drijfhout, S. S., Tuenter, E., Hilgen, F. J., and Lourens, L. J.: Response of the North African summer monsoon to precession and obliquity forcings in the ECEarth GCM, Clim. Dynam., 44, 279–297, https://doi.org/10.1007/s003820142260z, 2015a.
Bosmans, J. H. C., Hilgen, F. J., Tuenter, E., and Lourens, L. J.: Obliquity forcing of lowlatitude climate, Clim. Past, 11, 1335–1346, https://doi.org/10.5194/cp1113352015, 2015b.
Broecker, W. S.: Absolute dating and astronomical theory of glaciation, Science, 151, 299–304, https://doi.org/10.1126/science.151.3708.299, 1966.
Calov, R. and Ganopolski, A.: Multistability and hysteresis in the climatecryosphere system under orbital forcing, Geophys. Res. Lett., 32, L21717, https://doi.org/10.1029/2005GL024518, 2005.
Chalk, T. B., Hain, M. P., Foster, G. L., Rohling, E. J., Sexton, P. F., Badger, M. P. S., Cherry, S. G., Hasenfratz, A. P., Haug, G. H., Jaccard, S. L., MartinezGarcia, A., Palike, H., Pancost, R. D., and Wilson, P. A.: Causes of ice age intensification across the MidPleistocene Transition, P. Natl. Acad. Sci. USA, 114, 13114–13119, https://doi.org/10.1073/pnas.1702143114, 2017.
Chaudhuri, P. and Marron, J. S.: SiZer for exploration of structures in curves, J. Am. Stat. Assoc., 94, 807–823, https://doi.org/10.2307/2669996, 1999.
Cheng, H., Edwards, R. L., Sinha, A., Spotl, C., Yi, L., Chen, S. T., Kelly, M., Kathayat, G., Wang, X. F., Li, X. L., Kong, X. G., Wang, Y. J., Ning, Y. F., and Zhang, H. W.: The Asian monsoon over the past 640 000 years and ice age terminations, Nature, 534, 640–646, https://doi.org/10.1038/nature18591, 2016.
Clark, P. U. and Pollard, D.: Origin of the middle Pleistocene transition by ice sheet erosion of regolith, Paleoceanography, 13, 1–9, https://doi.org/10.1029/97pa02660, 1998.
Clark, P. U., Archer, D., Pollard, D., Blum, J. D., Rial, J. A., Brovkin, V., Mix, A. C., Pisias, N. G., and Roy, M.: The middle Pleistocene transition: characteristics. mechanisms, and implications for longterm changes in atmospheric pCO_{2}, Quaternary Sci. Rev., 25, 3150–3184, https://doi.org/10.1016/j.quascirev.2006.07.008, 2006.
Collis, W. B., White, P. R., and Hammond, J. K.: Higherorder spectra: The bispectrum and trispectrum, Mech. Syst. Signal Pr., 12, 375–394, https://doi.org/10.1006/mssp.1997.0145, 1998.
Crowhurst, S. J., Palike, H., and Rickaby, R. E. M.: Carbonate ions, orbits and Mg/Ca at ODP 1123, Geochim. Cosmochim. Ac., 236, 384–398, https://doi.org/10.1016/j.gca.2018.03.013, 2018.
Crucifix, M.: Oscillators and relaxation phenomena in Pleistocene climate theory, Philos. T. R. Soc. A, 370, 1140–1165, https://doi.org/10.1098/rsta.2011.0315, 2012.
Da Silva, A. C., Dekkers, M. J., De Vleeschouwer, D., Hladil, J., Chadimova, L., Slavik, L., and Hilgen, F. J.: Millennialscale climate changes manifest Milankovitch combination tones and Hallstatt solar cycles in the Devonian greenhouse world, Geology, 47, 19–22, https://doi.org/10.1130/G45511.1, 2018.
de Bakker, A. T. M., Tissier, M. F. S., and Ruessink, B. G.: Shoreline dissipation of infragravity waves, Cont. Shelf Res., 72, 73–82, https://doi.org/10.1016/j.csr.2013.11.013, 2014.
de Bakker, A. T. M., Herbers, T. H. C., Smit, P. B., Tissier, M. F. S., and Ruessink, B. G.: Nonlinear infragravitywave interactions on a gently sloping laboratory beach, J. Phys. Oceanogr., 45, 589–605, https://doi.org/10.1175/JPOD140186.1, 2015.
de Bakker, A. T. M., Tissier, M. F. S., and Ruessink, B. G.: Beach steepness effects on nonlinear infragravitywave interactions: A numerical study, J. Geophys. Res.Oceans, 121, 554–570, https://doi.org/10.1002/2015jc011268, 2016.
de Boer, B., van de Wal, R. S. W., Lourens, L. J., and Bintanja, R.: Transient nature of the Earth's climate and the implications for the interpretation of benthic δ^{18}O records, Palaeogeogr. Palaeocl., 335, 4–11, https://doi.org/10.1016/j.palaeo.2011.02.001, 2012.
deMenocal, P. B.: PlioPleistocene African climate, Science, 270, 53–59, https://doi.org/10.1126/science.270.5233.53, 1995.
Doering, J. C. and Bowen, A. J.: Parametrization of orbital velocity asymmetries of shoaling and breaking waves using bispectral analysis, Coast. Eng., 26, 15–33, https://doi.org/10.1016/03783839(95)00007X, 1995.
Elderfield, H., Ferretti, P., Greaves, M., Crowhurst, S., McCave, I. N., Hodell, D. A., and Piotrowski, A. M.: Evolution of ocean temperature and ice volume through the midPleistocene climate transition, Science, 337, 704–709, https://doi.org/10.1126/science.1221294, 2012.
Elgar, S. and Guza, R. T.: Observations of bispectra of shoaling surface gravitywaves, J. Fluid Mech., 161, 425–448, https://doi.org/10.1017/S0022112085003007, 1985.
Elgar, S.: Relationships involving third moments and bispectra of a harmonic process, IEEE T. Acoust. Speech, ASSP35, 1725–1726, https://doi.org/10.1109/TASSP.1987.1165090, 1987.
Elgar, S. and Sebert, G.: Statistics of bicoherence and biphase, J. Geophys. Res.Oceans, 94, 10993–10998, https://doi.org/10.1029/JC094iC08p10993, 1989.
Farmer, J. R., Honisch, B., Haynes, L. L., Kroon, D., Jung, S., Ford, H. L., Raymo, M. E., JaumeSegui, M., Bell, D. B., Goldstein, S. L., Pena, L. D., Yehudai, M., and Kim, J.: Deep Atlantic Ocean carbon storage and the rise of 100 000year glacial cycles, Nat. Geosci., 12, 355–360, https://doi.org/10.1038/s4156101903346, 2019.
Ganopolski, A. and Calov, R.: The role of orbital forcing, carbon dioxide and regolith in 100 kyr glacial cycles, Clim. Past, 7, 1415–1425, https://doi.org/10.5194/cp714152011, 2011.
Ganopolski, A., Winkelmann, R., and Schellnhuber, H. J.: Critical insolationCO_{2} relation for diagnosing past and future glacial inception, Nature, 529, 200–203, https://doi.org/10.1038/nature16494, 2016.
Hagelberg, T., Pisias, N., and Elgar, S.: Linear and nonlinear couplings between orbital forcing and the marine δ^{18}O record during the late Neogene, Paleoceanography, 6, 729–746, https://doi.org/10.1029/91PA02281, 1991.
Hagelberg, T. K., Bond, G., and deMenocal, P.: Milankovitch band forcing of submilankovitch climate variability during the Pleistocene, Paleoceanography, 9, 545–558, https://doi.org/10.1029/94pa00443, 1994.
Hasselmann, K.: Stochastic climate models. Part I. Theory, Tellus, 28, 473–485, https://doi.org/10.3402/tellusa.v28i6.11316, 1976.
Hasselmann, K., Munk, W., and MacDonald, G.: Bispectra of ocean waves, in: Proceedings of the Symposium on Time Series Analysis, edited by: Rosenblatt, M., John Wiley, 125–139, 1963.
Haug, G. H. and Tiedemann, R.: Effect of the formation of the Isthmus of Panama on Atlantic Ocean thermohaline circulation, Nature, 393, 673–676, https://doi.org/10.1038/31447, 1998.
Hays, J. D., Imbrie, J., and Shackleton, N. J.: Variations in the Earth's orbit – Pacemaker of the ice ages, Science, 194, 1121–1132, https://doi.org/10.1126/science.194.4270.1121, 1976.
Herbers, T. H. C. and Burton, M. C.: Nonlinear shoaling of directionally spread waves on a beach, J. Geophys. Res.Oceans, 102, 21101–21114, https://doi.org/10.1029/97jc01581, 1997.
Herbers, T. H. C., Russnogle, N. R., and Elgar, S.: Spectral energy balance of breaking waves within the surf zone, J. Phys. Oceanogr., 30, 2723–2737, https://doi.org/10.1175/15200485(2000)030<2723:SEBOBW>2.0.CO;2, 2000.
Hilgen, F. J., Lourens, L. J., Van Dam, J. A., Beu, A. G., Boyes, A. F., Cooper, R. A., Krijgsman, W., Ogg, J. G., Piller, W. E., and Wilson, D. S.: The Neogene Period, in: Geologic Time Scale 2012, vols. 1 and 2, edited by: Gradstein, F. M., Ogg, J. G., Schmitz, M. D., and Ogg, G. M., Elsevier, 923–978, 2012.
Huybers, P.: Glacial variability over the last two million years: an extended depthderived agemodel, continuous obliquity pacing, and the Pleistocene progression, Quaternary Sci. Rev., 26, 37–55, https://doi.org/10.1016/j.quascirev.2006.07.013, 2007.
Huybers, P.: Combined obliquity and precession pacing of late Pleistocene deglaciations, Nature, 480, 229–232, https://doi.org/10.1038/nature10626, 2011.
Huybers, P. and Curry, W.: Links between annual, Milankovitch and continuum temperature variability, Nature, 441, 329–332, https://doi.org/10.1038/nature04745, 2006.
Huybers, P. and Tziperman, E.: Integrated summer insolation forcing and 40 000year glacial cycles: The perspective from an icesheet/energybalance model, Paleoceanography, 23, PA1208, https://doi.org/10.1029/2007pa001463, 2008.
Huybers, P. and Wunsch, C.: A depthderived Pleistocene age model: Uncertainty estimates, sedimentation variability, and nonlinear climate change, Paleoceanography, 19, PA1028, https://doi.org/10.1029/2002pa000857, 2004.
Huybers, P. and Wunsch, C.: Obliquity pacing of the late Pleistocene glacial terminations, Nature, 434, 491–494, https://doi.org/10.1038/nature03401, 2005.
Kender, S., Ravelo, A. C., Worne, S., Swann, G. E. A., Leng, M. J., Asahi, H., Becker, J., Detlef, H., Aiello, I. W., Andreasen, D., and Hall, I. R.: Closure of the Bering Strait caused MidPleistocene Transition cooling, Nat. Commun., 9, 5386, https://doi.org/10.1038/s41467018078280, 2018.
Kennedy, A. B., Chen, Q., Kirby, J. T., and Dalrymple, R. A.: Boussinesq modeling of wave transformation, breaking, and runup. I: 1D, J. Waterw. Port C., 126, 39–47, https://doi.org/10.1061/(ASCE)0733950X(2000)126:1(39), 2000.
King, T.: Quantifying nonlinearity and geometry in time series of climate, Quaternary Sci. Rev., 15, 247–266, https://doi.org/10.1016/02773791(95)000607, 1996.
King Hagelberg, T. and Cole, J.: Determining the role of linear and nonlinear interactions in records of climate change: possible solar influences in annual to centuryscale record, in: Proceedings of the eleventh annual Pacific climate (PACLIM) workshop, edited by: Isaacs, C. M., and Tharp, V. L., Interagency ecological program, California department of water resources, 47–54, 1995.
Konijnendijk, T. Y. M., Ziegler, M., and Lourens, L. J.: On the timing and forcing mechanisms of late Pleistocene glacial terminations: Insights from a new highresolution benthic stable oxygen isotope record of the eastern Mediterranean, Quaternary Sci. Rev., 129, 308–320, https://doi.org/10.1016/j.quascirev.2015.10.005, 2015.
Köppen, W. and Wegener, A.: Die klimate der geologischen vorzeit., Gebrüder Borntraeger, Berlin, 1924.
Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A. C. M., and Levrard, B.: A longterm numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428, 261–285, https://doi.org/10.1051/00046361:20041335, 2004.
Laskar, J., Fienga, A., Gastineau, M., and Manche, H.: La2010: A new orbital solution for the long term motion of the Earth, Astron. Astrophys., 532, A89, https://doi.org/10.1051/00046361/201116836, 2011a.
Laskar, J., Gastineau, M., Delisle, J.B., Farrés, A., and Fienga, A.: Strong chaos induced by close encounters with Ceres and Vesta, Astron. Astrophys., 532, 1–4, https://doi.org/10.1051/00046361/201117504, 2011b.
Liebrand, D., de Bakker, A. T. M., Beddow, H. M., Wilson, P. A., Bohaty, S. M., Ruessink, G., Pälike, H., Batenburg, S. J., Hilgen, F. J., Hodell, D. A., Huck, C. E., Kroon, D., Raffi, I., Saes, M. J. M., van Dijk, A. E., and Lourens, L. J.: Evolution of the early Antarctic ice ages, P. Natl. Acad. Sci. USA, 114, 3867–3872, https://doi.org/10.1073/pnas.1615440114, 2017.
Lisiecki, L. E.: Links between eccentricity forcing and the 100 000year glacial cycle, Nat. Geosci., 3, 349–352, https://doi.org/10.1038/Ngeo828, 2010.
Lisiecki, L. E. and Raymo, M. E.: A PliocenePleistocene stack of 57 globally distributed benthic δ^{18}O records, Paleoceanography, 20, PA1003, https://doi.org/10.1029/2004PA001071, 2005a.
Lisiecki, L. E. and Raymo, M. E.: PliocenePleistocene stack of globally distributed benthic stable oxygen isotope records, PANGAEA, https://doi.org/10.1594/PANGAEA.704257, 2005b.
Lisiecki, L. E. and Raymo, M. E.: PlioPleistocene climate evolution: trends and transitions in glacial cycle dynamics, Quaternary Sci. Rev., 26, 56–69, https://doi.org/10.1016/j.quascirev.2006.09.005, 2007.
Lisiecki, L. E., Raymo, M. E., and Curry, W. B.: Atlantic overturning responses to Late Pleistocene climate forcings, Nature, 456, 85–88, https://doi.org/10.1038/nature07425, 2008.
Liu, Z., OttoBliesner, B. L., He, F., Brady, E. C., Tomas, R., Clark, P. U., Carlson, A. E., LynchStieglitz, J., Curry, W., Brook, E., Erickson, D., Jacob, R., Kutzbach, J., and Cheng, J.: Transient Simulation of Last Deglaciation with a New Mechanism for BollingAllerod Warming, Science, 325, 310–314, https://doi.org/10.1126/science.1171041, 2009.
Lourens, L. J. and Hilgen, F. J.: Longperiodic variations in the Earth's obliquity and their relation to thirdorder eustatic cycles and late Neogene glaciations, Quatern. Int., 40, 43–52, https://doi.org/10.1016/S10406182(96)000602, 1997.
Lourens, L. J., Becker, J., Bintanja, R., Hilgen, F. J., Tuenter, E., van de Wal, R. S. W., and Ziegler, M.: Linear and nonlinear response of late Neogene glacial cycles to obliquity forcing and implications for the Milankovitch theory, Quaternary Sci. Rev., 29, 352–365, https://doi.org/10.1016/j.quascirev.2009.10.018, 2010.
MartinezBoti, M. A., Foster, G. L., Chalk, T. B., Rohling, E. J., Sexton, P. F., Lunt, D. J., Pancost, R. D., Badger, M. P. S., and Schmidt, D. N.: PlioPleistocene climate sensitivity evaluated using highresolution CO2 records, Nature, 518, 49–54, https://doi.org/10.1038/nature14145, 2015.
Meyers, S. R. and Hinnov, L. A.: Northern Hemisphere glaciation and the evolution of PlioPleistocene climate noise, Paleoceanography, 25, PA3207, https://doi.org/10.1029/2009pa001834, 2010.
Milankovitch, M.: Canon of Insolation and the IceAge Problem, 1st ed., Zavod za udžbenike i nastavna sredstva & Muzej nauke i tehnike Srpske akademije nauka i umetnosti, Beograd, 1941.
Muller, R. A. and MacDonald, G. J.: Spectrum of 100 kyr glacial cycle: Orbital inclination, not eccentricity, P. Natl. Acad. Sci. USA, 94, 8329–8334, https://doi.org/10.1073/pnas.94.16.8329, 1997a.
Muller, R. A. and MacDonald, G. J.: Simultaneous presence of orbital inclination and eccentricity in proxy climate records from Ocean Drilling Program Site 806, Geology, 25, 3–6, https://doi.org/10.1130/00917613(1997)025<0003:SPOOIA>2.3.CO;2, 1997b.
Norheim, C. A., Herbers, T. H. C., and Elgar, S.: Nonlinear evolution of surface wave spectra on a beach, J. Phys. Oceanogr., 28, 1534–1551, https://doi.org/10.1175/15200485(1998)028<1534:Neosws>2.0.Co;2, 1998.
Paillard, D., Labeyrie, L., and Yiou, P.: AnalySeries, Macintosh program performs timeseries analysis, EOS Transactions AGU, 77, 379, https://doi.org/10.1029/96EO00259, 1996.
Pelletier, J. D.: The power spectral density of atmospheric temperature from time scales of 10^{−2} to 10^{6} yr, Earth Planet. Sc. Lett., 158, 157–164, https://doi.org/10.1016/S0012821x(98)00051X, 1998.
Pestiaux, P., Van der Mersch, I., Berger, A., and Duplessy, J. C.: Paleoclimatic variability at frequencies ranging from 1 cycle per 10 000 years to 1 cycle per 1000 years: evidence for nonlinear behaviour of the climate system, Climatic Change, 12, 9–37, https://doi.org/10.1007/Bf00140262, 1988.
Pisias, N. G., Mix, A. C., and Zahn, R.: Nonlinear response in the global climate system: evidence from benthic oxygen isotopic record in Core Rc13110, Paleoceanography, 5, 147–160, https://doi.org/10.1029/PA005i002p00147, 1990.
Raymo, M. E. and Nisancioglu, K.: The 41 kyr world: Milankovitch's other unsolved mystery, Paleoceanography, 18, 1011, https://doi.org/10.1029/2002pa000791, 2003.
Raymo, M. E., Lisiecki, L. E., and Nisancioglu, K. H.: Pliopleistocene ice volume, Antarctic climate, and the global δ^{18}O record, Science, 313, 492–495, https://doi.org/10.1126/science.1123296, 2006.
Rial, J. A. and Anaclerio, C. A.: Understanding nonlinear responses of the climate system to orbital forcing, Quaternary Sci. Rev., 19, 1709–1722, https://doi.org/10.1016/S02773791(00)000871, 2000.
Ridgwell, A. J., Watson, A. J., and Raymo, M. E.: Is the spectral signature of the 100 kyr glacial cycle consistent with a Milankovitch origin?, Paleoceanography, 14, 437–440, https://doi.org/10.1029/1999pa900018, 1999.
Roe, G.: In defense of Milankovitch, Geophys. Res. Lett., 33, L24703, https://doi.org/10.1029/2006gl027817, 2006.
Roe, G. H. and Allen, M. R.: A comparison of competing explanations for the 100 000yr ice age cycle, Geophys. Res. Lett., 26, 2259–2262, https://doi.org/10.1029/1999gl900509, 1999.
Rohling, E. J., Foster, G. L., Grant, K. M., Marino, G., Roberts, A. P., Tamisiea, M. E., and Williams, F.: Sealevel and deepseatemperature variability over the past 5.3 million years, Nature, 508, 477–482, https://doi.org/10.1038/nature13230, 2014.
Rutherford, S. and D'Hondt, S.: Early onset and tropical forcing of 100 000year Pleistocene glacial cycles, Nature, 408, 72–75, https://doi.org/10.1038/35040533, 2000.
Shackleton, N. J.: The 100 000year iceage cycle identified and found to lag temperature, carbon dioxide, and orbital eccentricity, Science, 289, 1897–1902, https://doi.org/10.1126/science.289.5486.1897, 2000.
Sun, Y., Chen, J., Clemens, S. C., Liu, Q. S., Ji, J. F., and Tada, R.: East Asian monsoon variability over the last seven glacial cycles recorded by a loess sequence from the northwestern Chinese Loess Plateau, Geochem. Geophy. Geosy., 7, Q12Q02, https://doi.org/10.1029/2006gc001287, 2006.
Tabor, C. R. and Poulsen, C. J.: Simulating the midPleistocene transition through regolith removal, Earth Planet. Sc. Lett., 434, 231–240, https://doi.org/10.1016/j.epsl.2015.11.034, 2016.
Tziperman, E., Raymo, M. E., Huybers, P., and Wunsch, C.: Consequences of pacing the Pleistocene 100 kyr ice ages by nonlinear phase locking to Milankovitch forcing, Paleoceanography, 21, PA4206, https://doi.org/10.1029/2005pa001241, 2006.
Viaggi, P.: δ^{18}O and SST signal decomposition and dynamic of the PliocenePleistocene climate system: new insights on orbital nonlinear behavior vs. longterm trend, Progress in Earth and Planetary Science, 5, 1–37, https://doi.org/10.1186/s406450180236z, 2018.
von Dobeneck, T. and Schmieder, F.: Using rock magnetic proxy records for orbital tuning and extended time series analyses into the super and subMilankovitch band, in: Use of proxies in paleoceanography: examples from the South Atlantic, edited by: Fischer, G. and Wefer, G., SpringerVerlag, Berlin, Heidelberg, 601–633, https://doi.org/10.1007/9783642586460_25, 1999.
Wara, M. W., Ravelo, A. C., and Revenaugh, J. S.: The pacemaker always rings twice, Paleoceanography, 15, 616–624, https://doi.org/10.1029/2000pa000500, 2000.
Werner, M., Heimann, M., and Hoffmann, G.: Isotopic composition and origin of polar precipitation in present and glacial climate simulations, Tellus B, 53, 53–71, https://doi.org/10.3402/tellusb.v53i1.16539, 2001.
Willeit, M., Ganopolski, A., Calov, R., Robinson, A., and Maslin, M.: The role of CO_{2} decline for the onset of Northern Hemisphere glaciation, Quaternary Sci. Rev., 119, 22–34, https://doi.org/10.1016/j.quascirev.2015.04.015, 2015.
Willeit, M., Ganopolski, A., Calov, R., and Brovkin, V.: MidPleistocene transition in glacial cycles explained by declining CO_{2} and regolith removal, Sci. Adv., 5, eaav7337, https://doi.org/10.1126/sciadv.aav7337, 2019.
Wu, P. L., Jackson, L., Pardaens, A., and Schaller, N.: Extended warming of the northern high latitudes due to an overshoot of the Atlantic meridional overturning circulation, Geophys. Res. Lett., 38, L24704, https://doi.org/10.1029/2011gl049998, 2011.
Yiou, P., Ghil, M., Jouzel, J., Paillard, D., and Vautard, R.: Nonlinear variability of the climatic system from singular and power spectra of late Quaternary records, Clim. Dynam., 9, 371–389, https://doi.org/10.1007/s003820050030, 1994.
Ziegler, M., Lourens, L. J., Tuenter, E., Hilgen, F., Reichart, G. J., and Weber, N.: Precession phasing offset between Indian summer monsoon and Arabian Sea productivity linked to changes in Atlantic overturning circulation, Paleoceanography, 25, PA3213, https://doi.org/10.1029/2009pa001884, 2010.