Modeling the Dynamics of Glacial Cycles

TL;DR

This paper analyzes a conceptual model of glacial cycles, demonstrating how simplified versions capture key dynamics like limit cycles and bifurcations, highlighting the role of atmospheric CO2 in glacial periodicity.

Contribution

It introduces a simplified two-dimensional model that retains essential features of the original three-dimensional system, elucidating the dynamics of glacial cycles.

Findings

The simplified model exhibits equilibrium states and limit cycles similar to the full model.

Identification of bifurcations and a Bogdanov-Takens point as organizing centers.

Symmetry breaking leads to splitting of the Bogdanov-Takens point with different local dynamics.

Abstract

This article is concerned with the dynamics of glacial cycles observed in the geological record of the Pleistocene Epoch. It focuses on a conceptual model proposed by Maasch and Saltzman [J. Geophys. Res.,95, D2 (1990), pp. 1955-1963], which is based on physical arguments and emphasizes the role of atmospheric CO2 in the generation and persistence of periodic orbits (limit cycles). The model consists of three ordinary differential equations with four parameters for the anomalies of the total global ice mass, the atmospheric CO2 concentration, and the volume of the North Atlantic Deep Water (NADW). In this article, it is shown that a simplified two-dimensional symmetric version displays many of the essential features of the full model, including equilibrium states, limit cycles, their basic bifurcations, and a Bogdanov-Takens point that serves as an organizing center for the local and global dynamics. Also, symmetry breaking splits the Bogdanov-Takens point into two, with different local dynamics in their neighborhoods.

Equations

$$\begin{array}{rl}{\tau }^{\prime }& =-\alpha I^{\prime }+\beta {\mu }^{\prime },\\ {\eta }^{\prime }& =e_I{I}^{\prime }-e_{\mu }{\mu }^{\prime }.\end{array}$$

$$\begin{array}{rl}\frac{d{I}^{\prime }}{d{t}^{\prime }}& =-s_1{\tau }^{\prime }-s_2{\mu }^{\prime }+s_3{\eta }^{\prime }-s_4{I}^{\prime },\\ \frac{d{\mu }^{\prime }}{d{t}^{\prime }}& =r_1{\tau }^{\prime }-r_2{\eta }^{\prime }-(r_3-b_3{N}^{\prime })N^{\prime }-(r_4+b_4{N}^{\prime 2}){\mu }^{\prime }-r_5{I}^{\prime },\\ \frac{d{N}^{\prime }}{d{t}^{\prime }}& =-c_0{I}^{\prime }-c_2{N}^{\prime }.\end{array}$$

$$\begin{array}{rl}\frac{d{I}^{\prime }}{d{t}^{\prime }}& =-a_0{I}^{\prime }-a_1{\mu }^{\prime },\\ \frac{d{\mu }^{\prime }}{d{t}^{\prime }}& =-b_0{I}^{\prime }+(b_1-b_4{N}^{\prime 2}){\mu }^{\prime }-(b_2-b_3{N}^{\prime })N^{\prime },\\ \frac{d{N}^{\prime }}{d{t}^{\prime }}& =-c_0{I}^{\prime }-c_2{N}^{\prime },\end{array}$$

$$a_0$$

$$b_0$$

$$t=a_0{t}^{\prime },$$

$$X=\frac{I^{\prime }}{\hat{I}},Y=\frac{{\mu }^{\prime }}{\hat{\mu }},Z=\frac{N^{\prime }}{\hat{N}},$$

$$\begin{array}{rl}\dot{X}& =-X-{\hat{a}}_{1}Y,\\ \dot{Y}& =({\hat{b}}_1-{\hat{b}}_4{Z}^2)Y-({\hat{b}}_2-{\hat{b}}_{3}Z)Z,\\ \dot{Z}& =-{\hat{c}}_{0}X-{\hat{c}}_{2}Z,\end{array}$$

$${\hat{a}}_1=\frac{a_1\hat{\mu }}{a_0\hat{I}},{\hat{b}}_1=\frac{b_1}{a_0},{\hat{b}}_2=\frac{b_2\hat{N}}{a_0\hat{\mu }},{\hat{b}}_3=\frac{b_3{\hat{N}}^2}{a_0\hat{\mu }},{\hat{b}}_4=\frac{b_4{\hat{N}}^2}{a_0},{\hat{c}}_0=\frac{c_0\hat{I}}{a_0\hat{N}},{\hat{c}}_2=\frac{c_2}{a_0}.$$

$$x=\left(({\hat{c}}_0/{\hat{c}}_2)\surd {\hat{b}}_4\right)X,y=\left({\hat{a}}_1({\hat{c}}_0/{\hat{c}}_2)\surd {\hat{b}}_4\right)Y,z=\left(\surd {\hat{b}}_4\right)Z,$$

$$\begin{array}{rl}\dot{x}& =-x-y,\\ \dot{y}& =ry-pz+s{z}^2-y{z}^2,\\ \dot{z}& =-qx-qz.\end{array}$$

$$p=\frac{{\hat{a}}_1{\hat{b}}_2{\hat{c}}_0}{{\hat{c}}_2},q={\hat{c}}_2,r={\hat{b}}_1,s=\frac{{\hat{a}}_1{\hat{b}}_3{\hat{c}}_0}{{\hat{c}}_2\surd {\hat{b}}_4}.$$

$$\overline{x}=\lim \underset{t\to \infty }{\sup }x(t).$$

$$\begin{array}{rl}\dot{x}& =-x-y,\\ \dot{y}& =ry+px-x^{2}y.\end{array}$$

$$\begin{array}{rl}\mathrm{O}& =\{(p,r)\in {\mathbb{R}}_{+}^2:r<p\text{mbox}forp<1,r<1\text{mbox}forp>1\},\\ \mathrm{I}& =\{(p,r)\in {\mathbb{R}}_{+}^2:1<r<p\text{mbox}forp>1\},\\ \mathrm{II}& =\{(p,r)\in {\mathbb{R}}_{+}^2:1<p<r\text{mbox}forr>1\},\\ \mathrm{III}& =\{(p,r)\in {\mathbb{R}}_{+}^2:p<r\text{mbox}forr<1,p<1\text{mbox}forr>1\}.\end{array}$$

$$\ddot{x}+g(x)\dot{x}+f(x)=0,$$

$$f(x)=x^3-(r-p)x,g(x)=x^2-(r-1).$$

$${\ell }^{\ast }=-\frac{{\omega }^{\ast }}{8}{\frac{d}{dx}\frac{g^{\prime }(x)}{f^{\prime }(x)}|}_{x=x^{\ast }},$$

$$\begin{array}{rl}\dot{x}& =y,\\ \dot{y}& =(r-p)x+(r-1)y-x^{2}y-x^3,\end{array}$$

$$\left(\begin{array}{cc}0& 1\\ r-p-3(x^{\ast }{)}^2& r-1-(x^{\ast }{)}^2\end{array}\right).$$

$$x(t)=\eta u(\tilde{t}),y(t)={\eta }^{2}v(\tilde{t}),\tilde{t}=\eta t.$$

$$\begin{array}{rl}\dot{u}& =v,\\ \dot{v}& =\mu u-u^3+\eta (\lambda -u^2)v.\end{array}$$

$$\lambda =\frac{r-1}{{\eta }^2},\mu =\frac{r-p}{{\eta }^2}.$$

$$(\lambda -\mu )(r-1)=\lambda (p-1).$$

$$\begin{array}{rl}\dot{u}& =v,\\ \dot{v}& =\mu u-u^3.\end{array}$$

$$\begin{array}{rl}\dot{u}& =v,\\ \dot{v}& =f(u)+\eta g(u,v).\end{array}$$

$$\begin{array}{rl}\dot{u}& =v,\\ \dot{v}& =f(u).\end{array}$$

