The Mid-Pleistocene Transition induced by delayed feedback and bistability
Courtney Quinn, Jan Sieber, Anna S. von der Heydt, Timothy M. Lenton
TThe MPT as a result of bistability
The Mid-Pleistocene Transition induced by delayedfeedback and bistability
Courtney Quinn, a) Jan Sieber, Anna S. von der Heydt, and Timothy M. Lenton College of Engineering, Mathematics and Physical Sciences,University of Exeter, Exeter EX4 4QE, United Kingdom Institute for Marine and Atmospheric Research, Centre for Extreme Matter andEmergent Phenomena, Utrecht University, Princetonplein 5, 3584 CC Utrecht,The Netherlands Earth System Science, College of Life and Environmental Sciences,University of Exeter, Exeter EX4 4QE, United Kingdom (Dated: 12 November 2018)
The Mid-Pleistocene Transition, the shift from 41 kyr to 100 kyr glacial-interglacial cycles that occurred roughly 1 Myr ago, is often considered as achange in internal climate dynamics. Here we revisit the model of Quaternaryclimate dynamics that was proposed by Saltzman and Maasch (‘Carbon cycleinstability as a cause of the late Pleistocene ice age oscillations: modellingthe asymmetric response’. Glob Biogeochem Cycle 1988; 2: 177-185-fromthis point, referred to as SM88). We show that it is quantitatively similarto a scalar equation for the ice dynamics only when combining the remainingcomponents into a single delayed feedback term. The delay is the sum of theinternal time scales of ocean transport and ice sheet dynamics, which is onthe order of 10 kyr.We find that, in the absence of astronomical forcing, the delayed feedbackleads to bistable behaviour, where stable large-amplitude oscillations of icevolume and an equilibrium coexist over a large range of values for the delay.We then apply astronomical forcing using the forcing data for integrated sum-mer insolation at 65 degrees North from Huybers and Eisenman (
IntegratedSummer Insolation Calculations.
NOAA/NCDC Paleoclimatology ProgramData Contribution δ O records’.
Paleoceanography a) [email protected] a r X i v : . [ phy s i c s . a o - ph ] N ov he MPT as a result of bistability 2 time (kyr BP) δ O ( pe r m il ) FIG. 1: Compilation of 57 globally distributed benthic δ O records for the last 2 Myr(Lisiecki and Raymo, 2005).
I. INTRODUCTION
The Pleistocene is characterised by its successive glacial-interglacial cycles whose docu-mented periodicities has been a subject of interest for decades. Hays, Imbrie, and Shackle-ton (1976); Paillard (1998), and Ganopolski and Calov (2011) have all presented evidencethat orbital forcing plays a large role in the onset and termination of glacial periods, butthe full extent of the variations cannot be explained solely by forcing (Ridgwell, Watson,and Raymo, 1999; Shackleton, 2000; Paillard, 2001). One actively researched aspectof the Pleistocene is the abrupt change in frequency of major glaciations known as theMid-Pleistocene Transition (MPT) (Pisias and Moore, 1981; Maasch, 1988; Maasch andSaltzman, 1990; Paillard, 1998; Clark et al. , 2006; Crucifix, 2012; Ashwin and Ditlevsen,2015). Fig. 1 shows a proxy record compiled by Lisiecki and Raymo (2005) for globaltemperature taken from ocean core sediments which can be directly related to global icecover. Spectral analysis has been performed on this time series, and it has been observedthat the signal of the 100-kyr cycle began to rise 1250 kyr BP and was established as thedominant cycle by 700 kyr BP (Clark et al. , 2006; Lisiecki and Raymo, 2007; Dijkstra,2013). Many researchers have tried to determine not only the cause of this switch, butalso the driving force behind the 100 kyr cycles themselves (Maasch and Saltzman, 1990;Saltzman and Maasch, 1991; Paillard, 1998; Gildor and Tziperman, 2001; Paillard andParrenin, 2004; Ganopolski and Calov, 2011; Ashwin and Ditlevsen, 2015). While the41-kyr oscillations have been attributed to external forcing (Paillard, 2001), the 100-kyrcycles have been proposed to be a result of nonlinear responses of the climate system itself(Imbrie et al. , 1993; Yiou et al. , 1994; Gildor and Tziperman, 2001).The major glacial-interglacial cycles were originally attributted to changes in insolationdue to orbital variations (i.e. the Milankovitch cycles). Milankovitch argued that summerinsolation at high northern latitudes determined the main pacing of glaciations, as thesechanges control the length of seasons and amount of solar energy being recieved in highlatitudes where the ice resides. The most prominent frequencies of these changes are 19and 23 kyr due to precession and 41 kyr due to obliquity (Paillard, 2001). Precession isthe orientation of the rotational axis of the earth (axial) and the rotation of the orbitalaxis (apsidal), and obliquity is the angle between the rotational axis and the orbital axis.Milankovitch (1941) argued that the main driving forces of climate fluctuations are obliq-uity and precession, which agrees with the records prior to the mid-Pleistocene transition.There is a smaller signal of orbital forcing that correlates with the 100-kyr oscillations (Im-brie et al. , 1993). Eccentricity, the amount the earth’s orbital ellipse deviates from a circle,shifts with a main period of 413 kyr, but there are components that vary with periodsbetween 95 kyr and 125 kyr. The amplitude of these forcing signal bands, however, are anorder of magnitude smaller than the signals of the 23- and 41-kyr bands (Hays, Imbrie, andShackleton, 1976).Since Milankovitch’s theory, there have been many other attempts to identify the causeof the MPT (see review paper by Crucifix (2012)). The starting point of our analysis isthe collection of Saltzman models from the late twentieth century. They were a collectionof attempts to model the transition as a bifurcation in low-order dynamical systems due toclimate feedbacks. The models typically only had three dynamic variables usually repre-senting some measure of global ice volume, global mean temperature, and ocean circulationstrength. Other modelling attempts, after the Saltzman models, include (in order of in-he MPT as a result of bistability 3creasing complexity) conceptual models (Paillard, 2017), threshold models (Tzedakis et al. ,2017), nonsmooth models (Paillard, 1998), relaxation oscillators (Ashwin and Ditlevsen,2015), box models of intermediate complexity (Gildor and Tziperman, 2001; Tzipermanand Gildor, 2003), and Earth system models of intermediate complexity (EMICs) (Ganopol-ski and Calov, 2011). While all of these studies were able to reproduce a transition similarto the MPT, they are all based on the hypothesis that the MPT occured because of somesignificant change in the dynamics of the system. Huybers (2009) proposed an alternativehypothesis that the MPT was a spontaneous response to the astronomical forcing. In thisstudy we focus on this alternative hypothesis, and provide more evidence using a previouslyestablished model by Saltzman and Maasch (1988).The paper is arranged as follows. Section 2 revisits and analyses the three-variable Saltz-man and Maasch 1988 model. We observe that two of the variables mostly act as delaysin feedback loops such that we may reduce the model to a scalar equation with delayedfeedback. We also identify a region of bistability at the parameter values of interest whichwas not explored in the original model analysis. Section 3 shows the model response in thebistable region when astronomical forcing is introduced, which includes MPT-like realisa-tions with no change in parameters. In section 4 we investigate the sensitivity of the model.We observe that the system forgets its initial condition and has low sensitivity whenever itexhibits a small-amplitude oscillatory response. It is much more sensitive to disturbanceswhen exhibiting a large-amplitude oscillatory response. This sensitivity affects mostly thetiming of major deglaciation events. Section 5 discusses our results.
II. ADAPTATION OF THE SALTZMAN AND MAASCH MODEL
We begin by considering the original model of Saltzman and Maasch (1988), which wewill refer to as SM88, ˙ X ( t ) = − X ( t ) − Y ( t ) , (1a)˙ Y ( t ) = − pZ ( t ) + rY ( t ) + sZ ( t ) − Z ( t ) Y ( t ) , (1b)˙ Z ( t ) = q ( − X ( t ) − Z ( t )) . (1c)Here X, Y, Z are scaled versions of global ice mass, atmospheric CO , and North AtlanticDeep Water (NADW) respectively. More precisely, the variables represent anomalies froma background state. The first term in (1a) represents the feedbacks of global ice mass.The combined effect of damping and negative feedback from NADW production is takento outweigh the positive feedback from ablation (melting and loss of ice through icebergs).The second term in (1a) is the direct effect of higher temperatures (caused by higher CO levels) leading to loss of ice mass. In equation (1c) the negative sign of ice anomaly X ismotivated by the reduction of NADW production with loss of ice mass. This equation has atime scaling q which is the ratio of ice sheet time constant (10 kyr, which equals t = 1 in thescaling of (1a)-(1c)) to the response time of the deep ocean. The equation for atmosphericCO (1b) has a negative term − pZ for the negative effect due to NADW- more NADW leadsto more ventilation of the deep ocean (stronger mixing between surface and deep waters) andan increased downdraw of atmospheric CO into the deep ocean. The term rY accounts forthe positive feedbacks of atmospheric CO , related to changes in sea surface temperatures,sea ice extent, and sea level, which outweigh negative feedbacks. The nonlinearity ( s − Y ) Z is motivated by the appearance of locally enhanced instabilities in the Southern Ocean dueto increased NADW production. The NADW meets colder, denser Antarctic Bottom Waterand induces vertical mixing. This brings CO -rich water to the surface which releases CO into the atmosphere depending on the current level of atmospheric CO ( Y ). Each of theseeffects (except damping of the nonlinear terms) have an associated parameter which controlstheir relative strengths. These were the feedbacks between global ice mass, atmosphericCO , and NADW deemed most important by Saltzman and Maasch (1988). Present daystudies continue to stress the interaction between these three variables as being essentialfor the MPT, including a study by Chalk et al. (2017) which uses a geochemical box modelhe MPT as a result of bistability 4to confirm that carbon cycle feedbacks related to ocean circulation are necessary to sustainthe long glacial cycles in the late Pleistoncene.A change of variable V ( t ) = − X ( t ), where V is the negative of the global ice massperturbations, replaces the two linear equations (1a) for X and (1c) for Z by a chain oflinear first-order filters˙ V ( t ) = Y ( t ) − V ( t ) , ˙ Z ( t ) = q [ V ( t ) − Z ( t )] . (2)A linear chain of filters is well known to approximate a delay (Smith, 2010): Z ( t ) ≈ Y ( t − τ ), where τ = 1 + 1 q . (3)Appendix A gives further comments on the linear chain approximation for delays. Usingthis connection between linear chains and delays, we can rewrite the model as a scalardelay-differential equation (DDE) for Y . To facilitate the comparison to data at a laterstage, we shift time and exploit that in the DDE, the ice mass anomaly V = − X is justthe delayed (by time 1) CO anomaly: X ( t ) = − Y ( t − X :˙ X ( t ) = rX ( t ) − pX ( t − τ ) − X ( t − τ ) [ s + X ( t )]. (4)In DDE (4) one of the parameters r , p , s , τ is redundant and could be removed by a rescalingof X and time. Thus, (4) is objectively simpler than the original SM88 model (1) as it hasfewer free parameters. For the purpose of model comparison, the DDE system will be leftas above. A. Dynamics of DDE, compared to SM88
We initially perform a systematic analysis of DDE (4) and compare it to the originalSM88. We find that models (1) and (4) have the same long-time behaviour (using therelation q = τ − implied by (3)).Since delays, whether in the form of chains as in SM88 or explicit, don’t affect locationof equilibria, i.e. X ( t ) = X ( t − τ ) (or X ( t ) = − Y ( t ) = − Z ( t )), equilibria for both models,(1) and (4), (called X eq ,j ) satisfy0 = − pX eq + rX eq − sX − X . (5)Thus, we have a trivial equilibrium X eq , = 0 for all parameter values and non-zero equilibria X eq , , for s > p − r ): X eq , = 12 (cid:104) − s ± (cid:112) s − p − r ) (cid:105) . (6)The stability of the equilibria may differ, however, since eigenvalues of the linearisation maycross the imaginary axis (leading to oscillations in a Hopf bifurcation (Kuznetsov, 2013))at different parameter values for (1) and (4).All other trajectories have to be studied numerically (for both, SM88 (1) and DDE(4)).Since the scaling between the two models are the exactly same, we will use values for p , r ,and s from Saltzman and Maasch (1988) for all numerical studies: p = 0 . r = 0 . s = 0 .
8, (7)while varying the timescale (or delay) of the feedback processes ( τ for DDE (4), q for SM88model (1)). A value of q = 1 . τ ≈ . q was q ≥ q was defined as the ratio of ice sheet time constant,10 kyr, to the slow response time of deep ocean). This would put realistic values of τ inthe range [1,2]. The value of q by Saltzman and Maasch (1988) would approximate a deephe MPT as a result of bistability 5 τ -4-202 X stable equilibriumunstable equilibriumsaddle equilibriumstable periodic orbitunstable peridic orbitHopf bifurcationHomoclinic bifurcation r ee r es r e r el r l (a) DDE model (4)The dotted black lines indicate values of τ used in section III: τ ref = 1 . τ = 1 . τ = 1 . τ = 1 . τ -4-202 X (b) SM88 ( τ = q + 1) FIG. 2: Bifurcation diagrams of both models for delay parameter τ . Other parameters: p = 0 . r = 0 . s = 0 . τ < . τ = 1 .
83 in SM88. The comments below(4) also show that the change τ → τ /α , p → αp , r → αr , s → √ αs leaves the systemunchanged such that all phenomena reported in this paper are observable for any delay andappropriately rescaled parameters.Figure 2 shows the bifurcation diagrams for varying τ in DDE (4) and SM88 model (1),respectively. Although the locations of bifurcations are slightly shifted, qualitatively thetwo figures agree nicely. In particular, even though the space of possible initial condi-tions for DDE (4) is infinite-dimensional (every possible history on [ − τ,
0] gives a differenttrajectory), the long-time behaviour of trajectories in the ( X ( t ) , X ( t − τ ))-plane follows atwo-dimensional ordinary differential equation (ODE) in the range we explored. The dia-grams in Fig. 2 are partitioned into five main regions of different global behaviour separatedby grey vertical lines can be seen in both figures. From left to right the attractors in eachregion are: • [ r ee ] two stable equilibria, • [ r es ] one stable equilibrium and one stable small-amplitude periodic solution, • [ r e ] one stable equilibrium, • [ r el ] one stable equilibrium and one stable large-amplitude periodic solution, and • [ r l ] one stable large-amplitude periodic solution.The region of interest in this study is the bistable region with large-amplitude periodicorbits [ r el ]. Note that [ r el ] can be split further into two sections: one containing unstableperiodic orbits and the other not. The effect of the unstable periodic orbits is a drasticallyreduced basin of attraction for the stable equilibrium. However, since there is no change inthe attractors, we will consider these as one region. Region [ r el ] was not discussed in a laterhe MPT as a result of bistability 6 time (kyr) -3-2-101 X ( t ) -2.5 -2 -1.5 -1 -0.5 0 0.5 1 X(t) -2.5-2-1.5-1-0.500.51 X ( t - τ ) FIG. 3: Example trajectories in the bistable region [ r el ] for τ = 1 . X ( t ) = − .
05 (blue)and X ( t ) = 0 .
05 (red) for t <
0. Top- time profiles, bottom - phase portraits.parameter study of the original model (Maasch and Saltzman, 1990). For the remainder ofthis paper we will discuss only the DDE model, but similar results can be seen for the ODEmodel.
B. The bistable region [ r el ] The focus of this study is the bistability seen for τ ∈ [1 . , . τ very close to the [ r e - r el ] boundary where theperiod approaches infinity). This cycle length agrees with what is seen in the data (evenmore so when one adds the external forcing; see Sec. III). In previous studies of SM88model (1), a transition between two stable states was enforced by a parameter shift througha Hopf bifurcation (Saltzman and Maasch, 1988; Maasch and Saltzman, 1990). We willdemonstrate in the next section that the bistability in region [ r el ] makes transitions betweenthe two states possible without any change of parameter when the model is subjected toexternal forcing. III. THE MID-PLEISTOCENE TRANSITION UNDER THE FORCED MODEL
In this section we show that when adding astronomical forcing to model (4), a transitiontypically occurs at the same time as the MPT is seen in recorded data without furtherparameter tuning. In section III A we describe the astronomical forcing considered and howwe include it in the model. We then examine the responses for different delays and differentforcing strengths in sections III B and III C.
A. Astronomical forcing
Hays, Imbrie, and Shackleton (1976) have provided evidence that the glacial cycles duringthe Pleistocene are driven primarily by variations in the earth’s orbital cycle. This includeshe MPT as a result of bistability 7FIG. 4: Normalised integrated July insolation at 65 ◦ N, adapted from Huybers andEisenman (2006).changes in precession (orientation of the rotational axis), obliquity (angle between the ro-tational axis and orbital axis), and eccentricity (orbital ellipse’s deviation from a circle).These modes vary approximately periodically, with cycle lengths of 19/23 kyr, 41 kyr, and100 kyr respectively (Milankovitch, 1941; Berger, 1978; Huybers and Eisenman, 2006).Fig. 4 shows a time series of average daily summer insolation at 65 ◦ N computed by Huybersand Eisenman (2006) based on the model introduced by Huybers (2006). To investigate theeffect of this forcing on DDE (4), we include the astronomical forcing in the same way asin Maasch and Saltzman (1990). We add forcing signal M ( t ) shown in Fig. 4 with negativeamplitude u , ˙ X ( t ) = rX ( t ) − pX ( t − τ ) − X ( t − τ ) [ s + X ( t )] − uM ( t ). (8)The precise procedure for extracting M ( t ) from the publicly available data source can befound in Appendix B. We are interested in how the bistable region responds to this externalforcing in dependence of two parameters: the delay τ and the forcing amplitude u .For small u the system is expected to exhibit two types of responses to forcing in thebistable region. Each of them is a perturbation of an attractor of the unforced system,namely the equilibrium and the large-amplitude periodic orbit, which persist for small u .We will refer to these responses as the small-amplitude and the large-amplitude response(compare red and blue time profiles in Fig. 5c). For increasing u we expect to observeincreasingly frequent transitions between these responses. For large forcing amplitudes u the internal dynamics of the model will be dominated by the forcing. B. Responses for different delays in the bistable region
We choose a moderate value of the forcing amplitude u = 0 .
25 and investigate the responseat four values of τ , labelled in Fig. 2a by τ ref , τ , τ , and τ . The first value τ ref = 1 .
25 isoutside of the region of bistability and is used as a reference trajectory to which the solutiontrajectories for the other delays τ , τ and τ are compared. The results are shown in Fig. 5(response for τ ref shown in red, for other delays in blue)The first comparison is made close to the boundary [ r e ]–[ r el ] at τ = 1 .
3. We see infigure 5a that both trajectories change in synchrony, exhibiting only the small-amplituderesponse. In the middle of the bistable region ( τ = 1 .
45, shown in Fig. 5b), the solutionshows a small-amplitude response in most of the first half of the time window and a large-amplitude response in the second half. Close to the Hopf bifurcation of the autonomousstable equilibrium ( τ = 1 .
6, shown in Fig. 5c) the solution exhibits primarily a large-amplitude response. This is expected due to the weakening attraction of the autonomousstable equilibrium and its shrinking basin of attraction. We will demonstrate in the followingsection that these transitions generate dynamics with time profiles that are qualitativelysimilar to the records of the MPT.
C. Variable forcing strength
For the exploration of the effect of different forcing strengths, we take the same values for τ as in Fig. 5, covering the range of the bistable region. For each τ we compute trajectorieshe MPT as a result of bistability 8 (a) τ = τ (blue), τ = τ ref (red)(b) τ = τ (blue), τ = τ ref (red)(c) τ = τ (blue), τ = τ ref (red) FIG. 5: Trajectories in bistable region (blue) compared to a reference trajectory at τ ref = 1 .
25 (red) for delays (a) τ = 1 .
30, close to homoclinic connection, (b) τ = 1 . τ = 1 .
60, close to Hopf bifurcation. Values of τ areindicated on bifurcation diagram in Fig. 2a. Other parameters: p = 0 . r = 0 . s = 0 . u = 0 .
25. Initial condition X ( t ) = − . t ∈ [2 + τ, (a) (b) (c) FIG. 6: Distance from reference solution at τ ref = 1 .
25, for different forcing strengths u .Averages taken over window length of size τ . (a) τ = 1 .
30, close to homoclinic connection.(b) τ = 1 .
45, middle of bistable region. (c) τ = 1 .
60, close to Hopf bifurcation. Initialcondition X ( t ) = − . t ∈ [2 + τ, in all cases. for different u , ranging from u = 0 to u = 0 .
6, with the initial history X ( t ) = − . t ∈ [2 + τ, τ ref = 1 .
25 and the same initial condition and forcing strength u ,averaged over a window of length τ . We note that differences to a reference trajectory at u = 0 and identical delay τ (using same initial condition) give qualitatively similar results.Remarkably, figures 6a and 6b show that there is a distinct period around 700-800 kyr BPwhere the solutions diverge from the reference trajectory in a large range of forcing strengths u . This suggests that some aspect of the forcing around this time kicks the trajectories intothe basin of attraction of the large-amplitude response. An additional area like this is seenin figure 6b near 1600 kyr BP.The second interesting feature in all three figures is a threshold behaviour in the forcingstrength u : below a certain value of u specific to each τ , transitions do not occur such thathe MPT as a result of bistability 9 time (kyr BP) -3-2-101 X FIG. 7: Example trajectory that exhibits an MPT-like transition with τ = 1 .
45 and u = 0 . τ we observe an additional larger thresholdvalue for the transition at 1600 kyr BP.Thresholds in u are also seen when periodic forcing is added to the model instead of theinsolation time series. We consider a new definition of the forcing function, M ( t ) = sin (cid:18) π . t + π (cid:19) . (9)This forcing has a period of 41 kyr. If u >
0, equation (8) has two stable solutions: asinusoidal signal with a 41 kyr period, or a large-amplitude asymmetrical quasiperiodicresponse. When recreating figure 6b with this periodic forcing, we see the same thresholdtransition at u = 0 .
08. We can show that this transition is due to moving basins of attractionfor the two stable solutions. The other sharp transition we saw in Fig. 6 would not bepossible with this forcing. Either the solution would start to transition immediately, or itwould not transition at all. This is also the case if an envelope of modulated amplitudedescribing changes of eccentricity at 1 . IV. MODEL SENSITIVITY TO NOISE: DESYNCHRONISATION AND INCREASEDROBUSTNESS OF TRANSITIONA. Finite-time Lyapunov Exponents
To analyse how trajectories depend on their history at specific instances of the forcing,we compute finite-time Lyapunov exponents (FTLEs) using a QR decomposition method(details in Appendix C). This method was previous used in De Saedeleer, Crucifix, andWieczorek (2013) to illustrate the desynchronisation of nearby trajectories in a van derPol-type oscillator model. We will apply the ideas presentated in that study to our forcedmodel.The difference between FTLEs and classical Lyapunov exponents is that FTLEs arerecorded for a family of time windows of finite length w rather than over the entire long timerun. Thus, FTLEs are time-dependent functions instead of real numbers. A positive FTLEalong a given trajectory X ( t ) at time t indicates that some nearby trajectories divergeexponentially from X ( t ) in the time window [ t − w, t ]. Therefore, X ( t ) is sensitive to smallperturbations in the time window [ t − w, t ]. While a trajectory could be asymptoticallystable, a positive FTLE at a time t indicates temporary amplification of perturbationsfrom the attracting trajectory (observed as temporary desynchronisation, see (De Saedeleer,Crucifix, and Wieczorek, 2013)). B. FTLE implications on MPT and timing of major deglaciations
We analyse one trajectory, showing a forcing-induced MPT ( τ = 1 . u = 0 . time (kyr BP) -202 X -0.200.2 FT L E FIG. 8: Finite-time Lyapunov exponents for window length w = 250 kyr computed frommodel run with τ = 1 .
45 and u = 0 . X ( t ) (dotted blue), FTLE([ t − w, t ]) (gold).FIG. 9: Trajectories of DDE model with noise: u = 0 . σ = 0 . τ = 1 . w = 250 kyr. This window length is chosen in order to filter outthe dominant frequencies of the forcing. Similar results are seen with any window lengths w from 150-500 kyr.Fig. 8 shows the time profile of the largest FTLE along the trajectory. Before thetransition around 800 kyr BP, the FTLE generally remains negative apart from a few shortexcusions above zero. At 1000 kyr BP the FTLE approaches zero and remains there forsome time. Just before 800 kyr BP it goes positive and then on average stays positive forthe remainder of the trajectory.The negative FTLEs leading up to the transition confirm that the trajectory forgetsits initial history and the effect of past disturbances. In particular, this implies that theinfinite-dimensional nature of the DDE’s possible initial conditions does not play a role forthe MPT. Whenever the system is exhibiting the small-amplitude response we observednegative FTLEs.The positive FTLEs indicate a sensitivity of the trajectory during the large-amplituderesponse. This sensitivity affects the precise timing of the deglaciations, as we now demon-strate by noise-induced desynchronisation of nearby trajectories.To explore further the desynchronisation phenomenon outlined in De Saedeleer, Crucifix,and Wieczorek (2013), we add noise to our system. The stochastic delay differential equation(SDDE) is then given as dX ( t ) =[ − pX ( t − τ ) + rX ( t ) − sX ( t − τ ) − X ( t − τ ) X ( t ) − uM ( t )] dt + σdW ( t ) , (10)Here, W ( t ) is standard white noise and σ is the noise amplitude. We will always consider σ as a fraction of the forcing amplitude u , i.e. σ = u .We set the deterministic forcing strength u = 0 .
15 and the noise amplitude σ = 0 . X ( t ) , X ( t − τ ))-plane along the unforced periodic orbit, and in Fig. 11. Figure 11shows the transition from synchronization to desynchronization in a density plot. At ap-proximately 800 kyr BP (indicated by the gray line) the probability density makes a sharptransition from a small-variance to a large-variance-low-maximum density. The timing ofhe MPT as a result of bistability 11 -2-101 -2-101 X ( t - τ ) -2 -1 0 1-2-101
400 kyr BP -2 -1 0 1
X(t)
200 kyr BP -2 -1 0 1
FIG. 10: Phase portraits at specific time instances of 500 model runs with noise: u = 0 . σ = 0 . τ = 1 .
45. The red curve represents the unforced periodic orbit.FIG. 11: Probability density in time (top) and its standard deviation (bottom) of 500model runs with noise: u = 0 . σ = 0 . τ = 1 .
45. The gray line indicates 800 kyr BPthis transition agrees with the first large positive excursion of the FTLE in the deterministiccase and with the MPT. The standard deviation shows this transition as well, but with a lagof about 50 kyr. Figures 10 and 11 also show that the distribution by desynchronisation isfar from uniform. There are still distinct deglaciation times favoured for many realisations.For this reason the desynchronisation observed in figures 10 and 11 does not contradictstudies that have argued for phase-locking to different components of the orbital forcing.One of the original hypotheses of phase-locking in the late Pleistocene was suggested to berelated to eccentricity, specifically associated to events of low eccentricity (Hays, Imbrie,and Shackleton, 1976; Paillard, 2015). More recent studies have looked at the possibilityof locking to precession or obliquity. In Ridgwell, Watson, and Raymo (1999) the authorsargue locking to every 4th or 5th precession cycle, while Huybers and Wunsch (2005) arguefor locking to every 2nd or 3rd obliquity cycle. Later, Huybers (2011) attributes locking toa combination of precession and obliquity, and states strongly that the pacing cannot beattributed solely to one of these components. Robust evidence for locking requires testingon a range of initial conditions or noise realizations (or long time series), but the short datarecord corresponds to only one realization of noise disturbance and one initial condition.Thus, data available may be insufficient to distinguish the higher-order locking proposed inthe literature from the level of desynchronization we report in figures 10 and 11.he MPT as a result of bistability 12FIG. 12: Trajectories of DDE model with noise: u = 0 . σ = 0 . τ = 1 . τ ref = 1 .
25, for different forcing strengths u and noise strength σ = u/
30. Averages taken over window length of size τ = 1 .
45. Initialcondition X ( t ) = 0 . t ∈ [2 + τ,
2] Myr BP.
C. Effect of noise with stronger forcing
Figure 12 shows that the addition of noise enhances the MPT for strong forcing. Whenadding noise, we observe transitions for forcing strengths u for which there was no transitionin the deterministic case (compare Figures 6b and 13). In Fig. 6b the last persistenttransition occurs at u = 0 .
33. Fixing the forcing strength at u = 0 .
45, we compute tenrealisations. Here, two realisations exhibit transitions that persist until the end of the run.Although the transitions don’t appear as commonly as in the weaker forcing case, with noiseit is possible to observe an MPT-like transition across a wider range of forcing amplitudes.Figure 13 shows a systematic overview of the transition enhancing effect of noise. Weadded noise of amplitude σ = u and compute the distance diagram as in Fig 6b. Weobserve transitions occuring above the maximal value of u for which transitions occured inthe deterministic case. V. DISCUSSION
We have presented an analysis of the well known Saltzman and Maasch (1988) model ofQuaternary climate dynamics. We have shown that it can be reduced to a scalar equationfor the ice mass anomaly X ( t ), by collecting all other variables into a delayed feedback term.Through this formulation we are able to explore the dependence on the time scale τ of thedelay in this feedback. We discovered a region in which a stable equilibrium and stablelarge-amplitude periodic orbit coexist for the unforced system. This parameter region ismore consistent with physically relevant timescales for the ice sheet and deep ocean responsethan those explored in SM88.We observe two different responses to astronomical forcing of the type propsed by Huy-bers and Eisenman (2006): small-amplitude response (tracking the equilibrium) and large-amplitude response (tracking the periodic orbit) In the deterministic bistable region themodel exhibits a switching behaviour between these two responses. The timings of theswitches are robust for different values of the delay and focing strength. The most promi-he MPT as a result of bistability 13nent switch occurs around 800 kyr BP which is within the range of when the Mid-PleistoceneTransition is believed to have occurred.One trajectory showing the MPT is analysed using finite-time Lyapunov exponents. Usinga window of 250 kyr we observe that the largest FTLE is negative up to 1000 kyr BP,remains around zero during the transition, and is mostly positive after the transition. Wedemonstrate that this positivity induces phase desynchronisation which is then confirmedby adding noise to the model. The added noise also allows for the transition to be seenwith stronger forcing.The novelty of this study is in the observation that transitions similar to the MPT occurrobustly without any additional climate forcing, such as a nonstationary background state.The transition arises from the bistability rather than a shift of a system parameter througha bifurcation. We also highlight that with noise the transition is possible but doesn’talways occur, and that the phase (or timings of the deglaciations) are variable (howeverstill influenced by the forcing frequencies). While our demonstration was done for thespecific low-order model SM88, we conjecture many of the effects (instability caused bydelayed feedback, transition induced but not caused by the forcing) to be present in othermodels with similar bistability.Many models assume a slowly varying background state, typically related to slowly de-clining CO levels, and attribute this as the primary mechanism for the MPT (Saltzmanand Maasch, 1991; Tziperman and Gildor, 2003; Ashwin and Ditlevsen, 2015; Tzedakis et al. , 2017; Paillard, 2017). Our proposed mechanism does not exclude the possibility ofthis changing background state. However, it does not require the slow change of backgroundstate to cause a bifurcation. Instead, we consider a new null hypothesis for the MPT whichinvolves only bistability and external (orbital and/or stochastic) forcing. Future work couldinclude an additional varying background state within the parameter regions mapped inFig. 2a and analyse the effects on the transition.We reported results for one specific astronomical forcing time series (65 ◦ N summer inte-grated insolation), which is a weighted average version of the typical insolation thresholdforcing (see Appendix B). Qualitatively identical results can also be found when choosinga particular insolation threshold. However the particular threshold at which the behaviouris observed changes for different delays. The weighted average over thresholds from Ap-pendix B does not depend on an additional parameter (forcing threshold) and the dynamicalbehaviour is preserved.Some questions that arise from this study will need to be addessed by deeper mathematicalanalysis. Why does the model prefer certain times to switch, and what causes the sharpthreshold in forcing strength where switching becomes possible, visible in Fig. 6? The lattercan be answered by considering simple harmonic forcing with the dominant period of 41 kyrseen in the astronomical forcing. This threshold should be analysed further in a follow-upstudy. Another question to consider is the effect of human activities on the system. Is itpossible anthropogenic contributions could knock the system out of the oscillations again,and perhaps to a completely different dynamical region? Addressing this question wouldbe useful for understanding current and future climate scenarios.
Appendix A: The linear chain approximation
A linear chain of first-order ODEs approximates a delay in the following sense (see Smith(2010) for a precise derivation). Here we adapt the derivation to our particular system.The 3-dimensional SM88 model (1), using (2), has a 2-dimensional linear chain and is ofthe following form: d Y d t = F ( Y ( t ) , Z ( t )) , (A1a)d V d t = ( Y ( t ) − V ( t )) , (A1b)d Z d t = q ( V ( t ) − Z ( t )) . (A1c)he MPT as a result of bistability 14Extracting the linear chain system we haved (cid:126)ydt = A(cid:126)y ( t ) + (cid:126)b ( t ) , (A2)where (cid:126)y = [ V, Z ] T . We have A and (cid:126)b ( t ) as follows: A = (cid:20) − q − q (cid:21) , (cid:126)b ( t ) = (cid:20) y ( t )0 (cid:21) . (A3)Since A is invertible, the fundamental solution matrix of (A2) is then given byΦ( t ) = P e Λ t P − (A4)where P is the matrix of eigenvectors associated with eigenmatrix Λ,Λ = (cid:20) − − q (cid:21) , (A5)We have the following basis of A and it’s inverse, P = (cid:20) q − q
01 1 (cid:21) , P − = (cid:20) qq − − qq − (cid:21) . (A6)From Perko (2013) we have the fundamental solution (cid:126)y ( t ) = Φ( t ) (cid:126)y (0) + (cid:90) t Φ( t )Φ − ( s ) (cid:126)b ( s )d s. (A7)Assuming the solution has existed arbitrarily far in the past and using a change of variableswe can extract an equation for Z ( t ), Z ( t ) = (cid:90) ∞ Y ( t − τ ) K ( τ )d τ , (A8)where the kernel in (A8) is K ( τ ) = qq − (cid:2) e − τ − e − qτ (cid:3) . (A9)Approximating the true distributed delay kernel K in (A8) with a delta distribution at itsexpected value E ( K ( τ )) = (cid:90) ∞ qq − (cid:2) e − τ − e − qτ (cid:3) τ d τ = q + 1 q (A10)we approximate (A1) by a DDE with a discrete delay dYdt = F (cid:16) Y ( t ) , Y (cid:16) t − q + 1 q (cid:17)(cid:17) . (A11) Appendix B: Extraction of integrated summer insolation forcing M ( t ) from data We use the Integrated Summer Insolation Calculations data set provided by Huybers andEisenman (2006) (found at ), last accessed:17 October 2018. The data is provided as daily average summer energy in GJ/m calculatedfrom number of days the insolation was above a given threshold (W/m ). We used thedata provided for the latitude of 65 ◦ North. Rather than selecting a particular threshold,he MPT as a result of bistability 15we took the average of all the thresholds. This has the effect of giving a linear weightingto the threshold reached in a day: days that reached higher insolation levels are given aproportionally larger weight. The precise expression is as follows, I agg ( t ) = (cid:88) i =1 I i ( t ) . (B1)Here I i ( t ) is the average daily summer insolation calculated by Huybers and Eisenman(2006) for the threshold 25( i −
1) in year t (column number i + 1 in the data file). We thennormalise this data through M ( t ) = I agg ( t ) − mean t I agg std t I agg , (B2)where mean t I agg and std t I agg are the mean and standard deviation of I agg ( t ) over all times t . Appendix C: Computation of Finite-Time Lyapunov Exponents (FTLEs)
Our algorithm for computing Lyapunov exponents for a delay differential equation (DDE)follows from Farmer (1982). a. Linearisation
We consider a trajectory X ( t ) of (8) and linearize (8) along this tra-jectory: ˙ x ( t ) = J ( t ) x ( t ) + J τ ( t ) x ( t − τ ), (C1)where J an J τ are the derivatives of the right-hand side of (8) with respect to X ( t ) and X ( t − τ ): J ( t )= r − X ( t − τ ) , (C2) J τ ( t )= − p − sX ( t − τ ) − X ( t ) X ( t − τ ) . (C3) b. Discretisation We consider an approximate solution of (C1) x i = x ( t i ) for t i = t + ih , where h is small step size, obtained using a second order trapezoidal numbericalsolver for DDEs. We then discretize (C1) into m = τ /h steps, (cid:126)y i +1 = M i (cid:126)y i ; (cid:126)y i = [ x i , ..., x i − m ] T , (C4)where M i is a square ( m + 1) × ( m + 1) matrix of the form M i = M i, . . . M i,m M i,m +1 . . . . . . with entries M i, =1 + h J ( t i +1 ) + J ( t i )) + h J ( t i +1 ) J ( t i ) ,M i,m = h J τ ( t i +1 ) ,M i,m +1 = h J τ ( t i ) + h J ( t i +1 ) J τ ( t i ) . he MPT as a result of bistability 16 c. QR method We use a continuous QR algorithm for computing the Lyapunov ex-ponents (see Dieci, Russell, and Van Vleck (1997)). The QR algorithm is based on thenumerical linear algebra factorization of a matrix M into an orthogonal matrix Q and anupper triangular matrix R . Although DDEs have an infinite number of Lyapunov expo-nents, this method allows us to compute up to m + 1 through the discretisation. We canchoose to compute only the largest l Lyapunov exponents by using a non-square Q . Weinitialise an arbitrary othogonal Q as an m × l matrix of form Q = (cid:20) I l ( m − l ) × l (cid:21) . Note that I l is the l × l identity matrix and 0 ( m − l ) × l is an ( m − l ) × l zero matrix. We define in-teratively Q i and R i by the QR decomposition of M i Q i − (matlab command qr( M i Q i − , Q i R i = M i Q i − , (C5)which produces a square l × l upper triangular matrix R i with eigenvalues R i,jj > j =1 , ..., l ). We store R i for each timestep. After N timesteps we have the relation Q N R N ...R = M N ...M Q . (C6)The infinite Lyapunov exponents can then be approximated by λ j = 1 N N (cid:88) i =0 ln R i,jj . (C7)For FTLEs we truncate the above sum after W timesteps, and view how this truncationchanges in time, i.e. λ j,n = 1 hW n (cid:88) i = n − W ln R i,jj . (C8)Here W = w/h , where w is the desired time window length. Note that n initialises at time W . ACKNOWLEDGMENTS
C.Q., J.S and T.L. have received funding from the European Union’s Horizon 2020research and innovation programme under the Marie Sklodowska-Curie grant agreementNo 643073. J.S. gratefully acknowledges the financial support of the EPSRC via grantsEP/N023544/1 and EP/N014391/1.
Ashwin, P. and Ditlevsen, P., “The middle pleistocene transition as a generic bifurcation on a slow manifold,”Climate Dynamics , 2683–2695 (2015).Berger, A. L., “Long-term variations of daily insolation and quaternary climatic changes,” Journal of theatmospheric sciences , 2362–2367 (1978).Chalk, T. B., Hain, M. P., Foster, G. L., Rohling, E. J., Sexton, P. F., Badger, M. P., Cherry, S. G.,Hasenfratz, A. P., Haug, G. H., Jaccard, S. L., et al. , “Causes of ice age intensification across themid-pleistocene transition,” Proceedings of the National Academy of Sciences , 13114–13119 (2017).Clark, P. U., Archer, D., Pollard, D., Blum, J. D., Rial, J. A., Brovkin, V., Mix, A. C., Pisias, N. G., andRoy, M., “The middle pleistocene transition: characteristics, mechanisms, and implications for long-termchanges in atmospheric pco 2,” Quaternary Science Reviews , 3150–3184 (2006).Crucifix, M., “Oscillators and relaxation phenomena in pleistocene climate theory,” Phil. Trans. R. Soc. A , 1140–1165 (2012).De Saedeleer, B., Crucifix, M., and Wieczorek, S., “Is the astronomical forcing a reliable and uniquepacemaker for climate? a conceptual model study,” Climate Dynamics , 273–294 (2013).Dieci, L., Russell, R. D., and Van Vleck, E. S., “On the compuation of lyapunov exponents for continuousdynamical systems,” SIAM journal on numerical analysis , 402–423 (1997).Dijkstra, H. A., Nonlinear climate dynamics (Cambridge University Press, 2013). he MPT as a result of bistability 17
Farmer, J. D., “Chaotic attractors of an infinite-dimensional dynamical system,” Physica D: NonlinearPhenomena , 366–393 (1982).Ganopolski, A. and Calov, R., “The role of orbital forcing, carbon dioxide and regolith in 100 kyr glacialcycles,” Climate of the Past , 1415–1425 (2011).Gildor, H. and Tziperman, E., “A sea ice climate switch mechanism for the 100-kyr glacial cycles,” Journalof Geophysical Research: Oceans , 9117–9133 (2001).Hays, J. D., Imbrie, J., and Shackleton, N. J., “Variations in the earth’s orbit: pacemaker of the ice ages,”Science , 1121–1132 (1976).Huybers, P., “Early pleistocene glacial cycles and the integrated summer insolation forcing,” Science ,508–511 (2006).Huybers, P., “Pleistocene glacial variability as a chaotic response to obliquity forcing,” Climate of the Past , 481–488 (2009).Huybers, P., “Combined obliquity and precession pacing of late pleistocene deglaciations,” Nature ,229–232 (2011).Huybers, P. and Eisenman, I., “Integrated summer insolation calculations. noaa/ncdc paleoclimatologyprogram data contribution ,491 (2005).Imbrie, J., Berger, A., Boyle, E., Clemens, S., Duffy, A., Howard, W., Kukla, G., Kutzbach, J., Martinson,D., McIntyre, A., et al. , “On the structure and origin of major glaciation cycles 2. the 100,000-year cycle,”Paleoceanography , 699–735 (1993).Kuznetsov, Y. A., Elements of applied bifurcation theory , Vol. 112 (Springer Science & Business Media,2013).Lisiecki, L. E. and Raymo, M. E., “A pliocene-pleistocene stack of 57 globally distributed benthic δ , 1–17 (2005).Lisiecki, L. E. and Raymo, M. E., “Plio–pleistocene climate evolution: trends and transitions in glacial cycledynamics,” Quaternary Science Reviews , 56–69 (2007).Maasch, K., “Statistical detection of the mid-pleistocene transition,” Climate dynamics , 133–143 (1988).Maasch, K. A. and Saltzman, B., “A low-order dynamical model of global climatic variability over the fullpleistocene,” Journal of Geophysical Research: Atmospheres , 1955–1963 (1990).Milankovitch, M., “History of radiation on the earth and its use for the problem of the ice ages,” K. Serb.Akad. Beogr (1941).Paillard, D., “The timing of pleistocene glaciations from a simple multiple-state climate model,” Nature , 378–381 (1998).Paillard, D., “Glacial cycles: toward a new paradigm,” Reviews of Geophysics , 325–346 (2001).Paillard, D., “Quaternary glaciations: from observations to theories,” Quaternary Science Reviews ,11–24 (2015).Paillard, D., “The plio-pleistocene climatic evolution as a consequence of orbital forcing on the carboncycle,” Climate of the Past , 1259–1267 (2017).Paillard, D. and Parrenin, F., “The antarctic ice sheet and the triggering of deglaciations,” Earth andPlanetary Science Letters , 263–271 (2004).Perko, L., Differential equations and dynamical systems , Vol. 7 (Springer Science & Business Media, 2013).Pisias, N. G. and Moore, T., “The evolution of pleistocene climate: a time series approach,” Earth andPlanetary Science Letters , 450–458 (1981).Ridgwell, A. J., Watson, A. J., and Raymo, M. E., “Is the spectral signature of the 100 kyr glacial cycleconsistent with a milankovitch origin?” Paleoceanography , 437–440 (1999).Saltzman, B. and Maasch, K. A., “Carbon cycle instability as a cause of the late pleistocene ice ageoscillations: modeling the asymmetric response,” Global biogeochemical cycles , 177–185 (1988).Saltzman, B. and Maasch, K. A., “A first-order global model of late cenozoic climatic change ii. furtheranalysis based on a simplification of co 2 dynamics,” Climate Dynamics , 201–210 (1991).Shackleton, N. J., “The 100,000-year ice-age cycle identified and found to lag temperature, carbon dioxide,and orbital eccentricity,” Science , 1897–1902 (2000).Smith, H., An introduction to delay differential equations with applications to the life sciences , Vol. 57(Springer Science & Business Media, 2010).Tzedakis, P., Crucifix, M., Mitsui, T., and Wolff, E. W., “A simple rule to determine which insolationcycles lead to interglacials,” Nature , 427–432 (2017).Tziperman, E. and Gildor, H., “On the mid-pleistocene transition to 100-kyr glacial cycles and the asym-metry between glaciation and deglaciation times,” Paleoceanography , 1–8 (2003).Yang, H. and Zhu, J., “Equilibrium thermal response timescale of global oceans,” Geophysical ResearchLetters , 1–5 (2011).Yiou, P., Ghil, M., Jouzel, J., Paillard, D., and Vautard, R., “Nonlinear variability of the climatic systemfrom singular and power spectra of late quaternary records,” Climate Dynamics9