The Middle Pleistocene Transition by frequency locking and slow ramping of internal period

TL;DR

This paper introduces the Ramping with Frequency Locking (RFL) mechanism, explaining the Middle Pleistocene Transition by slow internal period changes and frequency locking, accounting for abrupt glacial cycle length increases without requiring rapid external parameter shifts.

Contribution

The paper proposes the RFL mechanism as a novel explanation for the MPT, linking slow internal period ramping with frequency locking to produce abrupt cycle length changes.

Findings

RFL explains the MPT without rapid external changes.

Frequency locking can produce abrupt cycle shifts with slow parameter variation.

The scheme for detecting RFL in complex models is proposed.

Abstract

The increase in glacial cycle length from approximately $41$ to on average $100$ thousand years around $1$ million years ago, called the Middle Pleistocene Transition (MPT), lacks a conclusive explanation. We describe a dynamical mechanism which we call Ramping with Frequency Locking (RFL), that explains the transition by an interaction between the internal period of a self-sustained oscillator and forcing that contains periodic components. This mechanism naturally explains the abrupt increase in cycle length from approximately $40$ to $80$ thousand years observed in proxy data, unlike some previously proposed mechanisms for the MPT. A rapid increase in durations can be produced by a rapid change in an external parameter, but this assumes rather than explains the abruptness. In contrast, models relying on frequency locking can produce a rapid change in durations assuming only a slow change in an external parameter. We propose a scheme for detecting RFL in complex, computationally expensive models, and motivate the search for climate variables that can gradually increase the internal period of the glacial cycles.

Equations

$$\dot{x}$$

$$\displaystyle\text{ linearly decrease }x(t)\text{ to }0\text{ over {\color[rgb]{1,0,0}{time }}}T_{decay},\text{ repeat}.$$

$${\overline{D}}_{\tau }=\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{i=1}^n{D}_{i,\tau },$$

$$\frac{\Delta {\overline{D}}_{\tau }(\tau )}{\Delta \tau }\approx$$

$$\frac{\Delta {\overline{D}}_{\tau }(T_o,A)}{\Delta T_o}\frac{\Delta T_o(R)}{\Delta R}\frac{\Delta R(\tau )}{\Delta \tau }+$$

$$\frac{\Delta {\overline{D}}_{\tau }(T_o,A)}{\Delta A}\frac{\Delta A(\tau )}{\Delta \tau },$$

$${\overline{D}}_{\tau ,loc}=\frac{1}{|I(\tau )|}\sum _{i\in I(\tau )}{D}_{i,\tau },$$

$${\tau }_K=\frac{n_c-n_d}{\sqrt{n_0^2-n_0{n}_1}},$$

Full text

1 Introduction

Since the beginning of major Northern hemisphere glaciation $2.7$ million years (Myr) ago, Earth has undergone alternating epochs of icy and cold conditions on the one hand, and warm and ice-free conditions on the other (Fig. 1). While these glacial cycles are attested from geological records (Lisiecki and Raymo, 2005; Huybers, 2007; EPICA community members, 2004), there is no single conclusive theory of their origins. Historically, the focus has been to explain the approximately $100$ thousand years (kyr) long glacial cycles that dominate the past $800$ kyr (see (Imbrie and Imbrie, 1979) for a review). But as sediment records later revealed that these cycles were ca $40$ kyr long prior to $1.2$ Myr ago the question arose what caused the shift to approximately $100$ kyr long cycles, called the Middle Pleistocene Transition (MPT) (Clark et al., 2006).

The main strategy to address these questions has been to replicate the palaeoclimatic records using simple models of glacial cycles with few variables, referred to as conceptual models (see (Crucifix, 2012) for a review). One reason for this is that the rather regular and cyclic variations in data suggest that the the main dynamics can be captured by a system of few degrees of freedom, even as the full climate system obviously has a large number of degrees of freedom. None of these models describe the climate system in detail, but they are useful for understanding underlying dynamics. Virtually all models involve insolation variations due to changes to Earth’s orbital configuration relative to the sun, an idea heralded by Adhémar, Croll and Milankovitch (Imbrie and Imbrie, 1979). But the specific role played by insolation variations is still unknown and debated.

Several solutions to the cause of the MPT have been presented within the context of conceptual models. Some mechanisms rely on a bifurcation occurring in the unforced climate system which fundamentally changes how the system operates (Ashwin and Ditlevsen, 2015; Ditlevsen, 2009; Tziperman and Gildor, 2003; Maasch and Salzman, 1990; Huybers and Langmuir, 2017). Other mechanisms invoke a “spontaneous” change, such as a shift between attractors due to subtle changes in insolation (Quinn et al., 2018; Omta et al., 2015) or random fluctuations (Salzman and Verbitsky, 1993; Imbrie et al., 2011), or as a coincidence (Huybers, 2009). A third possible mechanism for the MPT assumes one essential mode of oscillation throughout the Pleistocene and relies crucially on the interaction between insolation variations and an increasing internal period. This mechanism, previously imprecisely referred to as phase/frequency locking and non-linear resonance — but here Ramping with Frequency Locking (RFL) — is the focus of this paper.

The main appeal of this mechanism is that nothing special had to occur in the climate system over the MPT (Huybers, 2007); it is only required that the internal period was ramped slowly — interactions with forcing are enough to cause an abrupt increase in durations between glacial terminations (Fig. 1, bottom panel).

