Pleistocene Glacial Variability as a Chaotic Response to Obliquity Forcing

. The mid-Pleistocene Transition from 40 ky to ∼ 100 ky glacial cycles is generally characterized as a singular transition attributable to scouring of continental regolith or a long-term decrease in atmospheric CO 2 concentrations. Here an alternative hypothesis is suggested, that Pleistocene glacial variability is chaotic and that transitions from 40 ky to ∼ 100 ky modes of variability occur spontaneously. This alternate view is consistent with the presence of ∼ 80 ky glacial cycles during the early Pleistocene and the lack of evidence for a change in climate forcing during the mid-Pleistocene. A simple model illustrates this chaotic scenario. When forced at a 40 ky period the model chaotically transitions between small 40 ky glacial cycles and larger 80 and 120 ky cycles which, on average, give the ∼ 100 ky variability.

This hypothesis of spontaneous switching will address several features of Pleistocene glaciation.First is that the mode of glaciation switches from 40 ky to ∼100 ky timescales in the absence of long-term (i.e.> 100 ky) trends in CO 2 -as observed for the late Pleistocene (Lüthi et al., 2008) and can be inferred for the middle Pleistocene (Medina-Elizalde and Lea, 2005)-or other controls upon glaciation.Second is to illustrate how the appearance of 40 ky and ∼100 ky glacial variability over the course of the Pleistocene is consistent with a fundamental obliquity pacing of both styles of glaciation (Huybers and Wunsch, 2005;Huybers, 2007) (Fig. 1a,b).Finally, the presence of longer than 40 ky glacial cycles during the early Pleistocene is addressed (Mudelsee and Schulz, 1997;Heslop et al., 2002).The hypothesis will not, however, address the initiation of glaciation nor the longer term trends in amplitude, duration, and skew of the glacial cycles observed over the course of the Plio-Pleistocene (e.g.Huybers, 2007;Lisiecki and Raymo, 2007).

Ice volume and its rate of change
A modification to the well-known simple glacial model of Imbrie and Imbrie (1980) is motivated by examination of the relationship between ice volume and its time rate of change.An average of globally distributed benthic δ 18 O records (Huybers, 2007) is used as a proxy for ice volume (Fig. 1a), where the record is first smoothed using a 7 ky 2 P. Huybers: Chaotic obliquity response window to focus on longer term variations.(Similar results are obtained if the record of Lisiecki and Raymo (2007) is instead used.) To begin, the ice volume proxy is compared to its rate of change.Deglaciation proceeds, on average, at twice the rate of ice growth, but is much more rapid following deep glaciations.Apparently, the rate of deglaciation is sensitive to prior amounts of ice volume (Fig. 2c), agreeing with the previous recognition that glacial terminations tend to follow after a critical build-up of ice volume (Raymo, 1997).
The dependence between ice volume and its rate of change can be quantified by computing the cross-correlation at different lags (Fig. 2d).For the early Pleistocene (∼2-1 Ma), the greatest magnitude cross-correlation occurs when ice volume rates of change are lagged 9 ky behind ice volume magnitude (r = −0.71,p < 0.01, the statistical test is discussed below), but this comes as little surprise because variations are primarily at 40 ky periods (i.e.d/dt sin(2πt/40ky) = − sin(2π(t − 10ky)/40ky)).More interestingly, the crosscorrelation retains its greatest magnitude at a 9 ky lag during the late Pleistocene (r = −0.55,p < 0.01, 1-0 Ma), when variations are primarily at the longer ∼100 ky periods.This consistency in the lagged correlation structure of the δ 18 O record suggests that the physics controlling glacial variability does not fundamentally change across the mid-Pleistocene, retaining what is fundamentally a 40 ky obliquity-like relationship.If no delay is included, the observed crosscorrelation is very nearly zero and statistically insignificant.
An autoregressive model is adopted as a null-hypothesis to test the significance of the delayed cross-correlation structure between ice volume and its rate of change.In particular, the fit obtained by Wunsch (2004) is used, involving normally distributed and uncorrelated noise innovations along with a second order auto-regressive relationship having coefficients of 1.181 and -0.1984.Individual realizations are made from this auto-correlation model, each a million years long, and a record is kept of the maximum and minimum lagged crosscorrelation between the time-series and its rate of change over a range of zero to 200 lags.This process is repeated 10,000 times, and p-values are derived from the maximum correlation observed for any lag across the ensemble of realizations.As reported above, the lagged cross-correlation observed between ice volume and its rate of change during both the early and late Pleistocene are found to have a p < 0.01, indicating that the lagged relationship is highly statistically significant.Experiments with a wide range of alternate auto-correlation models yield consistent results.Apparently, greater amounts of ice volume subsequently lead to a more rapid disintegration of that ice.
The relationship between ice volume and its rate of change appears most obvious during deglaciation (Fig. 1c), and raises the question of why, during deglaciation, would the climate system remember the amount of ice volume present during the preceding glacial?Isostatic adjustment seems to occur too quickly for a 9 ky lag, but there are other possi-bilities.One is that a thicker ice sheet will be more prone to basal melting because it better isolates an ice sheet's base from the cold overlying surface and has a lower pressuremelting point.Further, the past thickness of the ice may be remembered through the the ongoing conversion of potential energy into frictional heating and the associated maintenance of melt water at the bed, which could serve to destabilize the ice sheet (e.g.Oerlemans, 1984;Marshall and Clark, 2002).Another possibility involves the carbon cycle.Perhaps the rate and amount of CO 2 fluxed from the ocean into atmosphere, or from the solid Earth into the ocean/atmosphere system via volcanism (Huybers and Langmuir, 2009), depends upon the depth of the preceding glacial.Other possibilities doubtlessly exist.