$$M(\lambda ,x)={\int }_{I(x)}(\lambda -u^2(t))v^2(t)dt={\oint }_{\gamma (x)}(\lambda -u^2)v(u)du=\lambda I_0(x)-I_2(x),$$

$$I_0(x)={\oint }_{\gamma (x)}v(u)du,I_2(x)={\oint }_{\gamma (x)}{u}^{2}v(u)du.$$

$$\lambda =R(x),\text{mbox}whereR(x)=\frac{I_2(x)}{I_0(x)}.$$

$$\gamma(\surd 2)=\{(\surd{2}\,\mbox{$\mathrm{sech}$}\,t,-\surd{2}\,\mbox{$\mathrm{sech}$}\,t\tanh\,t),\;t\in\mathbb{R}\},$$

Full text

1 Introduction

Earth’s climate during the Pleistocene Epoch—the geological period from approximately 2.6 million years before present (2.6 Myr BP) until approximately 11.7 thousand years before present (11.7 Kyr BP)—is of great interest in the geosciences community. The geological record gives evidence of cycles of advancing and retreating continental glaciers and ice sheets, mostly at high latitudes and high altitudes and especially in the Northern Hemisphere.

To reconstruct the Pleistocene climate, geoscientists rely on geological proxies, particularly a dimensionless quantity denoted by $\delta^{18}\mathrm{O}$, which is measured in parts per mille. This quantity measures the deviation of the ratio ${}^{18}\mathrm{O}/{}^{16}\mathrm{O}$ of the stable oxygen isotopes ${}^{18}\mathrm{O}$ and ${}^{16}\mathrm{O}$ in a given marine sediment sample from the same ratio in a universally accepted standard sample. The relative amount of the isotope ${}^{18}\mathrm{O}$ in ocean water is known to be higher at tropical latitudes than near the poles, since water with the heavier oxygen isotope is slightly less likely to evaporate and more likely to precipitate first. Similarly, water with the lighter isotope ${}^{16}\mathrm{O}$ is more likely to be found in ice sheets and in rain water at high latitudes, since it is favored in atmospheric transport across latitudes. The global distribution of $\delta^{18}\mathrm{O}$ in ocean water therefore varies in a known way between glacial and interglacial periods. A record of these variations is preserved in the calcium carbonate shells of foraminifera, a class of common single cell marine organisms. These fossil records may be sampled from deep sea sediment cores, and their age and $\delta^{18}\mathrm{O}$ may be determined. Details are described in lisiecki2005pliocene.

Figure 1.1 shows the LR04 time series of $\delta^{18}\mathrm{O}$ over the past 5.3 million years, reconstructed from sediment core data collected at 57 geographically distributed sites around the globe [31]. As the observed isotope variations are similar in shape to the temperature variations reconstructed from ice core data for the past 420 Kyr at the Vostok Station in Antarctica, the values of $\delta^{18}\mathrm{O}$ (right scale) have been aligned with the reported temperature variations from the Vostok ice core (left scale) [40]. The graph shows a relatively stable temperature during the period preceding the Pleistocene and increasing variability during the Pleistocene.

The typical pattern throughout most of the Pleistocene resembles that of a sawtooth wave, where a slow glaciation is followed by a rapid deglaciation. In the early Pleistocene, until approximately 1.2 Myr BP, the period of a glacial cycle averages 41 Kyr; after the mid-Pleistocene transition, which occurred from approximately 1.2 Myr BP until approximately 0.8 Myr BP, the glacial cycles of the late Pleistocene have a noticeably greater amplitude, and their period averages 100 Kyr. Figure 1.2 shows the global mean temperature for the past 420 Kyr, reconstructed from Vostok ice core data. The 100 Kyr cycle and the sawtooth pattern are clearly visible.

These observations suggest a number of questions for climate science. What caused the glacial oscillations during the Pleistocene? Why were the periods of the glacial cycles during the early and late Pleistocene different? What could possibly have caused the transition from 41 Kyr cycles to 100 Kyr cycles during the mid-Pleistocene?

In this article, we discuss a conceptual model of the Pleistocene climate proposed by Maasch and Saltzman in [32] to explain the phenomenon of glacial cycles. The model is conceptual, in the sense that it describes the state of the climate in a few variables, ignoring most of the processes that go into a complete climate model, but still captures the essence of the phenomenon. It is based on sound physical principles and, as we will see, makes for an interesting application of dynamical systems theory.

The numerical continuation results for the bifurcation curves reported in this article were obtained using the software package AUTO [15]; see also [14, 16]. Some recent texts on issues of climate dynamics are Cook2013, Dijkstra2013, OT131, MarshallPlumb2007, SunBryan2013.

In Section 2, we present background information to motivate the particular choices underlying the Maasch–Saltzman model. In Section 3, we derive the Maasch–Saltzman model from physical principles and formulate it as a dynamical system in a three-dimensional state space with four parameters. In Section 4, we introduce two simplifications that render the Maasch–Saltzman model symmetric and reduce it to a two-dimensional dynamical system with two parameters that can be analyzed rigorously and completely. In Section 5, we introduce asymmetry into the simplified two-dimensional model and show the effects of symmetry breaking. In the final Section 6, we summarize our results.

2 Background

There is general agreement that the periodicity of the glacial cycles is related to variations in the Earth’s orbital parameters [30]. To first order, Earth’s climate is driven by the Sun. The Earth receives energy from the Sun in the form of ultraviolet (short wavelength) radiation. This energy is redistributed around the globe and eventually reemitted into space in the form of infrared (long wavelength) outgoing radiation. The amount of energy reaching the top of the atmosphere per unit area and per unit time is known as the insolation (incident solar radiation), which is measured in watts per square meter.

2.1 Orbital Forcing

The insolation varies with the distance from the Earth to the Sun and thus depends on Earth’s orbit around the Sun. This is the basis of the Milankovitch theory of orbital forcing [26, 35].

The Earth moves around the Sun in an elliptical orbit; its eccentricity varies with time but has dominant frequencies at approximately 100 Kyr and 400 Kyr. As the Earth moves around the Sun, it rotates around its axis. The axis is tilted with respect to the normal to the orbital plane; the tilt, known as obliquity, is also close to periodic, with a dominant frequency of approximately 41 Kyr. (This tilt is the main cause of the seasonal variation of our climate.)

In addition, the Earth is like a spinning top wobbling around its axis of rotation. This component of the Earth’s orbit is called precession; its period varies from 19 to 23 Kyr.

Given the three orbital parameters (eccentricity, obliquity, precession), one can compute the insolation at any latitude and at any time of the year. An example is given in Figure 2.1,

which shows the time series of $Q^{65}$—the average insolation at $65^{\mathrm{o}}$ North—during the month of July for the past 4.5 million years; other months show a similar behavior. A cycle with a period of approximately 400 Kyr is clearly visible. A spectral analysis reveals a dominant frequency around 21 Kyr coming from two clustered spikes in the power spectrum and another, smaller frequency component at approximately 41 Kyr.

2.2 Atmospheric Carbon Dioxide

The Pleistocene climate and, in particular, the mid-Pleistocene transition are topics of great interest in the geosciences community. The 41 Kyr glacial cycles of the early Pleistocene are commonly attributed to the 41 Kyr cycle of Earth’s obliquity; see, for example, [22, 43]. In contrast, there is less agreement on the origin of the 100 Kyr cycles of the late Pleistocene.

Some authors [18, 20, 25] attribute the 100 Kyr cycles to the eccentricity of Earth’s orbit. However, simple energy balance considerations imply that variations in eccentricity are too weak to explain the surface temperature variations that are observed in the paleoclimate record.

A possible way for orbital effects to influence the Earth’s surface temperature is suggested by the greenhouse effect and the almost perfect correlation between fluctuations in the atmospheric $\mathrm{CO}_{2}$ concentration and the surface temperature that is observed, for example, in the Vostok ice core data [40]; see Figure 2.2.