The first publications where RFL was used (Paillard, 1998; Paillard and Parrenin, 2004) did not explain why the durations between glacial transitions increased abruptly over the MPT. Ashkenazy and Tziperman (2004); Ashkenazy (2006) hinted how frequency locking (therein called phase locking) could produce an abrupt increase in duration, by showing diagrams of average duration as a function of a system parameter (Devil’s staircases). Huybers (2007) was first to both show a model trajectory of ice volume using the mechanism, and to attribute the effect to “skipping of obliquity cycles”, a frequency locking effect. Recently, Feng and Bailer-Jones (2015); Mitsui et al. (2015); Daruka and Ditlevsen (2015); Tzedakis et al. (2017) alluded to the mechanism, but neither emphasised that frequency locking can explain the MPT assuming only a slow linear change in a climate parameter. Instead, by ramping some parameter in a way that mimics the rapid change in durations over the MPT, they prescribe an abrupt increase in period over the MPT rather than explaining it. Here, we for the first time properly define RFL and emphasise its generality.

Following (Huybers, 2007), we question the common assumption that climate entered a stationary state in the late Pleistocene, and instead argue that the sequence of durations between glacial terminations is consistent with a slow increase of the internal period of the climate system until present (Fig. 1 b)). According to this view, the typical durations between transitions changed from ${\sim}40$ kyr to ${\sim}80$ kyr around $1200$ kyr ago, after which they increased gradually in the mean to present time, with the last duration being ${\sim}120$ kyr long.

Here, we first aim to explain RFL in a clear way, using a harmonically forced simple model. We use harmonic (pure sine) forcing because it makes frequency locking concepts clearer, while still producing qualitatively similar behaviour to astronomical forcing curves. We should not expect model runs with such simplified forcing to agree well with data, however. We use forcing with period $41$ kyr, corresponding to the main period of obliquity variations (Berger, 1978), which determine the total insolation integrated over the summer at Northern latitudes (Huybers, 2006).

We then define RFL, specify a class of models able to reproduce the MPT using the mechanism, and propose a decomposition of model components to understand the abruptness of the MPT. We consider evidence in data for a $40$ to $80$ kyr shift in durations between terminations and a subsequent gradual increase, and why this supports RFL in favour of some other mechanisms for the MPT. We then discuss how insights from harmonic forcing relate to non-harmonic forcing, how RFL can be detected in complex and computationally expensive models, and some climate variables that can cause an increase in the internal period of the glacial cycles.

2 The idea behind Ramping with Frequency Locking

We illustrate Ramping with Frequency Locking (RFL) using a deterministic and continuous time version of the H07 model (Huybers, 2007) (see Fig. 2). The model is arguably the simplest to represent alternating stages of intrinsic growth and decay of ice sheets, with the growth state ending abruptly as a critical ice volume is reached. It is is an integrate-and-fire threshold model conceptually very similar to the models in (van der Pol and van der Mark, 1927; Imbrie and Imbrie, 1980; Paillard, 1998; Ashkenazy and Tziperman, 2004; Huybers, 2007; de Saedeleer et al., 2013; Imbrie et al., 2011; Parrenin and Paillard, 2003, 2012; Glass and Mackey, 1979).

Physically, sudden and rapid deglaciation has been explained e.g. with isostatic rebound Oerlemans (1980), rapid $CO_{2}$ outgassing (Paillard and Parrenin, 2004) and rapid loss of Northern hemisphere sea ice cover (Gildor and Tziperman, 2000).

We assume that ice volume $x(t)$ grows at a constant rate $\mu$ in a glacial state until it reaches a threshold $\theta(t)$. Then deglaciation starts, whereby ice volume decays to [math] over a fixed time $T_{decay}=10$ kyr:

[TABLE]

Small perturbations to the model, such as having a constant rate of decay instead of a fixed time, does not qualitatively affect its behaviour.

We split $\theta(t)$ into a forcing term $A\cdot F(t)$ – a zero-mean sum of periodic components – and a ramping term $R(t)$: $\theta(t)=R(t)+A\cdot F(t)$.

In the limit of constant ramping $R(t)=R_{0}$ and zero forcing $A=0$, the system has a constant internal period of oscillation $T_{o}=\frac{R_{0}}{\mu}+T_{decay}$ (subscript $o$ for oscillator), see Fig. 2 a). But if the threshold increases slowly over time, for instance linearly $\theta(t)=R(t)=R_{0}+R_{1}t$ as in Fig. 3 a)i), then the internal period $T_{o}(t)=\frac{R_{0}}{\mu}+T_{decay}+\frac{R_{1}}{\mu}t$ also increases slowly (Fig. 3 b)i)). As there is no forcing, the durations between glacial terminations follow $T_{o}(t)$ closely.

However, with periodic forcing

$\theta(t)=R_{0}+R_{1}t+A\sin{(2\pi/T_{f})}$, durations $D_{i}$ are near multiples of the forcing period $D_{i}\approx NT_{f}$, $N\in\mathbb{N}$ (Fig. 3 b)i)). Roughly speaking, the multiple that is realised is the one closest to the internal period $T_{o}$. This phenomenon, called frequency locking (Pikovsky et al., 2001), has been studied extensively over the past century (e.g. van der Pol and van der Mark (1927); Cartwright and Littlewood (1945); Glass and Mackey (1979); Le Treut and Ghil (1983); Tziperman et al. (2006)).

In Fig 3 a)i) the durations $D_{i}$ change abruptly from $1\times T_{f}$ to $2\times T_{f}$ and finally $3\times T_{f}$. These abrupt changes in durations resulting from a gradual change in an underlying parameter is one possible dynamical mechanism behind the MPT.