An adaptation of Imbrie and Imbrie's simple model
The statistical and potential physical significance of the delayed relationship between the magnitude of ice volume and its rate of change prompts modification of the simple glacial model proposed by Imbrie and Imbrie (1980), who posited a zero-lag relationship between ice volume and its rate of change.Their model is reformulated to make the rate of deglaciation sensitive to prior ice volume, Here V is ice volume, t is time in ky, L is the time lag, T is a time constant, a represents accumulation, and P is an exponent whose value depends on whether the model is accumulating or ablating.The model is linearly dependent on the forcing when ice is accumulating, but when ice is ablating, it sensitively depends on the ice volume L ky ago, raised to a power p.A time-variable external forcing is introduced through θ t .Negative ice volumes, V t < 0, are reset to zero, V t = 0.For convenience, ice volume is treated as being nondimensional while units of time are retained.Deglaciations tend to occur during times of increased orbital obliquity (Huybers and Wunsch, 2005;Huybers, 2007) (Fig. 1a,b), and the model forcing, θ, is parametrized to have a 40 ky period, approximating that of obliquity.A simple 40 ky sinusoid is used for the forcing to avoid obscuring the internal model dynamics with the response to the amplitude and frequency modulations associated with the actual obliquity or insolation forcing (see Berger et al., 1998;Tziperman et al., 2006).A more realistic representation of the forcing term in Eq. 1 would also include some stochastic parametrization of the myriad climatic variations not explicitly represented in the model, including weather and longerperiod variations.Trials including stochastic terms decrease the fraction of energy at the obliquity periods but, if added in moderation, do not alter the basic structure of the variability described below.
Eq. 1 is discrete and is technically a map rather than a differential equation.As is common to all simple models (e.g.Maasch and Saltzman, 1990;Ghil, 1994;Paillard, 1998;Tziperman and Gildor, 2003), Eq. 1 is best interpreted as a schematic that serves to illustrate dynamical scenarios which the climate system may be capable of.While there is no assurance that the phenomena manifested by this schematic would be reproduced in more complete models, the model's simplicity facilitates a more thorough analysis of its nonetheless rich behavior.
Model parametrization are chosen in keeping with observations during the late-Pleistocene: an ice sheet growth time scale of T = 90 ky, accumulation rate a = 0.9, lag L = 9 ky, and exponent p = 9. (These parameters are also selected for their simplicity, though many choices for the four adjustable parameters in the model yield qualitatively similar behavior.)The exponent, p = 9, is large, but note that the time rate of change in the surface height of an ice sheet is proportional to its height to the fifth power and its surface slope to the third power, assuming the standard shallow-ice approximation and a Glen's flow law with an exponent of three (e.g.van der Veen, 1999).Instabilities within the ice may be represented by even higher power-law relationships.The forcing variance is set to one, and is not considered an adjustable parameter because adjustments to a and T can produce equivalent results.

