Asymptotic analysis of internal relaxation-oscillations in a conceptual climate model

TL;DR

This paper develops a simplified climate model demonstrating that internal relaxation-oscillations, akin to ice-age cycles, can occur naturally without external forcing, supported by asymptotic analysis and paleoclimatic data comparison.

Contribution

It introduces a new dynamical system based on the KCG climate model, analyzing its stability and oscillatory behavior, and derives a formula for oscillation periods.

Findings

The model exhibits stable limit cycles with relaxation-oscillations.

The period formula aligns with paleoclimatic data.

Oscillations occur without external forcing across various parameters.

Abstract

We construct a dynamical system based on the KCG (K\"all\'en, Crafoord, Ghil) conceptual climate model which includes the ice-albedo and precipitation-temperature feedbacks. Further, we classify the stability of various critical points of the system and identify a parameter which change generates a Hopf bifurcation. This gives rise to a stable limit cycle around a physically interesting critical point. Moreover, it follows from the general theory that the periodic orbit exhibits relaxation-oscillations which are a characteristic feature of the Pleistocene ice-ages. We provide an asymptotic analysis of their behaviour and derive a formula for the period along with several estimates. They, in turn, are in a decent agreement with paleoclimatic data and are independent of any parametrization used. Whence, our simple but robust model shows that a climate may exhibit internal relaxation-oscillations without any external forcing and for a wide range of parameters.

Figures

Tables

Table 1. Table 1: Symbols used in the paper. The typical values are based on [ 14 , 9 , 8 ] .

Symbol	Meaning	Typical value
$T$	Globally averaged temperature	$-$
$L$	Southward ice sheet extent	$-$
$t$, $\tau$	Time variables	$-$
$Q$	Solar constant	1361 W m-2
$\gamma$	Continent to ocean area ratio	$0.3$
$A$	Budyko constant in OLR flux	-267.96 W m-2
$B$	Budyko constant in OLR flux	1.74 W m-2 K-1
${\tau }_0$	Ice sheet yield stress	0.3$\times {10}^5$ Pa
${\rho }_i$	Ice density	0.92$\times {10}^3$ kg m-3
$H$	Ice sheet height scale	2.1 m${}^{\frac{1}{2}}$
$s$	Snow line slope	0.4 $\times {10}^{-3}$
$m$	Ablation rate	$0.5$ m year-1
$c$	Atmosphere thermal capacity	${10}^7$ J m-2 K-1
$h_0$	Height of the snow line over Arctic ocean	$1.2\times {10}^3$ m
$\kappa$	Nondimensional height of the snow line over Arctic ocean	0.1
${\alpha }_c$	Continental albedo	$-$
${\alpha }_o$	Oceanic albedo	$-$
$\xi$	Ratio of accumulation to ablation	$-$
${\lambda }_0$	The boundary between accumulation and ablation zones	$-$
${\alpha }_1$, ${\alpha }_2$	Parameters of the continental albedo	0.25 and 1, respectively
${\alpha }_{-}$, ${\alpha }_{+}$	limits of oceanic albedo	0.6 and 0.22, respectively
${\xi }_{-}$, ${\xi }_{+}$	limits of the ratio of accumulation and ablation	0.1 and 0.5, respectively
${\theta }_{\alpha }$, $\Delta \alpha$	translation and steepness parameters for oceanic albedo	1.4 and 0.15, respectively
${\theta }_{\xi }$, $\Delta \xi$	translation and steepness parameters for $\xi$	1.39 and 0.025, respectively
$\sigma$	An arbitrary sigmoid function	$-$

Equations

$$\theta =\frac{4B}{Q}T,\lambda =\frac{s^2}{H^2}L,\tau =\frac{2}{3}\frac{ms}{H^2}t\text{where}H:=\sqrt{\frac{4{\tau }_0}{3{\rho }_{i}g}},$$

$$\beta :=-\frac{4A}{Q},\mu :=\frac{3}{2}\frac{B{H}^2}{msc},\kappa :=\frac{s{h}_0}{H^2}.$$

$$\left\{\begin{array}[]{l}\dfrac{d\theta}{d\tau}=\mu\left(1-\beta-\gamma\alpha_{c}(\lambda)-(1-\gamma)\alpha_{o}(\theta)-\theta\right),\vspace{6pt}\\ \dfrac{d\lambda}{d\tau}=\sqrt{\lambda}\left(\left(1+\xi(\theta)\right)\lambda_{0}(\lambda,\theta)-1\right),\end{array}\right.\quad\text{for}\quad\theta,\lambda>0,$$

$${\lambda }_0(\lambda ,\theta ):=\frac{1}{\lambda }\left(-\left(\kappa (\theta )+\lambda +\frac{1}{2}\right)+\sqrt{\kappa (\theta )+2\lambda +\frac{1}{4}}\right).$$

$${\alpha }_c(\lambda ):=\frac{1}{2}\left({\alpha }_0+{\alpha }_1+({\alpha }_0-{\alpha }_1){\sigma }_c\left(\frac{\lambda -{\lambda }_{\alpha }}{\Delta \lambda }\right)\right),{\alpha }_0\le {\alpha }_1,$$

$${\alpha }_o(\theta ):=\frac{1}{2}\left({\alpha }_{+}+{\alpha }_{-}+({\alpha }_{+}-{\alpha }_{-}){\sigma }_o\left(\frac{\theta -{\theta }_{\alpha }}{\Delta \alpha }\right)\right),{\alpha }_{-}\ge {\alpha }_{+},$$

$$\xi (\theta ):=\frac{1}{2}\left({\xi }_{+}+{\xi }_{-}+({\xi }_{+}-{\xi }_{-}){\sigma }_{\xi }\left(\frac{\theta -{\theta }_{\xi }}{\Delta \xi }\right)\right),{\xi }_{-}\le {\xi }_{+}.$$

$$\begin{array}{r}{\lambda }_0(\lambda )\ge 0\text{(ice sheet not stagnant)},\\ \lambda \le 1\text{(maximal ice sheet size)},\\ 2\lambda +\kappa \ge 0\text{for}\kappa <0\text{(ice sheet sufficiently large)}.\end{array}$$

$$F(\theta ,\lambda )=\mu \left(1-\beta -\gamma {\alpha }_c(\lambda )-(1-\gamma ){\alpha }_o(\theta )-\theta \right),G(\theta ,\lambda )=\sqrt{\lambda }\left(\left(1+\xi (\theta )\right){\lambda }_0(\lambda )-1\right).$$

$$\begin{array}{rl}\frac{\partial F}{\partial \theta }& =-\mu \left((1-\gamma )\frac{d{\alpha }_o}{d\theta }+1\right),\frac{\partial F}{d\lambda }=-\mu \gamma \frac{d{\alpha }_c}{d\lambda },\\ \frac{\partial G}{\partial \theta }& =\sqrt{\lambda }\left(\frac{d\xi }{d\theta }{\lambda }_0+(1+\xi )\frac{\partial {\lambda }_0}{\partial \theta }\right),\frac{\partial G}{\partial \lambda }=\frac{1}{2\sqrt{\lambda }}\left(\left(1+\xi \right){\lambda }_0-1\right)+\sqrt{\lambda }(1+\xi )\frac{\partial {\lambda }_0}{\partial \lambda },\end{array}$$

$$\frac{\partial {\lambda }_0}{\partial \theta }=-\frac{1}{\lambda }\frac{d\kappa }{d\theta }\left(-1+\frac{1}{2\sqrt{\kappa +2\lambda +\frac{1}{4}}}\right),\frac{\partial {\lambda }_0}{\partial \lambda }=\frac{1}{{\lambda }^2}\left(\kappa +\frac{1}{2}+\frac{\kappa +\lambda +\frac{1}{4}}{\sqrt{\kappa +2\lambda +\frac{1}{4}}}\right).$$

$$\begin{array}{rl}h(\theta )& ={\alpha }_c^{-1}\left(-\frac{1}{\gamma }\left((1-\gamma ){\alpha }_o(\theta )+\theta +\beta -1\right)\right),\\ k^{\pm }(\theta )& =\frac{1}{2}\frac{(1+\xi )\xi }{(2+\xi {)}^2}\left(1-2\left(1+\frac{2}{\xi }\right)\kappa \pm \sqrt{1-4\left(1+\frac{2}{\xi }\right)\kappa }\right).\end{array}$$

$${\lambda }_0=\frac{1}{1+\xi },$$

$$\kappa \le \frac{1}{4}\frac{\xi }{2+\xi }.$$

$$\frac{dh}{d\theta }=-\frac{1}{\gamma {\alpha }_c^{\prime }}\left((1-\gamma )\frac{d{\alpha }_o}{d\theta }+1\right).$$

$$\Delta \alpha \text{ small enough for }\frac{dh}{d\theta }=0\text{ to have two solutions }{\theta }_{-}<{\theta }_{+}.$$

$${\mu }_c=-\sqrt{\lambda }\left(1+\xi \right)\frac{d{\lambda }_0}{\partial \lambda }{\left((1-\gamma )\frac{d{\alpha }_o}{d\theta }+1\right)}^{-1}\text{at}({\theta }_c,{\lambda }_c).$$

$$r_{\pm }=\frac{1}{2}\left(\text{tr}J\pm \sqrt{{\left(\text{tr}J\right)}^2-4\text{det}J}\right).$$

$$\begin{array}{rl}\text{tr}J(\mu )& =-\mu \left((1-\gamma )\frac{d{\alpha }_o}{d\theta }+1\right)+\sqrt{\lambda }\left(1+\xi \right)\frac{d{\lambda }_0}{\partial \lambda },\text{vspace}12pt\\ \text{det}J(\mu )& =\mu \gamma \sqrt{\lambda }\left(1+\xi \right)\frac{d{\alpha }_c}{\partial \lambda }\frac{\partial {\lambda }_0}{\partial \lambda }\left(\frac{dh}{d\theta }-\frac{d{k}^{\pm }}{d\theta }\right)\end{array}\text{at}(\theta ,\lambda )=({\theta }_c,{\lambda }_c),$$

$$\begin{array}{rl}\frac{dh}{d\theta }& =-\frac{\partial F}{\partial \theta }{\left(\frac{\partial F}{\partial \lambda }\right)}^{-1}=-\frac{1}{\gamma {\alpha }_c^{\prime }}\left((1-\gamma )\frac{d{\alpha }_o}{d\theta }+1\right),\\ \frac{d{k}^{\pm }}{d\theta }& =-\frac{\partial G}{\partial \theta }{\left(\frac{\partial G}{\partial \lambda }\right)}^{-1}=-\frac{1}{{\lambda }_0}\frac{{\lambda }_0^2\frac{dk}{d\theta }+\frac{\partial {\lambda }_0}{\partial \theta }}{(1+\xi )\frac{\partial {\lambda }_0}{\partial \lambda }}\end{array}\text{at}(\theta ,\lambda )=({\theta }_c,{\lambda }_c).$$

$${\left(\text{tr}J(\mu )\right)}^2-4\text{det}J(\mu ).$$

$$r_{+}=\frac{2\text{det}J}{\text{tr}J-\sqrt{{\left(\text{tr}J\right)}^2-4\text{det}J}},$$

$${\lambda }_0(\lambda )=1-4\lambda +O({\lambda }^2)\text{as}\lambda \to 0^{+}.$$

$$\left(1+\xi (\theta )\right){\lambda }_0(\lambda )-1=1-4\lambda +\xi (\theta )-1+O\left({\lambda }^2\right)+O\left(\xi {\lambda }^2\right)=\xi (\theta )-4\lambda +O\left({\lambda }^2\right)+O\left(\xi {\lambda }^2\right),$$

$${\alpha }_c(\lambda )={\alpha }_0+{\alpha }_1\lambda .$$

$$\left\{\begin{array}[]{l}\dfrac{d\theta}{d\tau}=\mu\left(1-\beta-\gamma\alpha_{0}-\dfrac{1}{2}(1-\gamma)(\alpha_{-}+\alpha_{+})-\gamma\alpha_{1}\lambda+\dfrac{1}{2}(1-\gamma)(\alpha_{-}-\alpha_{+})\sigma\left(\dfrac{\theta-\theta_{\alpha}}{\Delta\alpha}\right)-\theta\right),\vspace{6pt}\\ \dfrac{d\lambda}{d\tau}=\dfrac{1}{2}(\xi_{-}+\xi_{+})\sqrt{\lambda}\left(1+\dfrac{\xi_{+}-\xi_{-}}{\xi_{-}+\xi_{+}}\sigma\left(\dfrac{\theta-\theta_{\xi}}{\Delta\xi}\right)-\dfrac{8}{\xi_{-}+\xi_{+}}\lambda+O\left(\lambda^{2}\right)+O(\xi(\theta)\lambda^{2})\right),\end{array}\right.$$

$$\begin{array}{r}a=1-\beta -\gamma {\alpha }_0-\frac{1}{2}(1-\gamma )({\alpha }_{-}+{\alpha }_{+}),b=\frac{1}{8}\gamma {\alpha }_1({\xi }_{-}+{\xi }_{+}),\\ c=\frac{1}{2}(1-\gamma )({\alpha }_{-}-{\alpha }_{+}),d=\frac{{\xi }_{+}-{\xi }_{-}}{{\xi }_{-}+{\xi }_{+}},\nu =\frac{1}{8}\sqrt{2({\xi }_{-}+{\xi }_{+})}\gamma {\alpha }_1\mu ,x_{\alpha }={\theta }_{\alpha },x_{\xi }={\theta }_{\xi },\end{array}$$

$$x=\theta ,y=\frac{8}{{\xi }_{-}+{\xi }_{+}}\lambda ,t=\sqrt{2\left({\xi }_{-}+{\xi }_{+}\right)}\tau ,$$

$$\left\{\begin{array}[]{l}\dfrac{dx}{dt}=\nu\left(f(x)-y\right),\vspace{6pt}\\ \dfrac{dy}{dt}=\sqrt{y}\left(g(x)-y\right),\end{array}\right.$$

$$f(x)=\frac{1}{b}\left(a+c\sigma \left(\frac{x-x_{\alpha }}{\Delta \alpha }\right)-x\right),g(x)=1+d\sigma \left(\frac{x-x_{\xi }}{\Delta \xi }\right).$$

Full text

Keywords: dynamical system, conceptual climate model, relaxation-oscillations, matched asymptotics

1 Introduction

Conceptual climate models provide a feasible mean of grasping the most important mechanisms of the complex system as the climate itself [29, 4]. Focusing on only the most essential features they can help to understand the basic laws that govern the dynamics of various parameters of the climate. Of course, by construction, they are unable to predict the temperature changes in a great detail (as opposed to General Climate Models (GCMs)) but still are invaluable tool in understanding the Earth system evolution [22]. Their virtue is that they usually are low-dimensional dynamical systems that can be understood by analytical or numerical means without the use of the supercomputers. The analytical approach has the advantage of being able to simultaneously tackle an infinite number of initial conditions providing results that can then be verified with the use of GCMs.

One of the remarkable features of the paleoclimatic record is the emergence of relaxation-oscillations of the climate [17, 4]. These asymmetrical variations of the ice sheet extent indicate that there should exist a nonlinear mechanism based on various feedbacks that can drive such dynamics. More specifically, a slow growth of the ice sheet is followed by a rapid deglaciation [17]. Conceptual models that try to explain this phenomena are reviewed in [4].

In the literature there is a number of approaches to conceptual climate modelling with varying complexity and descriptive variables used. This work is influenced by Källén, Crafoord and Ghil’s model that analyses temperature and ice sheet evolution which essentially are the energy and mass balances [14, 10, 11]. The beginnings of formulating these balances as differential equations can be traced back to Budyko [3], Sellers [30] and Weertman [34]. The Budyko-Sellers model show that there is a hysteric response of the climate to the solar forcing which indicates that the climate can abruptly enter a completely remote cold or warm state (a tipping point). On the other hand, one of the main features of the KCG model indicates that the climate can oscillate without any external (astronomical) forcing as happens in the Milankovitch theory [23, 2, 12]. As a matter of fact, the periodic change in Earth’s axis tilt, precession and eccentricity proposed by Milankovitch to describe glacial episodes are nowadays treated as a pacemaker of the nonlinear oscillator being the climate itself [5, 31]. Similar considerations lead Fowler and his collaborators to incorporate the energy and mass balances into the chemical reaction framework following CO2 evolution [8, 7]. As authors show, the relaxation-oscillations of the climate can be realised as an inclusion of a mechanism responsible for formation of a progracial lakes.

Another conceptual climate framework was developed by Saltzman and collaborators [18]. They focus on describing the dynamics of total ice volume, CO2 concentration and deep-ocean temperature. This model has been investigated for many years and the whole theory is summarized in [29] (but see also [6]). The important feature of that model is the emergence of relaxation-oscillations which describe the Pleistocene ice ages.

More mathematically inclined considerations on conceptual climate models have been conducted by McGehee and collaborators. In [20, 32] they have investigated the full version of the Budyko-Sellers model stated as a partial differential equation. If its steady state is approximated by quadratic functions the whole dynamics can be captured into a low-dimensional dynamical system analysed in [21]. The model is constructed with a discontinuous vector field and possesses periodic orbits [33].

In this work, which is a continuation of [28], we investigate a generalization of the KCG model and show that for a reasonable parameter regime it predicts that the climate will undergo relaxation-oscillations. First, we state the model in its full generality and classify all its critical points according to the relative slopes of the nullclines. This is a reasonable geometrical approach since the model possesses a multitude of parameters which specify the position of the critical points. Their exact position is not necessarily needed to know as opposed to their stability or bifurcations. After investigating the general model we approximate it by focusing only on the physically most interesting critical point. It appears that for almost the whole range of the bifurcation parameter the system undergoes a relaxation-oscillation without any external forcing. With the use of matched asymptotic expansions we find the asymptotic formula for its period and amplitude. Further, we show that the former can be bounded by some elementary functions that produce very reasonable estimates on the glacial ages. The lower estimate appears to be completely independent from the particular choices of the functional forms of the albedo and accumulation/ablation. This makes it a universal bound. We close the paper with a numerical verification of our results and find that the asymptotic formula is decently accurate even for small values of the parameter.

2 Model statement and general results

2.1 Main dynamical system

In [28] we have provided a thorough derivation of the generalized version of the KCG model firstly proposed in [14]. Here, we will just summarize the main features of the model and an interested Reader is invited to learn about all the details in the original works.

The main constituents of the KCG model are the energy and mass balances. The energy equation is of Budyko-Sellers type being a competition between incoming short-wave and outgoing long-wave radiation (OLR). Earth acquires its energy from the Sun in means of the high quality light which some part is reflected and some absorbed by the planet. The ratio of the reflected to incident energy is called the albedo $\alpha$ and its present average value for Earth is around $0.3$ (see [9]). Moreover, within the naturally occurring features of the surface of our planet, the fresh snow has one of the highest values of $\alpha$ being $0.8-0.9$ while open ocean one of the lowest equalling $0.06$. Following KCG idea we will decompose albedo into two parts one representing the continents $\alpha_{c}$ and the other the oceans $\alpha_{o}$. The continental part describes the reflectivity of the land on which the ice sheet can form and advance. The oceanic albedo, in turn, depends mainly on the temperature since for cold climates the sea ice forms increasing the amount of light reflected. The variability of the albedo with respect to the amount of ice present is called the ice-albedo feedback.

The second ingredient of the model is the mass balance written as to describe the advances and retreats of the north hemisphere’s ice sheet. This has firstly been proposed in [34]. Assuming zonal, i.e. longitudinal, symmetry the ice sheet’s northern margin is in the Arctic ocean and the southward flow of the ice is taken to be perfectly plastic. The mass can change twofold: either by nourishment with a snowfall or by melting due to approach into the ablation zone at low latitudes. The imaginary boundary between these two regimes is called the snow line. It essentially is the $0^{\circ}$ C isotherm and it can be visualized as a slanted line which height over the surface increases southward (meridionally). That is, high altitude near the equator has a comparable temperature to the low altitudes near the North pole. The accumulation and ablation regions are thus demarcated by the intersection of the snow line and the surface of the ice sheet. The connection with the energy equation comes from the fact that the rate of the snowfall depends on the temperature. Higher temperature increases evaporation of the oceans and hence increases the amount of water vapour in the atmosphere. Some of it can fall as a precipitation at higher latitudes as a snow building up the ice sheet. This mechanism is called the precipitation-temperature feedback.

For the neatness of the notation we will write the model in the nondimensional form the very beginning and describe what scalings have been used. Let $\theta$ be the temperature, $\lambda$ the southward ice sheet extent and $\tau$ denote the time. They are related to the dimensional values by

[TABLE]

and all appearing constants are explained in Tab. 1. Moreover, define the following parameters

[TABLE]

Note that $\kappa$ can be negative and may depend on the temperature as an nondecreasing function. This is due to the fact that the height of the snow line over the Arctic Ocean $h_{0}$ changes according to the temperature. However, we will not specify the concrete form of $\kappa$ (in [11] it is taken to be a linear function). The nondimensionalized model can now be written as follows

[TABLE]

where the boundary between accumulation and ablation zones is given by

[TABLE]

Since $\alpha_{c}(\lambda)$, $\alpha_{o}(\theta)$ and $\xi(\theta)$ are monotone and bounded [14] it is useful to define the following well-known class of functions that can be used to represent those climatic quantities.

Definition 1**.**

A function $\sigma\colon\mathbb{R}\rightarrow\mathbb{R}$ is called sigmoid if it is bounded and differentiable with a non-negative derivative. As a normalization one can take $\lim\limits_{x\rightarrow\pm\infty}\sigma(x)=\pm 1$.

Then, it is useful to take

[TABLE]

where $0\leq\alpha_{0,1}\leq 1$ are limits of $\alpha_{c}(\lambda)$ for $\theta\rightarrow\pm\infty$, $\lambda_{\alpha}$ is the translation and $\Delta\lambda$ being the steepness parameter. Similarly, we define

[TABLE]

and

[TABLE]

Notice that $\alpha_{c}$ and $\xi$ are increasing while $\alpha_{o}$ is a decreasing function. These parametrizations of the climatic features are somewhat arbitrary since the physical processes that govern them are extremely complex and impossible to fully implement at the level of a conceptual model. As we will see in the sequel, many features of the relaxation-oscillations are independent of the particular choice of the sigmoid function. This makes the model more robust. Moreover, the various parameters can be tweaked in many ways to represent different climatic scenarios. It is not our aim to exactly describe reality but to show how aforementioned mechanisms give rise to an interesting dynamical behaviour which is also seen in the real world. Our choice of parameters is thus illustrative with a strong connection to the climate.

Lastly, we would like to remark about the physical validity of the model (3). The conditions for these are the following

[TABLE]

The physical derivation of the above is given in [28]. In the first and third of the above cases it is possible to provide different equations that build a meaningful model in these situations.

2.2 The phase plane

Here, we will provide a characterization of the dynamical behaviour of (3). First, denote the vector fields in (3) as

[TABLE]

Computing the derivatives we have

[TABLE]

where we have suppressed the explicit writing of arguments for the brevity of notation. Further, the derivatives of the snow line $\lambda_{0}$ have the form

[TABLE]

Note that thanks to the third assumption in (8) the function $\lambda_{0}$ is $\theta$-nonincreasing. On the other hand, by the analysis of $d\lambda_{0}/d\lambda$ we see that for fixed $\kappa$ the function $\lambda_{0}(\cdot,\kappa)$ has a global maximum for $\lambda=\kappa(1+4\kappa)/2$ when $\kappa\geq 0$ and $\lambda=-\kappa/2$ when $\kappa<0$.

The nullclines $h$ and $k$ of $\theta$ and $\lambda$ respectively, can be readily computed in a closed form giving

[TABLE]

The continental albedo is a monotone function of $\theta$ hence it has a well-defined inverse. Moreover, the $\lambda$-nullcline is a solution of the quadratic equation

[TABLE]

hence it has two branches denoted by $\pm$. They join at the common point for $\lambda=\kappa(1+4\kappa)/2$ when $\kappa\geq 0$ or $\lambda=-\kappa/2$ when $\kappa<0$. At that point $k^{\pm}$ has a singular derivative. Obviously, $k^{+}\geq k^{-}$ and the sufficeint and necessary condition for they to exist is

[TABLE]

since then the square root in (12) is real. Inverting the above yields an important result that no ice sheet can exist for sufficiently high temperatures [34].

The shape of the $\lambda$-nullcline depends on the particular form of $\kappa$. In [28] we have shown that $k^{-}(\theta)=O(\kappa^{2})$ when $\kappa\rightarrow 0$ hence it becomes physically irrelevant as soon as the snow line is low enough. On the other hand, the other branch $k^{+}$ converges in that case to $(1+\xi)\xi/(2+\xi)^{2}$ which is an increasing function. This is physically the most important branch at which the critical points can be situated. We will come back to this issue in the next section.

Since $\alpha_{c}^{-1}$ is increasing, the overall shape of the $\theta$-nullcline is not distorted by a composition with it. It is easy to observe that since $\alpha_{o}$ is a sigmoid function it is almost constant for arguments far to the left or far to the right of $\theta_{\alpha}$. In these regions, the dominant contribution comes from the decreasing linear function in the definition of $h$ in (12). On the other hand, in the vicinity of $\theta_{\alpha}$ the sigmoid dominates making the whole graph of $h$ to resemble the S-shape. More precisely, computing the derivative we have

[TABLE]

Wee see that the above has two zeros only if the equation $(1-\gamma)\alpha_{o}^{\prime}=-1$ has two solutions. Since the derivative of a sigmoid function has a single maximum which height is determined by the steepness parameter $\Delta\alpha$, the nullcline $h$ has two extrema for sufficiently small $\Delta\alpha$. Therefore, we further assume that

[TABLE]

Whence, $h$ increases on $(\theta_{-},\theta_{+})$ and decreases otherwise. Now, we are ready to state the stability result.

Theorem 1**.**

Assume (14) and (16). If $(\theta_{c},\lambda_{c})$ is a fixed critical point of the system (3) then the following holds.

• •

If $(\theta_{c},\lambda_{c})$ lies on $k^{-}$ then it is unstable.

• •

If $(\theta_{c},\lambda_{c})$ lies on $k^{+}$ then

• –

if $h$ is decreasing at that point then it is stable iff $dh/d\theta<dk^{+}/d\theta$,

• –

if $h$ is increasing at that point then it becomes unstable when $\mu$ increases through $\mu_{c}$. Moreover, if $dh/d\theta>dk^{+}/d\theta$ the Hopf bifrucation occurs. Here,

[TABLE]

Proof.

The stability of a particular critical point is determined from the eigenvalues $r_{\pm}$ of the Jacobi Matrix of the partial derivatives of $F$ and $G$ evaluated at the critical point $(\theta_{c},\lambda_{c})$. Denote this matrix by $J$ and then by elementary formulas we have

[TABLE]

Computing the trace and determinant at the critical point we obtain

[TABLE]

where we have used the Inverse Function Theorem to note that the derivatives of the nullclines are given by

[TABLE]

Notice that since $\lambda_{0}(\cdot,\kappa)$ has a maximum for fixed $\kappa$ and $\lambda$ on the upper branch of the nullcline, we have $\partial\lambda_{0}/\partial\lambda<0$ on $k^{+}$ and vice-versa.

Assume now that the critical point lies on the upper branch $k^{+}$. If $h$ is decreasing at that point, we have $\partial F/\partial\theta<0$ and $\partial\lambda_{0}/\partial\lambda<0$. Therefore, $\text{tr}J<0$ for every $\mu>0$. On the other hand, the square root in the determinant is pure real or imaginary according to the sign of the following quadratic

[TABLE]

Although the exact zeros of the above can be readily computed, they are not needed for the stability result. For if $\mu$ is chosen to yield a real square root we have $r_{-}<0$ and

[TABLE]

which due to (9) and (19) has its sign governed by the relative slope of the nullclines at the critical point. Therefore, by the Hartman-Grobman’s Theorem the critical point is stable if $\det J>0$ which is $dh/d\theta<dk^{+}/d\theta$ and we have a node. Moreover, if the slope of $h$ is larger than this of $k^{+}$ the critical point is a saddle. On the other hand, in the case when the square root in $\det J(\mu)$ is imaginary, the critical point is always stable because $\text{Re}\,r_{\pm}=\text{tr}J(\mu)/2<0$ (a focus).

Next, assume that at the critical point $h$ is increasing and consider the upper branch of $\lambda$-nullcline. Now, since $\partial F/\partial\theta>0$ the trace is positive for $\mu>\mu_{c}$. The reasoning essentially the same as above applies to that case yielding a change of stability at $\mu_{c}$. When $dh/d\theta>dk^{+}/d\theta$ at the critical point we have $\text{tr}J(\mu_{c})=0$ and the square root in (18) is purely imaginary. Therefore, as $\mu$ increases through $\mu_{c}$ the eigenvalues cross the imaginary axis on the complex plane with a non-zero speed because $\text{tr}J(\mu)$ is a linear function. The Hopf Bifurcation Theorem [27] yields the result.

Completely analogous reasoning can be applied to the lower branch $k^{-}$ and we omit the details. The important thing to keep in mind is that we cannot have $dh/d\theta>dk^{-}/d\theta$ at the critical point when $h$ is decreasing. For then the nullclines cannot intersect and there is no critical point. A careful sign counting helps to see that the lower branch is unstable. The proof is now complete. ∎

We have thus seen that since the lower branch of the $\lambda$-nullcline is completely unstable, the physically meaningful phenomena can happen only on $k^{+}$. Now, we will elaborate on this particular case.

3 Relaxation-oscillations

3.1 Model simplification

In order to facilitate the analysis of relaxation-oscillations present in (3) we make several simplifying and reasonable assumptions. First, since the nondimensional snow line height at the Arctic is usually small [8] we take the first order approximation with $\kappa=0$. As we have mentioned above, this forces the lower branch of the $\lambda$-nullcline to vanish leaving the upper one as a increasing function of $\theta$. Perhaps, taking a static snow line is not completely physically feasible but treating it as a first approximation is justified when the variations of the temperature are small. Therefore, our subsequent considerations can be thought as meaningful only when we are investigating vicinity of the critical point representing our climate. The model is not suitable for describing large excursions from it.

As was also shown in [28], the ice sheet extent $\lambda$ is a small number. Specifically, we always have $\lambda<2/9$ but in reality the upper bound is even smaller. In that case the snow line can be expanded into Taylor series for $\kappa=0$

[TABLE]

Moreover, from the definition of $\xi$ and the model data in Tab. 1 we can see that its magnitude is small letting us to approximate even further

[TABLE]

as $\lambda\rightarrow 0^{+}$ and $\xi\rightarrow 0^{+}$. Further, we take the simplest physically sensible form of the continental albedo

[TABLE]

Although it is not a sigmoid function, the small variations in it allow us to use the linear approximation. This is consistent with the assumption that $\kappa=0$. Now, we can write (3) as

[TABLE]

as $\lambda\rightarrow 0^{+}$ and $\xi\rightarrow 0^{+}$. If we neglect the higher order terms, define the following

[TABLE]

and rescale the variables by

[TABLE]

we obtain the following dynamical system

[TABLE]

with the nullclines

[TABLE]

Notice that now $x,y$ and $t$ are of order of unity. For simplicity and without any lose of generality we have taken the sigmoid functions $\sigma$ to be of the same family.

We have completed the derivation of the simplified model of climate under the assumption the the variations of the snow line are small enough for us to take $\kappa=0$. We can think of this description as a magnification of the vicinity of the present climate oscillations.

It is easy to note that $g$ is an increasing sigmoid function while $f$ has a typical S-shape having one minimum and one maximum (see Fig. 1) given by

[TABLE]

It is also convenient to decompose $f$ into three branches

[TABLE]

where subscripts denote the stable and unstable parts of $f$. On their respective domains the branches are monotone and hence invertible. More precisely

[TABLE]

where the vertical bars indicate the restriction of the domain. Moreover, for future convenience we define the following points

[TABLE]

Before we move to the analysis we have to make several assumptions concerning the parameters involved. Specifically, in order for the relaxation-oscillations to occur there should exist a unique unstable critical point which we denote by $x_{c}$. It suffices then to make the following assumptions.

[TABLE]

where the stability result can be proved essentially in the same way as in the previous section while the critical parameter is given by

[TABLE]

To get a quick peek on the reasonable size of the critical parameter we can make a crude estimate that a ice sheet oscillates between $70^{\circ}$ and $45^{\circ}$ of latitude when the temperature changes by at most $20^{\circ}$C (this is an overestimate). Then, assuming that $x_{\pm}$ bounds the amplitude of the oscillations and approximating $f^{\prime}(x_{c})$ by its finite difference we can estimate that $\nu_{c}\approx 0.15$. This can be considered as a small quantity. Conducting the exact calculations with our data it occurs that $\nu_{c}=0.1$.

Now, we can make some simple geometrical observations. From the Fig. 1 we can see that for a physically reasonable parameter choice the function $g$ is almost constant for $x<x_{-}$ and $x>x_{+}$. Moreover, near $x_{\xi}$ it is well-approximated by a linear function. The latter property holds also for $f$ near $x_{\alpha}$. These observations can help us to find a reasonable approximation of the critical point.

Proposition 1**.**

Let (35) be satisfied. Then

[TABLE]

which for $x_{\alpha}=x_{\xi}=x_{0}$ can be written as

[TABLE]

where $C=\frac{1}{3}\max\left\{|f^{\prime\prime\prime}(x)|,|g^{\prime\prime\prime}(x)|:x\in[x_{-},x_{+}]\right\}$.

Proof.

The proof goes by straightforward calculations. First, observe that if we expand $f$ and $g$ in Taylor series at the inflection points $x_{\alpha}$ and $x_{\xi}$ we have

[TABLE]

Equalling both expansions we obtain

[TABLE]

For the remainder it holds

[TABLE]

for the constant $C$ defined in the assertion. This finishes the proof. ∎

In practice the above approximation can be decently accurate and for our data yields an error of $3\times 10^{-4}$ which is smaller than necessary in this application.

3.2 Matched asymptotic analysis

Next, we proceed to the main topic. By using the general theory [24, 16] it follows that the dynamical system (29) possesses a stable limit cycle for all values of $\nu\geq\nu_{c}=0.1$. Our general considerations from previous section indicate that the Hopf bifurcation occurs at $\nu_{c}$ and a cycle is born with a period of $200\times 10^{3}$ years which is twice the data-estimated value of the Pleistocene ice age oscillations. This indicates that, if the model is correct, $\nu$ should be larger than $\nu_{c}$. With the increase of $\nu$ the amplitude of oscillations increases until they become relaxation-oscillations. It means that the climate oscillates for almost the whole parameter range, that is the relaxation-oscillations are a robust phenomenon. Below we will use the matched asymptotic analysis to find their leading order behaviour. These results can be compared with the Van der Pol oscillator (see [15] for construction of the approximation and [19] for a rigorous justification). The following method borrows from [26] while in [25] an algorithimc approach has been developed with rigorous proofs of convergence. A simplified version that captures the quintessence of our calculations is described in [13].

We will approximate the solution of (29) is several regions of the phase plane and the following material will be organized according to the particular layer. In order to track the progress it is helpful to consult Fig. 1. Moreover, as a notational convenience we will retain the same letter for denoting the exact solution with its leading order approximation. This is the only order to which we will match boundary layers and thus it will not cause any confusion. However, it will enhance the readability greatly. All the results of the following analysis are summarized in Theorem 2 below.

First, introduce the parameter $\epsilon=\nu^{-1}$ and write (29) as a singularly perturbed system

[TABLE]

where we assume that $\epsilon\ll 1$. From the phase plane we see that there are two asymptotically stable branches of the slow manifold $y=f(x)$ (in the sense of Tikhonov-Levinson theory [26]).

Suppose that we start our evolution on the left branch of the slow manifold where $g(x)<f(x)$ and we continue our oscillation sliding down-right. We will subsequently consider various time scales and the behaviour of $(x(t),y(t))$ on them.

Outer layer (left). This is the slow phase of the oscillations. The leading order behaviour of (45) emerges when we set $\epsilon=0$ and and then the equations constitute the reduced system

[TABLE]

where the subscript denotes the outer approximation. This is a differential-algebraic system saying that the evolution takes place on the slow manifold $y=f(x)$. By differentiating the first equation and using the second we arrive at

[TABLE]

Notice that this equation has a singularity at $x_{o}=x_{-}$ which has to be resolved. To this end, let $x_{o}(t)=x_{-}-\varphi(t)$ where $0<\varphi(t)\ll 1$. Then

[TABLE]

The leading order behaviour of $\varphi$ is thus given by

[TABLE]

Notice that we have incorporated the integration constant into the time variable setting $t=0$ precisely at the jumping point $x_{0}=x_{-}$. Since at that point the outer approximation becomes non-smooth, we have to find another set of equations describing the solution for $t>0$.

Transition layer. The singularity of the outer solution (49) suggests there should exist a stretching transformation that produces the distinguished limit. Hence, we introduce the transition layer variables

[TABLE]

where $p,q,r$ have to be found. Plugging into (45) gives

[TABLE]

Now, we have to find the appropriate balance between various terms in the above. For the outer solution the derivative was assumed to be small but eventually blew-up, hence we expect that the correct balance will be between all three terms in the equation for $x_{t}$. Since $f(x_{-}+\epsilon^{p}x_{t})-f(x_{t})=f^{\prime\prime}(x_{-})\epsilon^{2p}x_{t}^{2}/2+O(\epsilon^{3p})$ and we have to keep the derivative of $y_{t}$ it follows that the only choice of parameters is $1+p-r=2p=q$ and $q-r=0$. Therefore

[TABLE]

Whence, the dynamical system under this stretching becomes

[TABLE]

In the leading order the second equation integrates to a linear function and the first one becomes

[TABLE]

where the integration constant has been incorporated into $\delta_{t}(\epsilon)$. This nonlinear equation is of Riccati type and can be solved with a substitution

[TABLE]

In this new variable the differential equation has the form

[TABLE]

which is Airy’s equation having a solution

[TABLE]

with Ai and Bi being Airy’s functions. From the known asymptotic behaviour of these (see [1]) we can see that the only choice for $x_{t}$ to match $x_{i}$ for $t_{t}\rightarrow-\infty$ is to choose $F_{t}=0$. In that case we have

[TABLE]

If we had retained $F_{t}$ an unnecessary minus sign would appear in front of the above asymptotic expansion. In order to conduct the match with the outer approximation we write

[TABLE]

and $|t|\ll\epsilon^{-2/3}$. In order to satisfy (49) we take $\delta_{t}(\epsilon)=0$.

The matching with outer layer is complete and we proceed to investigate behaviour for large $t_{t}$ where a matching with the inner layer is anticipated. A potential difficulty arises since (58) has a singularity at the first zero of Airy’s function, that is to say at

[TABLE]

Moreover, from the properties of Ai we know that the aforementioned zero is simple and

[TABLE]

Now, we will see how does this behaviour agree with the next layer.

Inner layer. Naturally, we proceed to investigate the inner layer approximation. To this end introduce the fast time scale $t_{i}=\epsilon^{-1}\left(t-\delta_{i}(\epsilon)\right)$ where $\delta_{i}(\epsilon)$ is a yet unknown translation indicating the beginning of the layer. This leads to

[TABLE]

where now subscript denotes the inner approximation. The leading order behaviour emerges when we put $\epsilon=0$ leaving only the first equation

[TABLE]

where $Y(\epsilon)$ is the integration constant that may depend on $\epsilon$ with an order smaller than one. Now, observe that if the above equation were supposed to yield a leading order approximation to the solution of (45) we should have $x_{i}\rightarrow x_{-}$ as $t\rightarrow-\infty$ for $\epsilon=0$. Therefore we put $Y=f(x_{-})-\epsilon^{\frac{2}{3}}W$ and obtain

[TABLE]

Notice that we anticipated the first meaningful power of $\epsilon$ in the asymptotic expansion. This is due to the fact that the $2/3$ exponent emerges in the transition layer approximation to which inner has to match for $t_{i}\rightarrow-\infty$.

Next, focus only on the leading order case for which $\epsilon=0$. The above equation has two critical points: $x_{-}$ and $\widetilde{x}_{-}$ (see Fig. 1). The approach to the latter is exponential while to the former only algebraic since $f^{\prime}(x_{-})=0$. More specifically, setting $x_{i}=x_{-}+\varphi$ with $0<\varphi\ll 1$ lets us determine that

[TABLE]

In (65) the integration constant has been incorporated into $\delta_{i}(\epsilon)$.

In order to match inner and transition layers we use (65) with $t_{i}$ expressed in terms of $t_{t}$, that is to say $t_{i}=\epsilon^{-\frac{1}{3}}(t_{t}-\epsilon^{-\frac{2}{3}}\delta_{i}(\epsilon))$. And hence

[TABLE]

Remembering the stretching transformation (50) we can successfully match $x_{t}$ and $x_{i}$ if $\delta_{i}(\epsilon)=\epsilon^{2/3}t_{t}^{*}$ which fixes the time shift. The matching with the higher order requires retaining $\epsilon$ in (64) and plugging $x_{i}(t_{i})=x_{-}+\widetilde{\varphi}(t_{i})$ with $0<\widetilde{\varphi}\ll 1$. This yields

[TABLE]

This equation can be solved to give

[TABLE]

where $V$ is a time shift constant. The asymptotic expansion of the above requires some care. We are interested in letting $t_{i}\rightarrow-\infty$, i.e. when the solution leaves the inner layer to the left. On the other hand, when $\epsilon$ goes to zero the left asymptote of $\tan$ function approaches minus infinity. This suggests that we should expand as follows

[TABLE]

as $t_{i}\rightarrow-V-\epsilon^{-1/3}\pi/\left(2\sqrt{f^{\prime\prime}(x_{-})W/2}\right)$. Now, expressing the above in the transition layer time $t_{t}$ we have

[TABLE]

as $\epsilon\rightarrow 0$ and $|t_{t}|\ll\epsilon^{1/3}$, where we have conducted the irrelevant time shift with an appropriate choice of constant $V$ and used the previously determined value of $\delta_{i}(\epsilon)=\epsilon^{-2/3}t_{t}^{*}$. If we compare the above expansion with (61) we see that they are compatible if we choose

[TABLE]

We have thus completed the matching of inner layer with the transition approximation.

We move further and investigate the behaviour of $x_{i}$ for large $t_{i}$. Due to (64) and (71) we see that it approaches the following limit

[TABLE]

To determine the pace of the approach we set $x_{i}=x_{++}-\psi$ with $0<\psi\ll 1$ to have

[TABLE]

where the integration constant has been incorporated into $C_{i}$. Note that $f^{\prime}(x_{++})<0$. Next, we have to check if this behaviour matches the right branch of the outer layer.

Outer layer (right). It is convenient to introduce the following right outer variables

[TABLE]

Similarly as before in (49) we can show that $x_{p}$ has has a singular derivative at $x_{+}$. Moreover, by setting $x_{p}=x_{++}+\varphi$ we can show that

[TABLE]

where integration constant has been incorporated into $\delta_{c}(\epsilon)$. We can see that the above linear form cannot be matched to the exponential behaviour of the inner layer (73). Therefore, there exists a final layer that joins these two approximations.

Corner layer. In this layer we introduce the following scalings

[TABLE]

where we have explicitly included the previously determined time lag. Moreover, we can assume that $\delta_{c}(\epsilon)=o(\epsilon^{2/3})$ as $\epsilon\rightarrow 0$. The use of the same auxiliary constants $p,q$ and $r$ should not confuse the Reader with those used in the transition layer. The system (45) now becomes

[TABLE]

Since now $f^{\prime}(x_{++})\neq 0$ we can obtain the distinguished limit by requiring thet $1+p-r=p=q$ and $q-r=0$. Therefore,

[TABLE]

and the system (45) becomes

[TABLE]

To the leading order we have

[TABLE]

where the integration constant has been incorporated into $\delta_{c}(\epsilon)$. The above linear equation can be solved yielding

[TABLE]

where $C_{c}$ is integration constant. Since $f^{\prime}(\widetilde{x}_{-})<0$ we have the following asymptotic behaviours

[TABLE]

and

[TABLE]

Since $t_{c}=t_{i}-\epsilon^{-1}\delta_{c}(\epsilon)$ we thus have

[TABLE]

Remembering that $x=\widetilde{x}_{-}+\epsilon\,x_{c}$ in order to match corner and inner layers, i.e. the above equation with (73), we have to take

[TABLE]

We can see that $\delta_{c}(\epsilon)\ll\delta_{i}(\epsilon)$ and hence the solution spends less time in the corner layer than in the two previous ones.

On the other hand, for the matching inner with right outer layer we express the time scales as $t_{c}=\epsilon^{-1}(t_{p}+\delta_{p}(\epsilon)-\delta_{c}(\epsilon)-\epsilon^{2/3}t_{t}^{*})$ and write (82) in the form

[TABLE]

The consistency with (75) can be achieved by taking

[TABLE]

because $\delta_{c}(\epsilon)$ is of higher order.

So far the solution have travelled from $x_{-}$ on the left branch of the stable manifold to $x_{++}$ on the right one. It left outer layer and went through transition, inner and corner layers in time $\epsilon^{2/3}t_{t}^{*}+o(\epsilon^{2/3})$ as $\epsilon\rightarrow 0$. The further journey takes the solution to the larger fold $x_{+}$ and then the evolution is mirror reflected with a change of $x_{-}$ to $x_{+}$ and $x_{++}$ to $x_{--}$. Therefore, in the last part of the proof we have to calculate the time spent on the right branch of the stable manifold. Denote it by $T_{p}$ and integrate (47) with $x_{o}$ replaced by $x_{p}$. Therefore,

[TABLE]

where $W$ is defined in (71). Manipulating further we can separate the leading order behaviour

[TABLE]

The solution has now arrived at $x_{+}$ from where it can return to the left stable branch of the manifold. The time necessary for this travel can be found essentially as before giving (90). The amplitudes of the limit cycle is clearly given by $x_{++}-x_{--}$ and $f(x_{+})-f(x_{-})$. Whence, we have proved the following theorem.

Theorem 2**.**

Let (35) be satisfied. Then the period of the relaxation-oscillations has the following asymptotic expansion

[TABLE]

as $\nu\rightarrow\infty$, where

[TABLE]

and $-\zeta\approx-2.338$ is the largest zero of Airy’s function Ai. Moreover, the $x-$ and $y-$amplitudes of the cycle are

[TABLE]

and

[TABLE]

as $\nu\rightarrow\infty$.

From the above calculations we can immediately see that it can be conducted to any arbitrary $f$ and $g$ satisfying the usual geometric conditions. The specific forms (30) are irrelevant for all the considerations. Therefore, we have the following corollary.

Corollary 1**.**

The assertions of Theorem 2 hold for any twice-differentiable functions $f$ and $g$ where $f$ has a single non-degenerate minimum and maximum while $g$ is sigmoid.

Furthermore, it turns out that the period can be easily approximated by elementary functions.

Corollary 2**.**

Let (35) be satisfied. Then, the period of relaxation-oscillations for sufficiently large $\nu$ can be bounded

[TABLE]

where

[TABLE]

with

[TABLE]

Moreover, for sufficiently large $\nu$ the amplitude of $x=x(t)$ satisfies

[TABLE]

Proof.

We will find lower and upper estimates on the leading order term in the asymptotic formula for the period (90). Notice that in (90) the leading order term is written in terms of the integrals $I_{\pm}$ given in (91). Since $g$ is bounded due to the properties of sigmoid functions, the term $g\circ f_{s\pm}^{-1}$ can be bounded from below and above by $1\pm d$. This simplifies the integral yielding

[TABLE]

which follows by simple calculation. Therefore, since $g$ is increasing and $f_{s+}^{-1}(w)\geq x_{+}$ for $w\geq f(x_{-})$ we have

[TABLE]

and similarly

[TABLE]

Adding both above inequalities yields the assertion.

For the estimate on the amplitude observe that $f$ approaches its asymptotes from within, that is to say

[TABLE]

by the properties of sigmoid functions. Therefore, from (72) it follows that $x_{++}\leq f_{x_{+}}^{-1}(f(x_{-}))<\widehat{x}_{++}$ where $\widehat{x}_{++}$ is the abscissa of the point of intersection of the asymptote with a value $f(x_{-})$. Completely analogous reasoning concerning the left branch of $f$ implies the assertion. This ends the proof. ∎

The above corollary gives us some useful formulas for determining the period od the oscillations. Notice that in order to find the lower bound we need only four numbers: two local extrema of $f$ and two bounds of $\xi$. By Theorem 2 we know that the leading order term of $y-$amplitude of the oscillations is equal to the difference between these two extrema. Whence, knowing the ice sheet extent and its physical properties we are able to estimate the period od ice ages. Note that this result assumes only two essential mechanisms: ice-albedo and temperature-precipitation feedbacks. It is a universal formula independent of any parametrization! For our data we have

[TABLE]

which gives a very reasonable estimate of the oscillations during the last one million years where period is approximately $100\times 10^{3}$ years. Of course, by manipulating the parameters we can obtain other sensible outcomes of our model. However, the main point is that a simple two-dimensional model can provide a realistic leading-order approximation to the very complex natural system.

Notice that formulas for both period estimates in (94) contain logarithmic terms. It is a simple observation that $T_{+}$ can diverge to infinity when $g(x_{\pm})\rightarrow f(x_{\pm})$. However, this approach is logarithmic hence very slow and, therefore, is in principle resistant to small variations in the model parameters. Since not every parameter of our model is directly measurable, this is a valuable feature of the formula for the period that allows for certain error.

An exemplary illustration of the limit cycle of (29) is given on Fig. 2. We can see that the range of the ice sheet extent shows realistic values while the variation of the temperature is somewhat too large. This is probably due to the initial simplification of an immobile snow line, i.e. setting $\kappa=0$ (which will be removed in our future work). However, the period lies very close to the real-world value: for $\nu=10$ it is equal to $121\times 10^{3}$ y. The most noteworthy observation is the presence of asymmetric relaxation-oscillations. This is one of the most characteristic features of the Pleistocene ice-ages. As we mentioned above, our model exhibits such a behaviour for all $\nu>\nu_{c}=0.1$.

The numerical verification of the asymptotic formula (90) is depicted on Fig. 3. We have calculated the period of the limit cycle for $\nu$ changing from $1$ to $10^{7}$. The results are presented on two plots: the first concerning the small values of $\nu$ and the second in log-log scale for the large values of the bifurcation parameter. For the latter we have plotted the difference between the numerical value and the leading order term of (90), denoted by $\overline{T}_{ro}$, versus the $\nu^{-2/3}$ term. We can see that the accuracy is very good even for $\nu$ of order of unity. The power-law behaviour is clearly seen. Surprisingly, the formula is decently accurate even for small values of $\nu$.

4 Conclusion and future work

There are several important remarks concerning our generalization of the KCG model. As we have seen, assuming two climate feedback mechanisms, i.e. ice-albedo and precipitation-temperature, leads to self-sustained oscillations after a Hopf bifurcation which happens for very small parameter magnitude. The astronomical forcing is not necessary to produce the dynamics of a realistic period. Moreover, the mathematical formulations of the involved feedbacks can be arbitrary as long as they satisfy natural assumptions of monotonicity and boundedness. This frees the model from unnecessary ad-hoc choice of the complex system parametrizations and furnishes a more robust model.

In the vicinity of the Hopf-bifurcating critical point a stable limit cycle can exist provided there are no other critical points nearby. This let us to build an approximation of the periodic orbit of the model and analyse its relaxation-oscillations. Remarkably, a simple but physical model such as ours can exhibit a complex behaviour resembling the Pleistocene glaciations without any additional mechanism added. We have provided an asymptotic formula for the oscillation period and verified that it is decently accurate for a wide range of the bifurcation parameter. Although the exact closed-form formula for the leading order term in the period cannot be found, a simple estimates can be readily derived. The lower one is independent of the particular form of the parametrizations while the higher is insensitive to measurement errors. Both of them yield realistic estimates of the internal climate oscillations.

In our future work we will try to incorporate several other degrees of freedom and analyse the corresponding dynamical system. These may include a slowly varying $\nu$ used to model the Middle Pleistocene Transition where the period of glaciations changed from $40$ to $100$ thousands of years. Another point of interest will be to include the viscoelastic response of the lithosphere under the load of the ice sheet. We also plan to investigate the role of the carbon cycle.