Orbital variations have only a very weak effect on the composition of a planet’s atmosphere. Moreover, they would affect not only $\mathrm{CO}_{2}$, but other atmospheric components as well, and their effect would be negligible on the time scale of the glacial cycles. Hence, it is unlikely that orbital variations alone are the direct cause of the fluctuations in atmospheric $\mathrm{CO}_{2}$ concentrations, so one would need some feedback mechanism that reinforces the influence of orbital forcing on surface temperature through changes in atmospheric $\mathrm{CO}_{2}$ [44]. Changes in the carbon cycle and the climate system would then amplify each other to produce the glacial cycles, and atmospheric $\mathrm{CO}_{2}$ would have to play a central role in such a feedback mechanism.

Saltzman was the first to propose a conceptual climate model that highlighted the role of atmospheric $\mathrm{CO}_{2}$ in the dynamics of glacial cycles [45]. The model was further developed in joint work with Maasch in a series of articles [32, 46, 47]. In this article, we focus on the model proposed by Maasch and Saltzman in [32].

2.3 Other Models

There certainly is no unique way to explain the phenomenon of glacial cycles and the Pleistocene climate in a comprehensive manner. Other conceptual models can be found, for example, in [3, 4, 9, 18, 24, 36, 37, 38, 39, 42, 48, 51]. An interesting case was made by Huybers [24], who argued that the reconstruction of the temperature record from proxy data presented in Figure 1.1 relies on orbital assumptions and is therefore subject to bias. Huybers developed an unbiased age model which does not rely on orbital assumptions and showed that the late Pleistocene glacial terminations are paced by changes in Earth’s obliquity [23]. This theory would imply that the entire Pleistocene climate regime can be explained by obliquity alone. We refer the reader to [11] for a summary of the current state of the art.

As a final note, we caution that any model is a mathematical construct, and any phenomenon that results from its analysis is merely a manifestation of the assumptions underlying the model. The question whether the model reflects the true cause(s) of the glacial cycles lies outside the domain of mathematics.

3 The Maasch and Saltzman Model

The Pleistocene climate model proposed by Maasch and Saltzman [32] involves five state variables. They are, with their associated units in square brackets:

• •

the total global ice mass, $I$ [kg],

• •

the atmospheric $\mathrm{CO}_{2}$ concentration, $\mu$ [ppm],

• •

the volume of the North Atlantic Deep Water (NADW), $N$ [m3],

• •

the global mean sea surface temperature (SST), $\tau$ [K], and

• •

the mean volume of permanent (summer) sea ice, $\eta$ [m3].

The volume of the NADW is a measure of the strength of the global oceanic circulation (thermohaline circulation, THC). We can think of $N$ as a measure of the strength of the oceanic $\mathrm{CO}_{2}$ pump, since the oceanic $\mathrm{CO}_{2}$ pump is an integral part of the THC. The other variables are self-explanatory.

The state variables vary with time, albeit on rather different time scales. The total global ice mass, atmospheric $\mathrm{CO}_{2}$ concentration, and NADW vary on the order of thousands of years, while the SST and summer sea ice vary on the order of decades or centuries. Here the focus is on the slow time scale, where we assume that the fast variables equilibrate essentially instantaneously. That is, the long-term dynamics of the climate system are described in terms of $I$, $\mu$, and $N$ (the prognostic variables); $\tau$ and $\eta$ are diagnostic variables, which follow the prognostic variables in time.

3.1 Model Formulation

The climate model is formulated in terms of anomalies—deviations from long-term averages, which are indicated with a prime. The governing equations follow from plausible physical principles, which are detailed in [46, §2].

The global mean SST ($\tau$) and the mean volume of permanent sea ice ($\eta$) vary with the total global ice mass ($I$) and the atmospheric $\mathrm{CO}_{2}$ concentration ($\mu$) but are independent of the NADW ($N$); in particular, $\tau$ decreases as $I$ increases or $\mu$ decreases, while $\eta$ increases as $I$ increases or $\mu$ decreases. To leading order, the dependences are linear, so

[TABLE]

In the absence of external forces, the governing equations for $I^{\prime}$, $\mu^{\prime}$, and $N^{\prime}$ are

[TABLE]

Time $t^{\prime}$ is measured in units of one year [yr]. The coefficients in these equations are positive (or zero). Maasch and Saltzman also included an external forcing term related to $Q^{65}$ in the equation for $I^{\prime}$, but since we are interested in the internal dynamics of the system, we don’t include external forcing.

The physical assumptions underlying these equations are:

$\bullet$ The prognostic variables relax to their respective long-term averages, so their anomalies tend to zero as time increases; in particular, $I^{\prime}$ and $N^{\prime}$ decay at a constant rate, while the decay rate of $\mu^{\prime}$ increases quadratically with $N^{\prime}$  [46, §2.V].

$\bullet$ If the SST exceeds its mean value ($\tau^{\prime}>0$), the total global ice mass decreases and the atmospheric $\mathrm{CO}_{2}$ concentration increases (due to outgassing); if the SST is less than its mean value ($\tau^{\prime}<0$), the opposite happens. The coupling is linear to leading order.

$\bullet$ Since $\mathrm{CO}_{2}$ is a greenhouse gas, an increase in the atmospheric $\mathrm{CO}_{2}$ concentration leads to a warmer climate and thus a decrease in the total global ice mass.

$\bullet$ If the volume of permanent sea ice exceeds its mean value ($\eta^{\prime}>0$), the total global ice mass increases and the atmospheric $\mathrm{CO}_{2}$ concentration decreases; if the volume of permanent sea ice is less than its mean value ($\eta^{\prime}<0$), the opposite effect happens. The coupling is linear to leading order.

$\bullet$ A greater-than-average total global ice mass ($I^{\prime}>0$) negatively affects both the atmospheric $\mathrm{CO}_{2}$ concentration and the strength of the North Atlantic overturning circulation; a less-than-average total global ice mass ($I^{\prime}<0$) has the opposite effect. The coupling is linear to leading order.

$\bullet$ The atmospheric $\mathrm{CO}_{2}$ concentration decreases as the strength of the North Atlantic overturning circulation increases, but the coupling weakens (strengthens) as the strength of the NADW is above (below) average [46, §2.III a,b].

Upon substitution of the expressions (3.1), the governing equations (3.2) become

[TABLE]

where

[TABLE]

Following [32, 46], we take $b_{0}=0$ and assume that the remaining coefficients are all positive. Note that $a_{0}$ and $b_{1}$ involve positive as well as negative contributions, so the implicit assumption is that the positive contributions dominate.

3.2 Nondimensional Model

Next, we reformulate the system of equations (3.3) by rescaling time,

[TABLE]

and introduce dimensionless variables,

[TABLE]

where $\hat{I}$, $\hat{\mu}$, and $\hat{N}$ are reference values of $I$, $\mu$, and $N$, respectively. Since $a_{0}\approx 1.00\cdot 10^{-4}\,\mathrm{yr}^{-1}$, a unit of $t$ corresponds to (approximately) 10 Kyr.

The governing equations for $X$, $Y$, and $Z$ are

[TABLE]

where the dot $\dot{}$ indicates differentiation with respect to $t$. Recall that we have set $b_{0}=0$. The remaining coefficients are dimensionless combinations of the physical parameters in the system of equations (3.3),

[TABLE]

A rescaling of the variables $X$, $Y$, and $Z$,

[TABLE]

leads to the following dynamical system for the triple $(x,y,z)$:

[TABLE]

The coefficients $p$, $q$, $r$, and $s$ are combinations of the physical parameters,

[TABLE]

The coefficients are assumed to be positive, with $q>1$. The system of equations (3.8) is the model proposed by Maasch and Saltzman in [32, Eqs. (4)–(6)]. Note that the model considered here does not include external forcing, so it describes the internal dynamics of the climate system.