Chaotic 40 ky and ∼100 ky glacial cycles
For significant ablation to occur in the model, the obliquity forcing must be greater than a and V t−L must exceed one.This state-dependent sensitivity to external forcing is a common feature of what are broadly referred to as excitable systems (FitzHugh, 1961;Pikovsky et al., 2001).(In an early study of an excitable system, Van Der Pol (1927) focused on cardiac pacemaking by periodic electrical stimuli, in good analogy with insolation being described as the pacemaker of the ice ages (Hays et al., 1976).)As also found for some other excitable systems (Strain and Greenside, 1998;Othmer and Xie, 1999), Eq. 1 exhibits chaotic behavior.In particular, the amplitude of the model ice volume becomes chaotic under certain model parametrizations, whereas the phase is invariably synchronized with the forcing (Rosenblum and Pikovsky, 2003), in this case because some ablation always occurs near maximum obliquity.Saltzman andVerbitsky (1992, 1993) also describe the behavior of their model as chaotic, but use this term to describe mode switching in response to modulations of the orbital forcing or, in other cases, from stochastic forcing.Earth's orbital variations appear chaotic (Laskar, 1989) so that a model driven by the resulting insolation would presumably inherit this behavior.Here a distinction is made between a model being chaotic and a non-chaotic model being driven by a chaotic forcing.Following Lorenz (1963), the present model is called chaotic because, when forced by a periodic signal, it exhibits apparently deterministic, bounded, and nonperiodic solutions that are exquisitely sensitive to initial conditions.
The period of the model can be determined by counting the successive number of maxima which occur before returning to a previously visited state.If the amplitude of the forcing is made to increase starting from zero, the model exhibits a period-doubling route to chaos (Fig. 2a), in qualitative agreement with that of the logistic map and the Rossler system (e.g Strogatz, 1994).The model undergoes periodic 40 ky variations at a very small forcing amplitude, but as the amplitude increases, period doubling ensues, first giving 80 ky glacial cycles, and eventually very-long runs that do not repeatindicative of the infinite repeat times associated with chaos.At yet higher amplitudes of the forcing, periodic solutions reappear, generally with a 120 ky repeat time.Thus, the time scales inherent to Pleistocene glacial cycles-short 40 ky cycle and longer 80 ky to 120 ky cycles-are present within this simple model.
Another feature of chaotic systems is exponential divergence of model trajectories subject to small perturbations.Perturbations to V in Eq. 1 grow exponentially with a timescale of ∼300 ky.In the face of imperfect observations, the model trajectory will only be predictable for a few glacial cycles.
The model exhibits chaotic variability when it is forced by a 40 ky sinusoid with unit variance.The chaos is characterized by alternating strings of 40 ky and ∼100 ky cycles, independent of changing any model parameters.The origins of these distinct 40 ky and ∼100 ky modes of glacial variability can be understood by examining successive ice volume maxima against one another (Fig. 2b).It is useful to divide ice volume maxima into three states: interglacial (i), mild glacial (g), and full glacial (G).The model tends toward greater ice volume over most of a glacial cycle, so that if one ice volume maximum is in state i, the next is in g and then G.But the transition out of the G state, representing a deglaciation, can be into either an i or g state.These transitions between glacial states emerge as an intrinsic feature of the model, offering some explanation of the analogous three climate states that were hard-wired into the Paillard model (Paillard, 1998).
Cycling between G-g states gives 80 ky glacial cycles, and cycling between G-i-g gives 120 ky glacial cycle.80 and 120 ky cycles occur with roughly equal probability and, on average, give a ∼100 ky time-scale.(As will be noted again later, alternations between 80 and 120 ky cycles tends to generate a peak in associated spectral estimates that is centered at ∼100ky periods.)Thus, this purely obliquity paced model generates ∼100 ky variability, in keeping with observations (Huybers and Wunsch, 2005;Huybers, 2007).
The model also spontaneously generates strings of 40 ky glacial cycles.The 40 ky mode of glacial cycling results from the presence of an unstable fixed point at the boundary between the g and G states (Fig. 2b).If the model ever landed precisely on the g-G boundary it would be trapped, forever returning to the exact same state.While the probability of becoming permanently trapped is vanishingly small, the model intermittently happens to land near the unstable fixed point and then requires a number of g-G cycles to escape the fixed point's influence.This string of temporarilytrapped g-G cycles are of nearly the same amplitude, similar to the glacial cycles observed during the Pliocene and Early-Pleistocene.Notably, the model still generates an asymmetry during these 40 ky oscillations where ablation is more rapid than accumulation, in agreement with observations of the Pliocene and early Pleistocene glacial cycles (Hagelberg et al., 1991;Ashkenazy and Tziperman, 2004).