We call the mechanism Ramping with Frequency Locking (RFL), rather than non-linear resonance, phase locking or frequency locking as it has previously been called. This we do to emphasise both that an internal period must increase gradually over time (ramping), and that the internal oscillations must be locked to external forcing. This is opposed to e.g. the mechanism in (Omta et al., 2015), which realises the MPT through jumps between coexisting frequency locked solutions.

We note that RFL is a special case of “slow passage through bifurcation” (e.g. (Do and Lopez, 2012; Baer et al., 1989)), for which the bifurcations typically are saddle-node bifurcations of limit cycles marking transitions in and out of frequency locking regions (Pikovsky et al., 2001). (However, see e.g. Guckenheimer et al. (2003); Levi (1990) for other relevant bifurcations).

Finally, we note that H07 is an illustrative example of RFL, and not representative of all glacial cycle models. However, the rapid jumps between frequency locking regions occur generically in a broad class of models, defined next.

3 A formal description of RFL

H07 (Fig. 2) is just one particular model capable of realising the MPT through RFL. We could simply call these self-sustained oscillators, but we aim to be more precise and to establish notation.

First, we naturally require the model to be a dynamical system, such that there is an evolution rule $f(t,x)$ taking a state $x(t)$ forward in time $t$. We identify the model with the evolution rule and denote it $f$ (without arguments) for brevity. $f(t,x),x(t)$ and $t$ can be very general, for instance; $t$ can be continuous or discrete, $x$ can be of any dimension, and $f(t,x)$ can e.g. be a piecewise smooth ODE paired with a switching rule, as for H07.

We also require $f$ to be forced by a continuous zero-mean sum of periodic components $A(t)F(t)$ with an amplitude $A(t)$, called the forcing. We further require that $f$ is parametrised by a set of parameters $p(t)$, whose time-varying subset $R(t)$ is called the ramping. Thus we can write $f=f(t,x,R(t),A(t)F(t))$.

We define the frozen system f_{\tau}$$:=f(t,x,R(\tau),A(\tau)F(t)) as $f$ with parameters frozen at time $t=\tau$. Importantly, we require that $f$ is a self-sustained oscillator with internal period $T_{o}(R(\tau))$, meaning that every solution to $f_{\tau}$ with $A(t)\equiv 0$ tends asymptotically (as $t\to\infty$) to a periodic solution with period $T_{o}(R(\tau))$. For RFL to be relevant we require that $T_{o}(R(\tau))$ increases as a function of $\tau$. This is the ramping part of RFL.

The frequency locking part of RFL comes from the response of $f$ to non-zero but constant forcing $A(\tau)$. For small and medium size $A(\tau)$, asymptotic solutions to the frozen system $f(t,x,R(\tau),A(\tau)F(t))$, are generically periodic with periods related rationally to the forcing periods (Pikovsky et al., 2001). The oscillator period can for instance be twice that of the forcing period. If so, the oscillator period (and therefore frequency) remains constant on open sets of parameters and we say that solutions are frequency locked to the forcing (we return to this in Section 4).

The essence of RFL is that the period of the frozen system can change rapidly as function of $T_{o}(R)$ when a ramped parameter causes the system to switch between frequency locking regions.

However, some remarks are in place. Firstly, the system with time varying parameters $f(t,x(t),R(t),A(t)F(t))$ is not the same as the frozen system $f(t,x(t),R(t),A(\tau)F(\tau))$ since solutions to the former cannot equilibrate to solutions of the latter in finite time. We return to differences between the two systems in Section 5 but until then we focus on the frozen system.

Secondly, the period of an oscillator is not the same as the length of individual “cycles”. For instance, around $-1350$ kyr in Fig. 3, short and long “cycles” alternate. This makes the average time between terminations $61.5$ kyr, whereas the period (time until repetition, two large peaks) is $123$ kyr. Therefore, we instead characterise local behaviour with the average duration

[TABLE]

where $D_{i,\tau}$ denotes the $i$:th duration between successive crossings of a fixed threshold for the frozen system $f_{\tau}$.

For some models like H07 glacial terminations are natural such thresholds. For other models Poincaré sections can be considered instead (Pikovsky et al., 2001).

4 Breaking down the dependency of $\overline{D}_{\tau}$ on $\tau$

Comparing Fig. 3 b)i) and b)ii) shows that $\overline{D}_{\tau}$ can rise steeply both from frequency locking effects under a gradual change of parameter (Fig. 3 b)i)) and from ramping of a climate parameter rapidly (Fig. 3 b)ii)).

We wish to break down the contribution to the local change in average duration from these effects and do so by considering the change $\Delta\overline{D}_{\tau}$ under a small perturbation $\Delta\tau$:

[TABLE]

where e.g. $\frac{\Delta\overline{D}_{\tau}(T_{o},A)}{\Delta T_{o}}:=\frac{\overline{D}_{\tau}(T_{o}+\Delta T_{o},A)-\overline{D}_{\tau}(T_{o},A)}{\Delta T_{o}}$, and where we have neglected higher order terms. This approximation is generally better the smaller $\Delta\tau$ is. As $\Delta\tau\to 0$, (3) tends to the chain rule, but since $\frac{\Delta\overline{D}_{\tau}(T_{o},A)}{\Delta T_{o}}=0$ wherever differentiable (see Section 4.1) it is more appropriate to consider $\Delta\overline{D}_{\tau}$ over short intervals of time $\Delta\tau$. In what follows we restrict ourselves to $\frac{\Delta A(\tau)}{\Delta\tau}=\nobreak 0$.