3.3 Discussion

The system of equations (3.8) is what is known as a conceptual model. Its derivation involves physical arguments, but there is no guarantee that it corresponds to what actually happened in the climate system during the Pleistocene. Its sole purpose is to describe a possible mechanism that explains the observed behavior of the glacial cycles.

Loosely speaking, we identify $x$, $y$, and $z$ with the anomalies of the total amount of ice, the atmospheric $\mathrm{CO}_{2}$ concentration, and the volume of the NADW (the strength of the oceanic $\mathrm{CO}_{2}$ pump), respectively. Time is normalized and expressed in units of the characteristic time of the total global ice mass, typically of the order of 10 Kyr.

Because of the various transformations needed to get from the physical system (3.3) to the dynamical system (3.8), it is difficult to relate the parameters to actual physical processes. The best we can do is look at their effect on the possible solutions. For example, a nonzero value of the parameter $s$ renders the problem asymmetric, so $s$ is introduced to achieve the observed asymmetry of the glacial cycles. The coefficient $q$ is the characteristic time of NADW (expressed in units of the characteristic time of the total global ice mass). The assumption $q>1$ implies that NADW changes on a faster time scale than the total global ice mass, and as $q$ increases, this change occurs on an increasingly faster time scale. If we rewrite the second equation as $\dot{y}=(r-z^{2})y-pz-sz^{2}$, we see that the growth rate $r$ of the atmospheric $\mathrm{CO}_{2}$ concentration is balanced by the anomaly of NADW. Lastly, the coefficient $p$ expresses the sensitivity of the atmospheric $\mathrm{CO}_{2}$ concentration to NADW.

Conceptually, the following sequence of events hints at the possible existence of periodic solutions: (i) As the amount of $\mathrm{CO}_{2}$ in the atmosphere increases and $y$ becomes positive, the total amount of ice decreases and $x$ becomes negative (first equation); (ii) As $x$ becomes negative, the volume of NADW increases and $z$ becomes positive (third equation); (iii) As $z$ becomes positive, the amount of atmospheric $\mathrm{CO}_{2}$ decreases and $y$ becomes negative (second equation). This is the first part of the cycle. Once $y$ is negative, the opposite effects happen. (iv) The total global ice mass starts to increase again and $x$ becomes positive (first equation); (v) As a result, the volume of NADW decreases and eventually $z$ becomes negative (third equation); (vi) Once $z$ is negative, $y$ starts to increase again (second equation). This completes the full cycle and sets the stage for the next cycle. Of course, these arguments do not guarantee the existence of a periodic cycle and do not say anything about its period. The particulars will depend critically on the parameter values.

3.4 Computational Results

Maasch and Saltzman found computationally that the system (3.8) generates a limit cycle with a 100 Kyr period at the parameter values $p=1.0$, $q=1.2$, $r=0.8$, and $s=0.8$. The limit cycle is shown in Figure 3.1.

The three curves represent the total ice mass (black), the atmospheric $\mathrm{CO}_{2}$ concentration (red), and the volume of NADW (blue) in arbitrary units. Each cycle is clearly asymmetric: a rapid deglaciation is followed by a slow glaciation. Also, the three variables are properly correlated: as the concentration of atmospheric $\mathrm{CO}_{2}$  (a greenhouse gas) increases, the climate gets warmer and the total ice mass decreases; as the volume of NADW increases, the strength of the North Atlantic overturning circulation increases, more atmospheric $\mathrm{CO}_{2}$ is absorbed by the ocean and, consequently, the atmospheric $\mathrm{CO}_{2}$ concentration decreases.

More detailed numerical calculations show that the Maasch-Saltzman model possesses limit cycles in large portions of parameter space. We integrated the system (3.8) forward in time for a range of values of $(p,r)$, keeping $q$ and $s$ fixed at the values $q=1.2$ and $s=0.8$, starting from a randomly chosen initial point for each $(p,r)$, until there was a clear indication of either a limit cycle or a limit point. We then determined the quantity $\overline{x}$ for each pair $(p,r)$,

[TABLE]

Figure 3.2 shows the function $(p,r)\mapsto\overline{x}(p,r)$ as a color map. A limiting value 0 (light green) indicates convergence to the trivial state, a nonzero negative value (dark green) convergence to a nontrivial equilibrium state, and a nonzero positive value (orange or pink) convergence to a limit cycle with a finite amplitude. Limit cycles were observed in the entire orange-colored region of the $(p,r)$ plane.

These findings, as well as other results reported by Maasch and Saltzman in [32], especially when the effects of orbital forcing are included, suggest that the conceptual model (3.8) may indeed provide an explanation for the Pleistocene climate record.

4 Simplifying the Maasch–Saltzman Model

The system (3.8) has four positive parameters, $p$, $q$, $r$, and $s$, where $q>1$. As noted in Section 3.3, the parameter $s$ introduces asymmetry into the model. If $s=0$, the equations are invariant under reflection: if $(x,y,z)$ is a solution, then so is $(-x,-y,-z)$. Note furthermore that, as $q\to\infty$, the differential equation for $z$ reduces, at least formally, to the identity $z=-x$, so the system (3.8) becomes two-dimensional. These observations suggest that it may be helpful to analyze the dynamics of the Maasch–Saltzman model (3.8) in stages, where we first focus on the special case $q=\infty$ and $s=0$, and then consider the effects of finite values of $q$ and positive values of $s$.

If we set $q=\infty$ and $s=0$, the system (3.8) reduces formally to a two-dimensional system with $\mathbb{Z}_{2}$ symmetry,

[TABLE]

The state variables are $x$ and $y$, the state space is $\mathbb{R}^{2}$, $p$ and $r$ are parameters, and the parameter space is $\mathbb{R}_{+}^{2}$. Its dynamics can be analyzed rigorously and completely.

4.1 Equilibrium States and Their Stability

The origin $P_{0}=(0,0)$ is an equilibrium state of the system (4.1) for all values of $p$ and $r$. If $r>p$, there are two additional equilibrium states, $P_{1}=(x_{1}^{*},-x_{1}^{*})$ with $x_{1}^{*}=\sqrt{r-p}$, and $P_{2}=(x_{2}^{*},-x_{2}^{*})$ with $x_{2}^{*}=-\sqrt{r-p}$. A (local) linear stability analysis shows that $P_{0}$ is stable if $0<r<\min(p,1)$, unstable otherwise; $P_{1}$ and $P_{2}$ are stable if $0<p<\min(r,1)$, unstable otherwise. Thus, the parameter space is partitioned into four regions,

[TABLE]

The regions are shown in Figure 4.1, together with representative trajectories in regions O, I, and II.

The diagonal $r=p$ is the locus of pitchfork bifurcations, where $P_{1}$ and $P_{2}$ are created as stable nodes as $(p,r)$ crosses the diagonal from region O into region III or as unstable nodes as $(p,r)$ crosses the diagonal from region I into region II. The parabolic curves $C_{1}=\{p=\tfrac{1}{4}(r+1)^{2}\}$ and $C_{2}=\{r=\tfrac{1}{4}(p+1)^{2}\}$ (dashed curves shown in purple), which are tangent to the diagonal $r=p$ at the point $(1,1)$, mark the boundaries between spirals and nodes.

4.2 Hopf Bifurcation

The system (4.1) is equivalent with the Liénard equation

[TABLE]

where $f$ and $g$ are polynomial functions,

[TABLE]

A Hopf bifurcation occurs when $f(x^{*})=0$, $f^{\prime}(x^{*})>0$, and $x^{*}$ is a simple zero of $g$. The natural frequency of oscillations is $\omega^{*}=\sqrt{f^{\prime}(x^{*})}$, and the first Lyapunov coefficient at $x^{*}$ is

[TABLE]