Pleistocene-like realizations
The model episodically becomes trapped near the unstable fixed point, giving trains of 40 ky glacial cycles, and then transitions into the full g-G and i-g-G cycles, giving larger amplitude and longer-period ∼100 ky glacial variability.Given a long enough run of the model, a sequence of glacial variability similar to that observed during the Plio-Pleistocene is inevitable (e.g.Fig. 2a), but it is useful to consider the probability of actually obtaining such a realization.
The requirements for obtaining a long string of 40 ky glacial cycles can be approximated from the nearly linear relationship between successive ice volume maxima in the vicinity of the unstable fixed point.The relationship can be expressed as z n+1 − z f = a(z n − z f ), where z n is the nth glacial maximum, z f is the location of the fixed point and equals 1.34, and a is the slope and equals -2.0 (see Fig. 2b).Thus, an initial glacial maximum, z n=0 , near z f will have successors that hop to alternate sides of the fixed point, moving progressively further away according to z n −z f = (z n=0 −z f )a n .This linearization around the fixed point permits exploration of the probability of obtaining long strings of 40 ky glacial cycles.
In order to be specific, 40 ky glacial cycles are defined as occurring when successive ice volume maxima, z n , have values ranging between 1.3 and 1.4.To obtain a string of forty 40 ky glacial cycles-roughly the number identifiable during the early Pleistocene-would require the magnitude of z n=0 − z f to be only 5 × 10 −13 .If z n=0 is considered as being uniformly distributed between 0.85 and 1.35, the odds of the model yielding a string of 40 ky glacial cycles commensurate with those observed during the early Pleistocene seems exceedingly small at one in 10 12 , at least given the parametrizations considered here.
But the situation may not actually be this improbable.The early Pleistocene δ 18 O stratigraphy leading up to the mid-Pleistocene Transition exhibits irregular glacial variability, with strings of three to eight 40 ky glacial interspersed with longer glacial cycles occurring near 1.8 Ma, 1.6 Ma, and 1.2 Ma (Mudelsee and Schulz, 1997;Heslop et al., 2002;Huybers, 2007) (Fig. 1a).To obtain a string of five 40 ky glacial cycles, the distance between z n=0 and z f must equal 0.02, which occurs with an approximately 4% probability.Indeed, when run over a few million years, the model tends to regularly produce strings of four to six 40 ky glacial cycles.The suggestion, then, is that during the early Pleistocene, the glacial/climate system episodically became trapped near a 40 ky fixed point and that it escaped into ∼100 ky glacial cycles during the late Pleistocene.No specific event or change in forcing is required to explain transitioning from a string of 40 ky glacial cycles into ∼100 ky glacial variability.
Of course, model parameters could be adjusted to force or nudge the model toward generating a long string of 40 ky glacial cycles followed by an interval of ∼100 ky glacial cycles.For example, changing the exponent from p = 9 to p = 6 gives purely 40 ky variations, and p = 7 or 8 gives longer glacial cycles.Changes to the forcing variance, accumulation, or time constant yields similar results.More generally, the possibility outlined here is not exclusive of changes in atmospheric CO 2 , continental regolith, or other features of the background climate changing so as to influence the course of glaciation.A weaker version of this hypothesis is that a small shift in the background conditions governing glacial variability could be sufficient to shift the model's trajectory from one of primarily 40 ky to ∼100 ky glacial cycles.
If the underlying dynamics of the glacial system are chaotic, or in a regime near chaos, the qualitative behavior of the glacial system could be exquisitely sensitive to small changes in atmospheric CO 2 , conditions at the base of an ice sheet, or other features of the climate.But the topic of forced mid-Pleistocene transitions has been focused on elsewhere (e.g.Saltzman and Sutera, 1987;Saltzman and Verbitsky, 1992;Ghil, 1994;Raymo, 1997;Mudelsee and Schulz, 1997;Paillard, 1998;Berger et al., 1999;Clark et al., 1999;Tziperman and Gildor, 2003;Raymo et al., 2006;Clark et al., 2006), and will not be further pursued here.The main point is that the specific cause of the mid-Pleistocene Transition could be intrinsic to the glacial/climate system, without a need to call upon changes in background conditions.