Eq. (3) says that the rate of change (abruptness) in time of the average duration $\overline{D}_{\tau}$ is approximately the product of the rates at which $R(t)$ changes with time, $T_{o}$ changes with $R$, and $\overline{D}_{\tau}$ changes with $T_{o}$. Our point is that each of $\overline{D}_{\tau}$, $T_{o}$ and $R$ can contribute to an abrupt change of $\overline{D}_{\tau}$ at the MPT, but they have different interpretations from a modelling perspective. We discuss these factors next.

4.1 $\overline{D}_{\tau}(T_{o},A)$: Arnold tongues and Devil’s staircases

The average duration $\overline{D}_{\tau}(T_{o},A)$ as a function of internal period $T_{o}$ and forcing amplitude $A$ describes the frequency locking contribution to changes to $\overline{D}_{\tau}$ over time $\tau$.

Frequency locking can be visualised in Arnold tongue diagrams; Fig. 4 a) reveals regions of constant average duration $\overline{D}_{\tau}$ in $(T_{o},A)$ space called Arnold tongues (Pikovsky et al., 2001; Crucifix, 2013; de Saedeleer et al., 2013). Inside major 1:N tongues, solutions are periodic with period $N$ times the forcing period $T_{f}=\nobreak 41$ kyr, as evidenced in (Fig. 4 a)). Minor tongues emanate at $A=0$ from other rationals of $T_{f}$, and in between them are quasiperiodic solutions. (We show only $M:N,\leavevmode\nobreak\ M=\nobreak\{1,2\}$ Arnold tongues, defined numerically as sets for which $|\overline{D}_{\tau}-T_{f}\frac{N}{M}|<\nobreak 0.5$. $\overline{D}_{\tau}$ is estimated over $6$ million years.)

A change in $\tau$ that in turn leads to a change in $T_{o}(R(\tau))$ traces out a path in $(A,T_{o})$ space (black and magenta lines in Fig. 4 a). Such a path represents the change in system state as one or more parameters change in time over the MPT. The paths in Fig. 4 a) pass through the major 1:1, 1:2 and 1:3 locking tongues, in which there are respectively $1$, $2$ and $3$ forcing periods per oscillator period. We learn that for larger $A$, a larger portion of the path stays inside the major 1:N tongues, an observation also made in (Ashkenazy, 2006).

Another way of visualising the change in average duration $\overline{D}_{\tau}$ as a function of $T_{o}$ are Devil’s staircases (Pikovsky et al., 2001) (Fig. 4 b)), in which the forcing amplitude $A$ is fixed. We see that the average duration $\overline{D}_{\tau}$ is constant within Arnold tongues and that the staircase for larger $A$ contains longer steps of constant duration, as predicted from Fig. 4 a). Hence, stronger forcing influence tends to cause more abrupt changes to the average period.

4.2 $T_{o}(R)$ and $R(\tau)$

The function $T_{o}(R)$, if continuous and monotonic, stretches and squeezes Arnold tongues by scaling the independent variable $T_{o}$ of $\overline{D}_{\tau}(T_{o})$. In the Ashkenazy model (Ashkenazy, 2006), for instance, a faster-than-linearly increasing $T_{o}(R)$ makes the 1:2 and 1:3 Arnold tongues, as a function of ice volume threshold, narrower and more closely spaced than the 1:1 tongue.

Ramping $R(\tau)$ continuously and monotonically, similarly stretches and squeezes Arnold tongues. For instance, in Fig. 5 a sigmoidal ramping $R(t)$ makes the 1:2 Arnold tongue narrower compared to a linear change of $R(t)$. Fig. 3 further illustrates this, showing model runs for either a sigmoidally ($R(t)=\nobreak 26+50\times\nobreak(\tanh(\frac{t+1000}{300})+1)$) or a linearly $R(t)=\nobreak 26+\nobreak 0.05\times(t+2000)$ ramped threshold. The sigmoidal ramping accelerates the increase in average duration around $-1000$ kyr, making the transition more abrupt. Note that the parameters in the functions $R(\tau)$ in Fig. 5 and Fig. 3 are different.

4.3 The roles of $\overline{D}_{\tau}(T_{o},A)$, $T_{o}(R)$ and $R(\tau)$ in reproducing the MPT

All of $\overline{D}_{\tau}(T_{o},A)$, $T_{o}(R)$ and $R(\tau)$ govern the average duration $\overline{D}_{\tau}$ and are able to cause an abrupt change of it, like the one observed at the MPT. From a modelling point of view, however, the functions carry different assumptions and are compatible with different hypotheses.

A model having an abrupt change due to $\overline{D}_{\tau}(T_{o},A)$ relies on frequency locking properties, and assumes only slowly varying functions $T_{o}(R)$ and $R(\tau)$. Hence, the internal period is assumed to change slowly with model parameters, and parameters are assumed to change slowly in time. Such a model, relying on few assumptions about the climate system, makes full use of the RFL mechanism. The models in (Paillard, 1998; Paillard and Parrenin, 2004; Huybers, 2007) and H07 are of this kind.

A model which relies predominantly on $T_{o}(R)$ for an abrupt change in average duration $\overline{D}_{\tau}$ is also consistent with a slowly changing external parameter $R(\tau)$, but a particular function $T_{o}(R)$ requires a physical explanation.