see [28, §3.5]. If $\ell^{*}<0$, the Hopf bifurcation is supercritical; if $\ell^{*}>0$, it is subcritical. The sign of $\ell^{*}$ is the same as that of $f^{\prime\prime}(x^{*})g^{\prime}(x^{*})-f^{\prime}(x^{*})g^{\prime\prime}(x^{*})$, which in the case of the polynomial functions $f$ and $g$ given in (4.4) is the same as that of $3(x^{*})^{2}+r-p$.

The red line segments in Figure 4.1 are loci of Hopf bifurcations. On the horizontal segment at $r=1$, which is associated with $P_{0}$, we have $x^{*}=0$. The sign of $\ell^{*}$ is the same as that of $1-p$, which is negative, so the Hopf bifurcation is supercritical. The natural frequency is $\omega^{*}=\sqrt{p-1}$. On the vertical segment at $p=1$, which is associated with $P_{1}$ and $P_{2}$, we have $(x^{*})^{2}=r-p$. The sign of $\ell^{*}$ is the same as that of $4(r-1)$, which is positive, so the Hopf bifurcation is subcritical. The natural frequency is $\omega^{*}=\sqrt{2(r-1)}$.

4.3 Organizing Center

The point $(p,r)=(1,1)$ plays a pivotal role in understanding the complete dynamics of the system (4.1). To see why, rotate the coordinate system by the transformation $(x,-(x+y))\mapsto(x,y)$. In the new coordinates, the system (4.1) is

[TABLE]

and the equilibrium points are $P_{0}=(0,0)$, $P_{1}=(x_{1}^{*},0)$, and $P_{2}=(x_{2}^{*},0)$.

Let $P=(x^{*},0)$ be any of the equilibrium points, with $x^{*}=0$, $x_{1}^{*}$, or $x_{2}^{*}$. The Jacobian of the vector field at $P$ is

[TABLE]

The matrix has a double-zero eigenvalue at the point $(p,r)=(1,1)$, so the system (4.6) undergoes a Bogdanov–Takens (BT) bifurcation. Holmes and Rand [21] refer to such a point as an organizing center. Specifically, given the $\mathbb{Z}_{2}$ symmetry of the system (4.6), the point $(p,r)=(1,1)$ is a BT point with $\mathbb{Z}_{2}$ symmetry; examples of such points are discussed in [7, Ch. 4.2] and [27, § 8.4].

The system (4.6) can be analyzed in the neighborhood of the organizing center by the unfolding procedure outlined in the original papers by Bogdanov [5] and Takens [50, pp. 23–30] (reprinted in [6, Chapter 1]) and described in the textbooks of Guckenheimer and Holmes [19, §7.3] and Kuznetsov [28, §8.4]. Here, we summarize the results; the details are given in Section 4.5 below.

Near the point $(p,r)=(1,1)$, region III of Figure 4.1 decomposes into three subregions; see Figure 4.2.

In region IIIa, there is one stable limit cycle, with a pair of unstable limit cycles in its interior, one around each of the equilibrium states $P_{1}$ and $P_{2}$. As $(p,r)$ transits from region IIIa into region IIIb, the two unstable periodic solutions merge to become a pair of unstable homoclinic orbits to the saddle $P_{0}$. This homoclinic bifurcation curve (shown in blue) is tangent to the line $r-1=-4(p-1)$ at the organizing center $(1,1)$. In region IIIb, there is one stable limit cycle with an unstable limit cycle in its interior. As $(p,r)$ transits from region IIIb into region IIIc, there is a curve of saddle-node bifurcations of limit cycles, along which the stable and unstable limit cycles disappear. This curve (shown in black) is tangent to the line $r-1\approx-3.03(p-1)$ at the organizing center $(1,1)$. In region IIIc, only the three equilibrium states remain, $P_{0}$ as an unstable saddle, $P_{1}$ and $P_{2}$ as stable spirals or nodes.

4.4 Computational Results

To complement the analysis, we performed an integration of the system (4.1) forward in time for a range of values of $(p,r)$, following the procedure described in Section 3.4 for Figure 3.2, to determine the quantity $\overline{x}=\lim\sup_{t\to\infty}x(t)$, as in (3.10). Figure 4.3 shows the function $(p,r)\mapsto\overline{x}(p,r)$ as a color map, together with the bifurcation curves obtained with AUTO.

A limiting value 0 (light green) indicates convergence to the trivial state, a nonzero negative value (dark green) convergence to a nontrivial equilibrium state, and a nonzero positive value (orange or pink) convergence to a limit cycle with a finite amplitude.

The stability region O of the trivial state $P_{0}$ is clearly visible, as is its Hopf bifurcation curve. As one crosses the Hopf bifurcation curve in the direction of increasing $r$, the color changes slowly toward increasing values of $\overline{x}$, indicating a supercritical Hopf bifurcation and periodic orbits with amplitudes $\mathcal{O}(\sqrt{r-1})$.

Throughout regions I and II, the color map changes from orange to pink as $r$ increases, indicating that the solutions all approach the stable limit cycle that surrounds $P_{0}$; cf. Figure 4.1.

In region IIIa, there is a similar shift to pink as $r$ increases. One also sees some green and orange patches in region IIIa, indicating that some of the randomly chosen initial conditions lie in the basins of attraction of the stable equilibrium states $P_{1}$ and $P_{2}$. Next, in region IIIb, the color map has largely the same characteristics as in region IIIa, corresponding to the fact that solutions with initial conditions that lie inside the large unstable limit cycle approach one of the stable equilibria (green or orange), and those with initial conditions outside the unstable limit cycle approach the large stable limit cycle (pink). Finally, in region IIIc, the color map consists entirely of green and orange, indicating that all of the solutions are attracted either to $P_{1}$ or to $P_{2}$, as expected since there are no stable limit cycles in region IIIc.

4.5 Bogdanov–Takens Unfolding

In Section 4.3, we presented the results of a bifurcation analysis of the system (4.6) in a neighborhood of the organizing center at $(p,r)=(1,1)$. In this section, we present the details of the Bogdanov–Takens unfolding procedure used to establish these results. The section is somewhat technical, but since it is self-contained, it can be skipped at first reading.

The unfolding is achieved by introducing a small positive parameter $\eta$ and rescaling the dependent and independent variables,

[TABLE]

If $(x,y)$ is a solution of the system (4.6), then $(u,v)$ must satisfy the system

[TABLE]

Here, the dot $\dot{\ }$ stands for differentiation with respect to the variable $\tilde{t}$, and $\lambda$ and $\mu$ are parameters, which are defined in terms of $p$ and $r$,

[TABLE]

Note that $\mu$ is negative in region I and positive in regions II and III (Figure 4.1). Henceforth, we omit the tilde and write $t$, instead of $\tilde{t}$.

Remark. The definition (4.10) of $\lambda$ and $\mu$ generates a linear relation between $p$ and $r$,

[TABLE]

This is the equation of a pencil through the organizing center $(1,1)$ parameterized by $\lambda$. Referring to the regions labeled I, II, and III in Figure 4.1, we note that $\lambda$ increases from 0 to infinity as one rotates counterclockwise from the horizontal line $\{p>1,r=1\}$ through region I, then decreases as one continues to rotate through region II, until $\lambda=1$ at the vertical line $\{p=1,r>1\}$, and decreases further as one rotates through region III, until $\lambda=0$ at the horizontal line $\{0<p<1,r=1\}$.

The results of the local analysis of Sections 4.1 and 4.2 may be recovered directly from the system (4.9), as follows. The origin $(0,0)$ is an equilibrium state of (4.9) for all $\lambda$ and $\mu$, and if $\mu>0$, there are two additional equilibrium states, $(\pm\surd\mu,0)$. A linearization of (4.9) with $\mu<0$ about $(0,0)$ shows that the real parts of the two eigenvalues pass through zero at $\lambda=0$, which corresponds to the line of supercritical Hopf bifurcations $\{p>1,r=1\}$. Similarly, a linearization of (4.9) with $\mu>0$ about $(\pm\surd\mu,0)$ shows that the real parts of the two eigenvalues pass through zero at $\lambda=\mu$, which corresponds to the line of subcritical Hopf bifurcations $\{p=1,r>1\}$.