Comparison of a specific realization with observations
It is instructive to analyze a particular realization of the model to show how spectral energy can be concentrated alternately at 40 ky and ∼100 ky periods.A model realization similar to the δ 18 O record over the last 2 Ma is deliberately selected from a very long run of the model (Fig. 2c).Spectral analysis of the nominally early-(2 to 1.1 Ma) and late-Pleistoecene (1.1 to 0 Ma) parts of the model realization indicate that energy concentrated at the obliquity period (taken as bands between 1/41 ± 1/150 ky ) decreases from 63% to 35% of the total, while energy at ∼100 ky period (between 1/100 ± 1/150 ky) increases from 23% to 51% (Fig. 2d).
By comparison, the analogous early-and late-Pleistocene time intervals in the δ 18 O record shows a decrease from 51% to 19% at the obliquity period, and an increase from 35% to 71% at the ∼100 ky period (Fig. 1e).Thus, the redistribution of spectral energy from 40 ky to ∼100 ky periods generated spontaneously by this realization of the model is similar to the observed redistribution.Note that while individual model glacial cycles have 80 ky and 120 ky periods, these periods do not appear in the spectral estimates because concatenation of glacial cycles of alternating lengths, unlike superposition, tends to concentrate spectral energy at the average period.
Currently available methods (e.g.Gottwald and Melbourne, 2004) are not able to distinguish whether El Niño or the glacial cycles are truly chaotic, given the finite and noisy data which is available.Similarly, the simple structure relating successive glacial maxima in the model (Fig. 2b) is not identifiable in the δ 18 O data, as might be expected given the influence of temperature, noise, and the various ice sheets on the marine δ 18 O signal.
These model results share some commonality with a simple chaotic model proposed for El Niño (Tziperman et al., 1994) -both models are driven by periodic forcing, incorporate a time delay, and exhibit a period doubling route to chaos.It has been suggested that obliquity period variability in the Tropical Pacific is responsible for glacial variability (Liu and Herbert, 2004;Fedorov et al., 2006;Huybers and Molnar, 2007), but establishing such a connection requires further theoretical and observational study.

Concluding remarks
Implications are that the climate system may spontaneously transition from 40 to ∼100 ky modes of glacial variability.Changes in background climate conditions, to the extent they occur, are expected to influence the glacial cycle trajectory but, in this view, are not strictly required to explain the transition from 40 to ∼100 ky modes of glacial variability, or episodes of long-period glacial cycles prior to the transition.If the glacial cycles are chaotic, or almost so, the implied sensitivity to initial conditions-in conjunction with the inevitable stochastic influences upon glaciation-renders predictions of future glacial cycles highly uncertain.

Fig. 1 .Fig. 2 .
Fig. 1.Glacial variability over the last 2 Ma.(a) Marine δ 18 O averaged across many sediment cores (see Huybers, 2007, for more detail) and (b) variations in Earth's obliquity.Maxima in obliquity that are not accompanied by a deglaciation are shaded in yellow.Deglaciations are identified as decreases in δ 18 O exceeding one standard deviation of the record, the length of which is indicated by the bar at right.(c) δ 18 O plotted against its time rate of change (in ‰/ky).Arrows indicate the direction and average rate of change.Deglaciation is more rapid than glaciation and is most rapid after the largest positive excursions.(d) The lagged cross-correlation between δ 18 O and its rate of change for the early Pleistocene (solid line) and the late Pleistocene (dashed line).The largest magnitude cross-correlation occurs when rates of change are lagged by 9 ky for both the early (r = −0.7,p < 0.01) and late (r = −0.6,p < 0.01) Pleistocene.Note that the lagged auto-correlation relationship reaches extrema (and tends to exceed the 99% confidence interval indicated by the yellow lines) at integer multiples of the obliquity period: lags of 9, 50, 90, and 130 ky for negative excursions and 30, 70, 110, and 150 ky for positive excursions during the early Pleistocene, and lags of 9 and 130 ky as well as 70 and 150 ky for the late Pleistocene.(e) Spectrogram of the marine δ 18 O record over the last 2 Ma.The shading indicates spectral energy that has been normalized to a maximum value of one, truncated below values of 1/100, and reported in base-10 logarithmic units.