A model relying on a rapidly changed external parameter $R(\tau)$ does not need frequency locking properties of $\overline{D}_{\tau}(T_{o},A)$ or a non-linear response of internal dynamics to the parameter $T_{o}(R)$. However, such a model prescribes the abrupt change in average duration at the MPT rather than explaining the dynamics behind it. Therefore, such an explanation requires justification for the rapidly changed external parameter. The models in (Tzedakis et al., 2017; Mitsui et al., 2015; Ashkenazy and Tziperman, 2004; Daruka and Ditlevsen, 2015) can be said to fall under this category, although they also achieve some abruptness through $\overline{D}_{\tau}(T_{o},A)$.

5 Validity of the quasi-static approximation $f\sim f_{\tau}$

The quasistatic approximation is the approximation that parameters $R(\tau)$ change so slowly that the local average duration of $f$ at time $t=\nobreak\tau$

[TABLE]

is equal to $\overline{D}_{\tau}$. $I(\tau)$ is the set of indices of durations $D_{i,\tau}$ within a time interval $[\tau-\tau_{0},\tau+\tau_{1}]$ around $\tau$, with $\tau_{0},\tau_{1}>0$. If $I(\tau)=\emptyset$, then we define $\overline{D}_{\tau,loc}=0$.

If the quasistatic approximation holds, then the average duration, Arnold tongue diagrams and Devil’s staircases calculated for the frozen system $f_{\tau}$ provide accurate information about local dynamics of $f$.

However, if $R(t)$ and/or $T_{o}(R(t))$ change rapidly around $t=\tau$, then there are two sources of discrepancy between $\overline{D}_{\tau,loc}$ and $\overline{D}_{\tau}$.

The first comes from that the length of the interval of time needed for a good average may be long relative to the local change of $\overline{D}_{\tau}$ for the system $f_{\tau}$. Fig. 6 illustrates that for a $9\times 41=\nobreak 369$ kyr-periodic solution, a long interval is needed to get a local average duration $\overline{D}_{\tau,loc}$ in agreement with $\overline{D}_{\tau}$. At the same time, a long averaging interval fails to capture abrupt changes to $\overline{D}_{\tau}$.

The second is that solutions to $f$ may fail to track solutions to $f_{\tau}$. This occurs if the “frozen” attractor of $f_{\tau}$ changes (in some sense) at a fast rate, and if solutions attract to the frozen attractor at a slow rate. Quantifying these rates in a coordinate- and model-independent way seems difficult, however.

A candidate measure of rate of attraction is the maximal Lyapunov exponent of the return map mapping one transition time to another (Pikovsky et al., 2001). This can be normalised to a common time scale between models and is coordinate independent. However, since it is only a local measure it neglects the time it takes to enter a small neighbourhood of the attractor. This time can in practice dominate, as is the case in the standard circle model (not shown, model described in (Pikovsky et al., 2001)).

The local change in average duration $|\frac{\Delta\overline{D}_{\tau}}{\Delta\tau}|$ is a candidate measure of rate of change of an attractor of $f_{\tau}$, since it exists in all models $f$ and is coordinate-independent. It is ambiguous how large $\Delta\tau$ should be, however. Furthermore, the average duration $\overline{D}_{\tau}$ is only a proxy for the position of an attractor in phase space; an attractor can move even if $\frac{\Delta\overline{D}_{\tau}}{\Delta\tau}=\nobreak 0$. This explains the consistent deviation of single durations from the predicted and locally constant $\overline{D}_{\tau}=\nobreak 82$kyr in the inset of Fig. 3 b) i).

6 Is there a 100 kyr world?

The late Pleistocene (${\sim}800$ – [math] kyr) is sometimes referred to as the “$100$ kyr world”, carrying the implicit notion that the Earth system has settled in a stationary mode with a dominant time scale of $100$ kyr (Fig. 1 a)). This view, originating from the closeness to the $100$ kyr component of eccentricity (an astronomical parameter), is supported by the rate of increase of mean ice volume seemingly levelling off (Clark et al., 2006; Mudelsee and Schulz, 1997), and that the Fourier spectrum over the last ${\sim}800$ kyr is centred around $100$ kyr.

We propose on the contrary, following (Huybers, 2007), that the glacial period increased gradually from ${\sim}80$ kyr around $-1200$ kyr to ${\sim}120$ kyr at present day. The change from ${\sim}40$ to ${\sim}80$ kyr long cycles at $-1200$ kyr can be a shift from $1\times 41$ to $2\times 41$ kyr obliquity frequency locking, and/or $2\times 21$ to $4\times 21$ kyr precession locking. We base this claim on durations between major glacial terminations and a wavelet spectrum of the LR04 stack (see Fig. 1 b)); both quantities increase rather rapidly around $-1200$ kyr and show a steady but irregular increase towards present time.

6.1 Identifying the shift to longer periods

While Huybers (2007) observed that the mean period of global ice volume variations increases over time, we make the stronger claim that an abrupt shift from $40$ to $80$ kyr long durations occurred around $-1200$ kyr. We base this claim on our identification of major glacial terminations, which unlike spectral decomposition ignores glacial cycle shape and is unaffected by time-frequency resolution.

A disadvantage of using glacial termination events is that it is unclear what constitutes a major termination, and whether it is meaningful to characterise glacial cycles by termination events. Nevertheless, we believe that our identification of major terminations is sufficiently robust to support the claim that the duration shifted abruptly from ${\sim}40$ to ${\sim}80$ kyr around $-1200$ kyr.

6.2 Testing for trend after the MPT

It appears that the durations between successive glacial terminations are increasing over time starting at the onset of the MPT around $-1200$ kyr.