4.5.1 Hamiltonian Structures

The equilibrium states of the system (4.9) are independent of $\eta$, so they persist as $\eta\to 0$. In the limit, (4.9) reduces to the Hamiltonian system

[TABLE]

Closed orbits of this system are level curves of the Hamiltonian, $H(u,v)=\tfrac{1}{2}v^{2}-\tfrac{1}{2}\mu u^{2}+\tfrac{1}{4}u^{4}$. If $\mu<0$, $H(u,v)$ reaches its minimum value 0 at the origin, so the closed orbits are nested and surround the origin. The more interesting case is $\mu>0$, where $H(u,v)$ reaches its minimum value at $(\pm\surd\mu,0)$, and $H(u,v)=0$ at the saddle point at the origin. We will analyze the case $\mu>0$ in detail and return to the case $\mu<0$ in Section 4.5.5. Recall that $\mu>0$ implies that $r>p$, so the following analysis applies to the regions II and III in Figure 4.1.

Figure 4.4 shows the phase portrait of the Hamiltonian system (4.12) with $\mu=1$. (The phase portrait for other positive values of $\mu$ is similar.)

We see that there are several types of closed orbits. There is a pair of homoclinic orbits to the origin, there are periodic orbits in the interior of each of the homoclinic orbits surrounding the equilibrium states $(\pm 1,0)$, and there are large-amplitude periodic orbits external to the homoclinic orbits. The question is, which of these closed orbits persist as the Hamiltonian system (4.12) is perturbed to the system (4.9). Because of the $\mathbb{Z}_{2}$ symmetry, it suffices to consider the homoclinic orbit in the right half of the $(u,v)$ plane and the periodic orbits in its interior; the results for the closed orbits in the left half of the $(u,v)$ plane follow by reflection. Of course, we also need to consider the large-amplitude periodic orbits that are external to the homoclinic orbits. Notice the clockwise orientation of all these orbits.

4.5.2 Melnikov Function

The persistence of closed (homoclinic or periodic) orbits of Hamiltonian systems under perturbations can be analyzed by means of the Melnikov function. This function dates back at least to Poincaré [41]; it features in articles by Melnikov [34] and Arnold [2] and in the book by Andronov et al. [1]. A definitive discussion can be found in the book of Guckenheimer and Holmes [19, §4.5]. The Melnikov function and the associated theory apply to a range of different systems, but for our purposes it suffices to summarize the results for the general system

[TABLE]

(The functions $f$ and $g$ are not to be confused with those in Section 4.2.) In the limit as $\eta\to 0$, this system reduces to the Hamiltonian system

[TABLE]

Let $\Gamma_{0}=\{t\mapsto(u(t),v(t)),\,t\in I\}$ be any closed orbit of (4.14). The Melnikov function associated with $\Gamma_{0}$ is the integral $\int_{I}g(u(t),v(t))v(t)\,dt$. Thus, the Melnikov function measures the cumulative effect of the projection of the perturbed component of the vector field, $[0\;g]^{\mathrm{t}}$, on the normal vector, $[-f\;v]^{\mathrm{t}}$, of the unperturbed vector field along $\Gamma_{0}$. If the Melnikov function vanishes on $\Gamma_{0}$, then there exists—under suitable nondegeneracy conditions—a family of closed orbits $\Gamma_{\eta}$ of the perturbed system (4.13) which are $\mathcal{O}(\eta)$-close to $\Gamma_{0}$ as $\eta\to 0$. Moreover, if the Melnikov function vanishes on $\Gamma_{0}$ and is positive (negative) on nearby orbits that are to the right as $\Gamma_{0}$ is traversed, then the $\Gamma_{\eta}$ are locally stable (unstable).

4.5.3 Dynamics in Regions II and III

We apply the general results of the previous section to the closed orbits of (4.9) identified at the end of Section 4.5.1: either the homoclinic orbit in the right half of the $(u,v)$ plane, or one of the periodic orbits in its interior, or one of the large-amplitude periodic orbits external to the homoclinic orbits (Figure 4.4). We assume, without loss of generality, that $\mu=1$, so $f(u)=u-u^{3}$ and $g(u,v)=(\lambda-u^{2})v$.

Let $\gamma=\{(u(t),v(t)):t\in\mathbb{R}\}$ be any of these closed orbits. To indicate a particular orbit, we label $\gamma$ by the maximum value of its first coordinate, $u(t)$, on its trajectory. We consider this label as a variable and denote it by $x$ (not to be confused with the dependent variable $x$ in the Maasch-Saltzman model). Thus, the function $\gamma:x\mapsto\gamma(x)$ is defined for all $x\in(1,\infty)$. Specifically, $\gamma(x)$ is the homoclinic orbit in the right half plane if $x=\surd 2$, a periodic orbit inside this homoclinic if $x\in(1,\surd 2)$, and a large-amplitude periodic orbit that is external to the double homoclinic if $x>\surd 2$.

Remark. There are several ways to choose an identifier for $\gamma$. For example, we could equally well have chosen the level-set value $h=H(\gamma)$, as was done in [8].

Consider the closed orbit $\gamma(x)=\{(u_{x}(t),v_{x}(t)):t\in I(x)\}$ for any $x>1$, where $I(x)=\mathbb{R}$ if $x=\surd 2$, and $I(x)$ is a period interval otherwise. The Melnikov function associated with $\gamma(x)$ is

[TABLE]

where $I_{0}$ and $I_{2}$ are defined by

[TABLE]

Here, we have used the relation $v=\dot{u}$ to convert the time integral to a contour integral on the closed orbit $\gamma(x)$.

Recall that $\gamma(x)$ is oriented clockwise; hence, Green’s theorem yields the identity $I_{0}(x)=\iint_{D(x)}\,du\,dv$, where $D(x)$ is the domain enclosed by $\gamma(x)$. Consequently, $I_{0}(x)>0$, so the condition $M(\lambda,x)=0$ is satisfied if and only if

[TABLE]

If, given $\lambda$, $M(\lambda,x)=0$ at $x=\tilde{x}$, then $M(\lambda,x)>0$ for nearby orbits that are to the right of $\gamma(\tilde{x})$ if $R^{\prime}(\tilde{x})<0$, and $M(\lambda,x)<0$ for such nearby orbits if $R^{\prime}(\tilde{x})>0$. Hence, the local stability of closed orbits is determined by the sign of $R^{\prime}(\tilde{x})$.

Since the components of the homoclinic and periodic orbits are known in terms of hyperbolic and elliptic functions, respectively, $R(x)$ can be evaluated explicitly. The computations were first done by Carr [7] and subsequently refined by Cushman and Sanders [12]. For example, for the homoclinic orbit,

[TABLE]

and $M(\lambda,\surd 2)=\tfrac{4}{3}\lambda-\tfrac{16}{15}$. It follows that $M(\lambda,\surd 2)=0$ if and only if $\lambda=\tfrac{4}{5}$. Moreover, $\lambda=\tfrac{4}{5}$ is a simple zero. Therefore, for all sufficiently small $\eta$ there exists a $\lambda(\eta)=\tfrac{4}{5}+\mathcal{O}(\eta)$ and a homoclinic orbit near $\gamma(\surd 2)$, with a symmetric result in the left half of the $(u,v)$ plane. In the $(p,r)$ plane, the homoclinic bifurcation curve is, to leading order, tangent to the line $r-1=-4(p-1)$ at the organizing center $(1,1)$; see (4.11).

The figure below shows the graph of $R:x\mapsto R(x)$ for $1<x<2$.