We evaluate whether this trend is statistically significant, using a variation on the Mann-Kendall test (Mann, 1945; Kendall, 1955). Our null hypothesis H0 is that the sequence of thirteen durations from $-1126$ kyr until present is generated by a process with stationary mean, and that any observed monotonicity is by chance. Since $\tilde{D}_{i}=D_{i}-D_{mean}$, successive deviations from the mean duration $D_{mean}=91$ kyr are correlated, we immediately reject a white noise process as assumed in the standard Mann-Kendall test. Instead, we model them as an AR(1) process, such that $\tilde{D}_{i+1}=\alpha\tilde{D}_{i}+\sigma_{d}\xi_{i}$, where $\xi_{i}$ are independent Gaussian zero mean and unit variance elements. The parameters $\alpha=0.6$ and $\sigma_{d}=14.5$kyr are the standard estimates of lag 1 and 0 autocorrelation coefficients respectively.

We test the hypothesis using the Kendall $\tau$ test statistic for monotonicity $\tau_{K}$, based on the number of ordered and disordered pairs in a sequence. $\tau_{K}=1$ for a perfectly ordered sequence and $\tau_{K}=0$ for sequence with equally many ordered and disordered pairs (see Appendix A for a definition of $\tau_{K}$). We evaluate $\tau_{K}$ for $2\cdot 10^{4}$ samples of the AR(1) process. As indicated in Fig. 7, it is unlikely $(p<0.05)$ to observe the test statistic in durations from data, assuming that the durations follow an AR(1) process. Therefore, we reject the null hypothesis of no trend.

Adding age model uncertainty to the Monte Carlo sequences of durations only makes it more difficult to reject H0. Furthermore, slightly different choices of major glacial terminations, or the use of an untuned record, does not influence the conclusion of the test.

6.3 Consequences for modelling the MPT

Some explanations for the MPT do not reproduce the sequence of successively longer durations between glacial terminations in data as naturally as RFL. Instead, they produce long period cycles at the onset of the MPT which shorten towards the present as a parameter is ramped.

The Maasch and Salzman 1990 model in Fig. 8 is one such model (Maasch and Salzman, 1990). The inconsistency with data is evident when comparing the model durations with those in the LR04 stack (Fig. 1). Another such model is the Tziperman and Gildor 2003 model (Tziperman and Gildor, 2003).

Although different dynamical mechanisms are at play in these models, they have in common that a long period limit unforced cycle emerges near a region of slow motion in phase space. As a parameter is varied, the limit cycle moves farther from this region, shortening the internal period.

RFL on the other hand naturally explains both a sudden shift from $40$ kyr to $80$ kyr cycles and a gradual increase towards longer cycles, since the system can respond both smoothly and abruptly to an increasing internal period, due to the Devil’s staircase structure (e.g. Fig. 4). We interpret the progression of durations as evidence against models like Maasch and Salzman 1990 and Tziperman and Gildor 2003, and for mechanisms that naturally produce increasing glacial cycle length, such as RFL.

7 Non-harmonic forcing

RFL is not restricted to harmonic forcing, but occurs also for astronomical, non-harmonic forcing. This is for instance the case for the Paillard and Parennin 2004 model (Fig. 10), forced by summer solstice insolation at $65$ degrees North (65Nss, Fig. 9 a)). As a parameter is increased linearly, durations first cluster around $41$ kyr, then shift abruptly to cluster around $82$ kyr at $-1000$ kyr, after which they increase gradually until present. The shift to $80$ kyr durations is later than in proxy data (Fig. 1) and there are some short and long durations not clear in the proxy record, but overall the glacial terminations coincide well.

Multi-frequency forcing generally produces Devil’s staircases with shorter steps of constant duration, making them look “smooth” (e.g. Fig. 11). This is apparently a problem for RFL since it relies on rapid jumps in durations. However, RFL can still be relevant as demonstrated by H07 forced by an equal amount of obliquity and precession (Fig. 9 b), Fig. 12 and Fig. 11). Such forcing corresponds well to e.g. caloric summer insolation or integrated insolation above a threshold (Huybers, 2011; Tzedakis et al., 2017). The median and mode of the distribution of durations change more abruptly than the mean, which reflects that the gradual increase in average duration is caused by a gradual redistribution of durations between clusters, rather than a gradual increase of the most typical durations. In a simulation with time-dependent ramping parameter, the local-in-time distribution of durations cannot be sampled well. Therefore the majority of the realised durations come from the dominant clusters of durations, which can give the impression that durations shift rapidly, in spite of the average duration changing gradually (Fig. 12 and Fig. 11)).

Multi-frequency forcing gives rise to many interesting phenomena regarding predictability of solutions, see for instance (Tziperman et al., 2006; Crucifix, 2013; Grebogi et al., 1984; Mitsui et al., 2015; de Saedeleer et al., 2013; Le Treut and Ghil, 1983; Ashwin et al., 2018; Imbrie and Imbrie, 1980). Importantly, however, these phenomena are not essential to RFL. Whether solutions are truly frequency locked or depend on initial conditions is irrelevant, as long as durations undergo abrupt change and tend to cluster.

We conclude that RFL, clearly understood under periodic forcing, also is relevant for astronomical forcing. Indeed, recent studies provide evidence for the long-standing hypothesis that a combination of precession and obliquity paces the glacial cycles (Feng and Bailer-Jones, 2015; Huybers, 2011; Tzedakis et al., 2017). Differences between periodic and multi-frequency forcing exist, but are not crucial for modelling the MPT with RFL.

8 Relevance for complex models and physical mechanisms

We see two practical uses of our description of RFL: To guide modelling of the MPT in complex models, and to drive the search for slowly changing climate variables.

8.1 Relevance for complex models

While climate physics are highly simplified in conceptual models like H07, their dynamics are well understood. The opposite holds true for Earth System Models (ESMs), which resolve multiple processes of climate in detail. To learn if the dynamical mechanism of RFL applies to such a model, we could in theory produce an Arnold tongue diagram as for H07 (Fig. 4). However, since running ESMs is computationally expensive this is presently not possible. Nevertheless, it might be possible to detect signatures of RFL from only few model runs.

First, one should investigate whether the glacial cycles are self-sustained by fixing model parameters at plausible values and fixing the insolation field at its mean value. This is the case if, after a transient time, variations in ice volume on the order of $10$-$100$ kyr persist.

The next step is to sparsely sample an Arnold tongue diagram. First, a ramping parameter must be chosen. This does not have to be a scalar, but can be a function like a parametrisation, as long as its change over time is well defined. The parameter should be one that feasibly could influence the internal period of glacial cycles.

If changing the parameter changes the internal period, then one can compare the average duration of (insolation variation) forced and unforced solutions. If the average duration of the forced solutions is close to either 40 or 80 kyr and remains close even under parameter perturbations that change the internal period, then this is an indication that the system is frequency locked to insolation in a way relevant for the glacial cycles. In that case, there is good reason to research RFL more closely in the model.

Ashkenazy (2006) suggested that synchronisation can be detected by running the system from multiple initial conditions and see if solutions converge. This procedure is not enough for us; we need to know if the internal period can be shifted appropriately with a change in parameter, and we need to know if the durations can robustly cluster on $40$ and $80$ kyr.

8.2 Ramped climate variables

To evaluate whether RFL caused the MPT one must identify slowly changing climate variables. Two such candidate variables are atmospheric $CO_{2}$ and atmospheric or oceanic temperatures. Since less $CO_{2}$ leads to a generally cooler atmosphere, it can be viewed as a proxy for global average atmospheric temperature. Local cooling can occur for other reasons, however.

There are currently no direct measurements of atmospheric $CO_{2}$ across the MPT, but a recent reconstruction back to $-2000$ kyr suggests that the mean $CO_{2}$ did not change in the mean until at least $-1300$ kyr (Hönisch et al., 2009). Since the reconstruction implies that $CO_{2}$ fell $31$ ppm by $-700$, the decrease in $CO_{2}$ must either have been rapid and driving the MPT, or a consequence of it. A rapid change in $CO_{2}$ is still consistent with RFL, but in that case RFL does not explain the abrupt increase in cycle length at the MPT, and instead one must find an explanation for the rapid increase in $CO_{2}$. However, the planned European BEOIC deep ice core drilling in Antarctica can hopefully improve estimates of $CO_{2}$ across the MPT.

There is evidence of a gradual deep ocean cooling since the onset of northern Hemisphere glaciation $2.7$ Million years ago (Lisiecki and Raymo, 2005). How much of this cooling occurred across the MPT is not known, however. The reconstruction of deep water temperatures by (Elderfield et al., 2012) indicates a gradual cooling in the mean from $-1300$ kyr until present, but also a puzzling warming from $-1500$ kyr to $-1300$ kyr. Therefore, glacial cycle length does not appear to have a direct relation with mean deep ocean temperature. However, it may be that sea surface temperatures in the vicinity of major ice sheets are more relevant for glacial dynamics. If so, detailed and reliable reconstruction of such temperatures is necessary to evaluate whether they act as ramped climate variables in RFL.

Another slowly varying parameter could be the erosion of regolith. According to this hypothesis soft material under ice sheets eroded throughout the Pleistocene, enabling them to grow larger before collapsing Clark and Pollard (1998). The hypothesis is difficult to test empirically, however.

In addition to the candidate ramping climate variables mentioned, there may be others that are relevant for the MPT. RFL motivates the search for other such climate variables. These might not only be relevant for RFL, but for any mechanism of the MPT invoking deterministic bifurcation.

9 Criteria for RFL

Having demonstrated RFL in H07 and Paillard and Parennin 2004, we ask in which models RFL is most likely to be relevant.

We expect RFL in all models similar to H07, that is, models with a critical threshold of deglaciation (explicit or not), two intrinsic growth and decay states, additive forcing, and a climate variable that naturally controls the internal period.

Furthermore, RFL is facilitated by dynamics focussed on a single strongly attracting limit cycle. This is because solutions to the system $f$ with ramped parameters then track frozen solutions of $f_{\tau}$ well, and because it is difficult for perturbations to bring solutions away from the neighbourhood of the attractor.

Crucially, a model using RFL needs a parameter that can increase the internal period by $100$ kyr. The models in e.g. (Le Treut and Ghil, 1983) and (Maasch and Salzman, 1990) are therefore difficult to reconcile with RFL since the internal periods are on the order of $10$ and $100$ kyr respectively, and do not change much within the physical range of model parameters.

Lastly, we note that e.g. excitable systems and dissipative resonant oscillators (Crucifix, 2012) also can undergo a rapid change in durations due to frequency locking related phenomena, although they are not self-sustained oscillators. Self-sustained oscillators are distinguished by having an internal period $T_{o}(R)$ through which we can define Arnold tongue diagrams and Devil’s staircases; for non-self-sustained oscillators we have to define these through parameters $R$ directly. Furthermore, it has been argued that the term frequency locking should be restricted to self-sustained oscillators (Pikovsky et al., 2001; Marchionne et al., 2018), why it makes sense to define RFL for self-sustained oscillators only.

10 Conclusions