Starting from the values $R(1)=1$ and $R^{\prime}(1)=0$, $R(x)$ and $R^{\prime}(x)$ decrease monotonically as $x$ increases. At the homoclinic orbit (marked by a black dot), $x=\surd 2$, $R(\surd 2)=\tfrac{4}{5}$ and $\lim_{x\to\surd 2}R^{\prime}(x)=-\infty$. Beyond the homoclinic orbit, $R(x)$ decreases further, while $R^{\prime}(x)$ increases until $R^{\prime}(x)=0$; at that point (marked by a black cross), $x=x^{*}\approx 1.471$ and $R^{\prime\prime}(x^{*})>0$, so $R(x)$ reaches its minimum value $R_{\min}(x^{*})=\lambda^{*}\approx 0.752$. Beyond this point, $R(x)$ increases monotonically; $R(x)\sim x^{2}$ as $x\to\infty$.

Proofs of these statements, which do not use elliptic functions, can be found, for example, in [8]. It follows that closed orbits exist only for $\lambda>\lambda^{*}$, and they are locally stable only for $x>x^{*}$.

4.5.4 Limit Cycles in Regions II and III

The properties of the Melnikov function listed in the previous section lead to the following results for the dynamics of the perturbed system (4.9) with $\mu=1$. In all statements, it is assumed that $\eta$ is sufficiently small positive.

• •

For $\lambda>1$ (region II), there are only stable large-amplitude limit cycles which encircle the origin and pass through points $(x,0)$ with $R(x)>1$;

• •

For $\frac{4}{5}<\lambda<1$ (region IIIa), there are stable large-amplitude limit cycles which encircle the origin and unstable limit cycles in their interior which encircle the equilibrium points $(\pm 1,0)$;

• •

For $\lambda=\frac{4}{5}+\mathcal{O}(\eta)$, there is a symmetric pair of unstable homoclinic orbits, one in each half plane;

• •

For $0.752\ldots<\lambda<\tfrac{4}{5}$ (region IIIb), there is a stable large-amplitude limit cycle and an unstable large-amplitude limit cycle in its interior, both encircling the origin;

• •

At $\lambda=0.752\ldots+\mathcal{O}(\eta)$, the stable and unstable large-amplitude limit cycles join in a saddle-node bifurcation;

• •

For $\lambda<0.752\ldots$ (region IIIc), there are no limit cycles.

Thus, in addition to the homoclinic bifurcation curve found earlier, which is tangent to the line $r-1=-4(p-1)$ at the organizing center, there is a curve of saddle-node bifurcations of limit cycles, which, to leading order, is tangent to the line $r-1\approx-3.03(p-1)$ at the organizing center. This follows from (4.11), with $\mu=1$ and $\lambda^{*}\approx 0.752$. The two bifurcation curves partition the region III of Figure 4.1 into the three regions IIIa, IIIb, and IIIc, as sketched in Figure 4.2 and superimposed on the color map of Figure 4.3.

4.5.5 Limit Cycles in Region I

It remains to investigate the dynamics of the system (4.9) for $\mu<0$ ($r<p$, region I in Figure 4.1). This case is considerably simpler than the case $\mu>0$. Without loss of generality, we may assume that $\mu=-1$. The Hamiltonian is $H(u,v)=\tfrac{1}{2}v^{2}+\tfrac{1}{2}u^{2}+\tfrac{1}{4}u^{4}$. The closed orbits can again be identified by the maximum value, $x$, of its first coordinate $u(t)$, which in this case ranges over all positive values, $x>0$. The Melnikov function is given by the same expression (4.15) and vanishes if $\lambda=R(x)$, as in (4.17). In this case, both $R(x)$ and $R^{\prime}(x)$ increase as $x$ increases, so the Melnikov theory establishes that, for each $x>0$, there exists a value of $\lambda$ (given by the simple zero of the Melnikov function) such that, for some $\lambda(\eta)$ that is $\mathcal{O}(\eta)$ close to this value, the perturbed system has a unique limit cycle.

5 The Asymmetric Two-Dimensional Model

Having a complete understanding of the dynamics of the symmetric two-dimensional model (4.1), we are in a position to study the effects of symmetry breaking. The asymmetric two-dimensional model is derived formally from the Maasch–Saltzman model (3.8) by setting $q=\infty$ ($z=-x$),

[TABLE]

We will see that this system has two nondegenerate Bogdanov–Takens points, which act as organizing centers in the $(p,r)$ parameter space. The geometry of these organizing centers and the bifurcation curves emanating from them may be understood naturally as a result of the breaking of the lone $\mathbb{Z}_{2}$-symmetric Bogdanov-Takens point studied in Section 4.3.

5.1 Equilibrium States and Their Stability

The origin $P_{0}=(0,0)$ is an equilibrium state of (5.1) for all (positive) values of $p$, $r$, and $s$. If $r>p-\tfrac{1}{4}s^{2}$, there are two additional equilibrium states, namely $P_{1}=(x_{1}^{*},-x_{1}^{*})$ and $P_{2}=(x_{2}^{*},-x_{2}^{*})$, where

[TABLE]

We refer to the line $r=p-\tfrac{1}{4}s^{2}$ as the shifted diagonal (marked ‘sd’ in Figure 5.1). Note that $x_{2}^{*}<x_{1}^{*}<0$ if $p-\tfrac{1}{4}s^{2}<r<p$, and $x_{2}^{*}<0<x_{1}^{*}$ if $r>p$.

Let $P=(x^{*},-x^{*})$ be any of the equilibrium states, with $x^{*}=0$, $x_{1}^{*}$, or $x_{2}^{*}$. The Jacobian of the vector field at $P$ is

[TABLE]

The system is linearly stable at $P$ if the trace is negative, $-1+r-(x^{*})^{2}<0$, and the determinant is positive, $p-r+3(x^{*})^{2}+2sx^{*}>0$.

The stability results are illustrated in Figure 5.1.

We see that the parameter space $\mathbb{R}_{+}^{2}$ is partitioned into six regions, which depend on $s$. Referring to the labels in Figure 5.1, these regions are

[TABLE]

Summarizing the results of the stability analysis, we find that

• •

$P_{0}$ is stable in regions Oa and Ob, undergoes a supercritical Hopf bifurcation on ‘e0’ with natural frequency $\omega_{0}^{*}=\sqrt{p-1}$;

• •

$P_{1}$ is stable in region III, undergoes a subcritical Hopf bifurcation on ‘e1’ with natural frequency $\omega_{1}^{*}=\sqrt{2(r-1)+s\sqrt{r-1}}$; and

• •

$P_{2}$ is stable in regions Oa, III, and IIIo, undergoes a subcritical Hopf bifurcation on ‘e2’ with natural frequency $\omega_{2}^{*}=\sqrt{2(r-1)-s\sqrt{r-1}}$.

The introduction of asymmetry ($s>0$) results in two changes. The vertical line $\{p=1,r>1\}$, which is the locus of Hopf bifurcations for $P_{1}$ and $P_{2}$ in the symmetric case (Figure 4.1), unfolds into a parabola. The vertex of this parabola is at the point $(1,1)$, and the parabola is tangent to the shifted diagonal $r=p-\tfrac{1}{4}s^{2}$ at the point $(1+\tfrac{1}{2}s^{2},1+\tfrac{1}{4}s^{2})$. As we will see in Section 5.2, both these points are organizing centers. For convenience, we label them

[TABLE]

As $(p,r)$ moves across the shifted diagonal in the direction of decreasing $p$, the equilibrium states $P_{1}$ and $P_{2}$ emerge in a saddle-node bifurcation. If the crossing occurs above $Q_{2}$, $P_{1}$ and $P_{2}$ are both unstable; if it occurs below $Q_{2}$, $P_{1}$ is unstable while $P_{2}$ is stable. In the region Oa, $P_{0}$ and $P_{2}$ co-exist as stable equilibria, while $P_{1}$ is an unstable equilibrium. On the diagonal for $p<1$, $P_{0}$ and $P_{1}$ exchange stability in a transcritical bifurcation.