The glacial cycles did not enter a stationary $100$ kyr world at the MPT; instead, durations between glacial terminations shifted abruptly from approximately ${\sim}40$ to ${\sim}80$ kyr around $-1200$ kyr, followed by a gradual increase (Fig. 1). The dynamical mechanism Ramping with Frequency Locking (RFL) naturally explains this progression of durations. As the internal period of a model glacial cycle model increases gradually, frequency locking to insolation variations causes the durations between glacial terminations to increase sometimes abruptly and sometimes gradually.

The RFL mechanism is rather general and explains the behaviour of a range of models describing glacial cycles and the MPT (Feng and Bailer-Jones, 2015; Huybers, 2007; Crucifix et al., 2011; Paillard, 1998; Paillard and Parrenin, 2004; Tzedakis et al., 2017; Ashkenazy, 2006; Mitsui et al., 2015).

Here we described how RFL can be understood in terms of a dynamical system $f$ and a frozen system $f_{\tau}$ with parameters $R(t)$ fixed at times $t=\tau$. The average duration $\overline{D}_{\tau}$ defined for $f_{\tau}$ provides some information about single durations in solutions $x(t)$ to $f$ around $t=\tau$, but since $\overline{D}_{\tau}$ is defined asymptotically, one must interpret solutions to $f$ in terms of $f_{\tau}$ with care.

Model behaviour can be understood from considering parameter paths through Arnold tongue diagrams and corresponding Devil’s staircases (Fig. 3, 4 and 5). These diagrams as functions of frozen time $\tau$ depend on

• •

the change in average duration $\overline{D}_{\tau}(T_{o})$ as function of $T_{o}$,

• •

the change in internal period $T_{o}(R)$ as function of $R$,

• •

the change in parameters $R(\tau)$ as function of time $\tau$, and

• •

the amplitude $A$ of the forcing.

This decomposition clarifies different ways in which the average duration can change abruptly in models of the class $f$. For instance, the abruptness of the change in $\overline{D}_{\tau}$ can be adjusted either by changing the forcing amplitude $A$ or the ramping of $R(t)$. While the effects of changing $A$ or the ramping of $R(t)$ are typically easy to guess, we are not aware of any general rules dictating the widths of particular Arnold tongues. Such understanding may be researched further.

RFL is relevant also for multi-frequency astronomical forcing. Multi-frequency forcing tends to make Devil’s staircases less abrupt, but durations can still increase rapidly when a model parameter is slowly ramped.

The RFL mechanism provides an explanation for the MPT without the climate system entering a new mode of operation. A shift from ${\sim}40$ kyr long to ${\sim}80$ kyr long cycles due to frequency locking to obliquity and precession, is consistent with data (Fig. 1), and is used in models (Paillard and Parrenin, 2004; Huybers, 2007). This warrants further study of frequency locking characteristics of models throughout the model hierarchy, as well as a search for gradually increasing climate parameters. Some models use a rapidly ramped parameter to accelerate the increase in durations between glacial terminations at the MPT, but such a ramping begs for justification that a model relying solely on frequency locking does not require.

Acknowledgements

We thank Peter Ashwin for valuable discussions.

This research has been funded by the European Union’s Horizon 2020 innovation and research programme for the ITN CRITICS under the Marie Skłodowska-Curie grant agreement No. 643073.

Appendix A Kendall’s tau

Kendall’s tau, here denoted $\tau_{K}$, when testing for monotonicity of a sequence $\{D_{i}\}_{i=1}^{n}$ is defined as

[TABLE]

where $n_{c}-n_{d}=\sum_{i<j}\operatorname{sign}{(D_{j}-D_{i})}$ is the number of pairs $(D_{i},D_{j})$ that are ordered ($D_{j}>D_{i}$) minus the number that is disordered, $n_{0}=n(n-1)/2$ is the total number of pairs and $n_{1}=\sum_{k}t_{k}(t_{k}-1)/2$ is the sum of the number of tied elements $t_{k}$ in the $k$:th group of tied elements. For example, the sequence $\{1,2,2\}$ has two ordered pairs $(1,2)$ and $(1,2)$, zero disordered pairs, and one tied pair $(2,2)$. Hence, $k=1$ such that $n_{c}-n_{d}=2$, $n_{0}=3$ and $n_{1}=1$, giving $\tau_{K}=2/\sqrt{6}\approx 0.82$.

Appendix B Wavelets

Wavelet spectra are estimated with the MATLAB®function cwt, using Morlet basis functions with bandwidth parameter $\omega_{0}=6$ (Torrence and Compo, 1998). Contours show wavelet amplitude (square root of variance) relative to the maximum, incremented in evenly spaced percentage units. The cone of influence marks the $e$-folding time of the amplitude of a discontinuity at the edge of the time interval. Inside the cone of influence edge effects are negligible (Torrence and Compo, 1998).

Appendix C Glacial terminations in LR04

Major glacial terminations in the LR04 stack (Fig. 1) are identified at times $t=$-[$1948$, $1900$, $1863$, $1795$, $1748$, $1708$, $1655$, $1575$, $1535$, $1496$, $1456$, $1412$, $1372$, $1336$, $1290$, $1248$, $1198$, $1126$, $1038$, $964$, $876$, $794$, $718$, $630$, $536$, $434$, $341$, $252$, $140$, $18]$ kyr. We assume conservatively an age model uncertainty with constant standard deviation $6$ kyr over the past $2000$ kyr (Lisiecki and Raymo, 2005), which gives a standard deviation $\sqrt{2}\cdot 6$kyr on the durations between terminations, assuming somewhat wrongly that errors are independent and normally distributed.