5.2 Organizing Centers

We make the change of variables $(x,-(x+y))\mapsto(x,y)$ as in Section 4.3. In the new coordinate system, (5.1) becomes

[TABLE]

The equilibrium points are $P_{0}=(0,0)$, $P_{1}=(x_{1}^{*},0)$, and $P_{2}=(x_{2}^{*},0)$, where $x_{1}^{*}$ and $x_{2}^{*}$ are again given by (5.2).

Let $P=(x^{*},0)$ be any of the equilibrium points. The Jacobian of the vector field at $P$ is

[TABLE]

If $x^{*}=0$, the Jacobian has a double-zero eigenvalue at $Q_{1}$ for any $s$, and if $x^{*}=x_{1}^{*}$ or $x^{*}=x_{2}^{*}$, it has a double-zero eigenvalue at $Q_{2}$. Hence, the introduction of asymmetry causes the organizing center to unfold into a center at $Q_{1}$ associated with $P_{0}$ and a center at $Q_{2}$ associated with $P_{1}$ and $P_{2}$.

Figure 5.2 shows the bifurcation curves emanating from the two organizing centers for $s=0.8$ (the value chosen by Maasch and Saltzman). They were computed with the AUTO continuation package.

There are three Hopf bifurcation curves (shown in red), two emanating from $Q_{1}$ and one emanating from $Q_{2}$: (i) a parabolic curve to the left of $Q_{1}$, where $P_{1}$ undergoes a subcritical Hopf bifurcation, (ii) a horizontal line $\{r=1,p>1\}$ to the right of $Q_{1}$, where $P_{0}$ undergoes a supercritical Hopf bifurcation; (iii) a parabolic curve to the right of $Q_{2}$, along which $P_{2}$ undergoes a subcritical Hopf bifurcation. These three curves are the same as the ones identified in the local analysis, cf. Section 5.1. In addition, there are three homoclinic bifurcation curves (shown in blue), two emanating from $Q_{1}$ (solid blue) and one emanating from $Q_{2}$ (dashed blue). These curves are identified by an unfolding procedure, as in Section 4.5. Since the points $Q_{1}$ and $Q_{2}$ are both nondegenerate Bogdanov–Takens points, there are no other bifurcation curves besides the Hopf and homoclinic bifurcation curves emanating from them.

Figure 5.3 shows the phase portraits at the six small black diamond markers along the vertical line $p=1.55$ in Figure 5.2. The color scheme is as follows: The flow of (5.1) is shown as blue streamlines. The black and red curves correspond to the stable and unstable limit cycles, respectively. The equilibrium states $P_{0}$, $P_{1}$, and $P_{2}$ are indicated by black, red, and green markers, respectively.

• •

Frame (a), $r=1.2$: a stable limit cycle surrounds the unstable equilibrium point $P_{0}$.

• •

Frame (b), $r=1.45$: the unstable equilibrium points $P_{1}$ and $P_{2}$ exist but are unstable, a stable limit cycle surrounds $P_{1}$ and $P_{2}$.

• •

Frame (c), $r=1.6$: the equilibrium points $P_{1}$ and $P_{0}$ have switched positions, $P_{2}$ has become stable, an unstable limit cycle surrounds $P_{2}$.

• •

Frame (d), $r=2.0$: the unstable limit cycle has disappeared in the homoclinic bifurcation.

• •

Frame (e), $r=2.5$: a large-amplitude unstable limit cycle exists inside the large-amplitude stable limit cycle.

• •

Frame (f), $r=3.0$: the large-amplitude stable and unstable limit cycles have disappeared in a saddle-node bifurcation; $P_{2}$ is the only attractor.

6 Summary

The purpose of this chapter is to show an interesting application of dynamical systems theory to a problem of climate science. The object of investigation is a model of the Pleistocene climate, proposed by Maasch and Saltzman in 1990. The model is a conceptual model designed specifically to explain the persistence of glacial cycles during the Pleistocene Epoch. The Milankovitch theory of orbital forcing establishes a correlation between glacial cycles and periodic oscillations in the Earth’s orbit around the Sun, but orbital forcing by itself is insufficient to explain the observed temperature changes.

The Maasch–Saltzman model incorporates a feedback mechanism that is driven by greenhouse gases, in particular atmospheric $\mathrm{CO}_{2}$. The model is based on plausible physical arguments, and preliminary computational experiments indicate that it reproduces several salient features of the Pleistocene temperature record. In this article, the emphasis is on the internal dynamics of the model. We focused on the prevalence and bifurcation properties of limit cycles in the various parameter regimes in the absence of external forcing.

The Maasch–Saltzman model is formulated in terms of the anomalies of the total global ice mass, the atmospheric $\mathrm{CO}_{2}$ concentration, and the volume of the North Atlantic Deep Water (NADW, a measure of the strength of the North Atlantic overturning circulation). It consists of three differential equations with four parameters and is difficult to analyze directly. Our results indicate how one can obtain fundamental insight into its complex dynamics by first considering a highly simplified two-dimensional version. The dimension reduction is achieved by (formally) letting one of the parameters—representing the rate of change of the volume of NADW relative to that of the total global ice mass—tend to infinity. The approximation is justified by the observation that the NADW changes on a much faster time scale than the total global ice mass.

The two-dimensional model has two primary parameters, $p$ and $r$, which are both positive, and one secondary parameter $s$, which reflects the asymmetry between the glaciation and deglaciation phases of the glacial cycles. By first ignoring the asymmetry $(s=0)$, we obtained a complete understanding of the dynamics and the persistence of limit cycles.

Figure 4.2, which is a sketch of the various bifurcation curves in the $(p,r)$ parameter space, summarizes the main results. The origin is an equilibrium state $P_{0}$ for all $(p,r)$; $P_{0}$ is stable in region O. In addition, there are two equilibrium states $P_{1}$ and $P_{2}$ in regions II and IIIa-c; they are generated in a pitchfork bifurcation along the diagonal $r=p$ and are stable in region IIIa-c. Stable limit cycles exist in regions O, I, II, IIIa, and IIIb. They are created in supercritical Hopf bifurcations along the boundary between regions O and I. In region I, they surround the unstable equilibrium state $P_{0}$, and in the other regions they surround all three equilibrium states. Along the boundary between regions II and IIIa, a pair of unstable limit cycles, each surrounding one of the two stable equilibrium states $P_{1}$ and $P_{2}$, are created in a subcritical Hopf bifurcation. These newborn limit cycles grow in amplitude until they become homoclinic orbits at the boundary between regions IIIa and IIIb; in region IIIb, they have merged into one large-amplitude unstable limit cycle, which surrounds the three equilibrium states and sits just inside the large-amplitude stable limit cycle. Finally, the stable and unstable large-amplitude limit cycles merge and disappear in a saddle-node bifurcation as $(p,r)$ crosses the lower boundary of region IIIb into region IIIc. All these results have been confirmed computationally; see Figure 4.3.

Having thus obtained a complete understanding of the dynamics of the symmetric simplified model, we then re-introduced asymmetry ($s>0$). A comparison of Figure 5.2 with Figure 4.2 shows the effects of symmetry breaking. The main change is that the single organizing center, which governed the dynamics in the symmetric case, splits into two organizing centers. Also, the curves of homoclinic bifurcations and saddle-node bifurcations of limit cycles become more complex.

The complexity of the bifurcation diagram of Figure 5.2 gives an indication why it is difficult to analyze the dynamics of the Maasch–Saltzman model directly. It also justifies our approach of first analyzing the highly simplified model and then gradually relax the constraints that were imposed to derive the simplified model [17].