An analysis of the periodically forced PP04 climate model, using the theory of non-smooth dynamical systems
AAn analysis of the periodically forced PP04 climatemodel, using the theory of non-smooth dynamicalsystems.
Kgomotso S. Morupisi a , Chris Budd* a a Department of Mathematical Sciences, University of Bath, BA2 7AY, UK
Abstract
In this paper we perform a careful analysis of the forced PP04 model for climatechange, in particular the behaviour of the ice-ages. This system models thetransition from a glacial to an inter-glacial state through a sudden release ofoceanic Carbon Dioxide into the atmosphere. This process can be cast in termsof a Filippov dynamical system, with a discontinuous change in its dynamicsrelated to the Carbon Dioxide release. By using techniques from the theory ofnon-smooth dynamical systems, we give an analysis of this model in the casesof both no insolation forcing and also periodic insolation forcing. This revealsa rich, and novel, dynamical structure to the solutions of the PP04 model. Inparticular we see synchronised periodic solutions with subtle regions of existencewhich depend on the amplitude and frequency of the forcing. The orbits can becreated/destroyed in both smooth and discontinuity induced bifurcations. Westudy both the orbits and the transitions between them and make comparisonswith actual climate dynamics.
Keywords:
Climate models, ice ages, PP04 model, non-smooth dynamics,Filippov systems
1. Introduction
Reduced climate models (RCMs), see for example [EKKV17], [SM90][SM91],[Pai01], [AD15], [WWHM16], [Cru12],[KE13],[Dij13], have been used exten-sively to study various forms of climate dynamics. Whilst not in any way asubstitute for general climate models (GCMs) for an accurate simulation ofclimate dynamics from which predictions can be made, they are nonethelessvery useful for investigating certain types of qualitative climate phenomena,particularly those that occur over time scales which are too long for a realistic ∗ Corresponding author
Preprint submitted to Elsevier May 19, 2020 a r X i v : . [ n li n . C D ] M a y alculation on a GCM. RCMs can be used both to give insights into the macro-scopic behaviour of certain types of climate phenomena, and also as ways oftesting the predictions of the GCMs. In this paper we will in particular performa careful analysis of the RCM usually called the PP04 model [PP04] which hasbeen used both to gain insight into the past behaviour of the Earth’s glacialcycles (ice-ages) and also to predict future glacial events [Cru13], [ADCvdH18].Glacial cycles themselves show very subtle dynamics, with an interplay of varia-tions of ice, temperature and of Carbon Dioxide, all coupled together by variousfeedback loops, and with external forcing from the Sun. These cycles have ledto periodic, and significant, variations in the temperature, ice cover and CarbonDioxide levels of the Earth. In the most recent glacial periods (over the lasthalf a million years) these cycles have a roughly 100 kyr periodicicty, whereasbefore that period the cycles has a shorter period of around 40 kyr. It is gen-erally believed [IBB + The emphasis of this paper will the study of the solutions of the PP04 modelwhen it is driven by (quasi-)periodic insolation forcing.If the magnitude of the insolation forcing is zero we will show that the modeladmits periodic solutions, associated with natural internal time-scales for thegrowth and retreat of ice sheets coupled to Carbon Dioxide levels in the atmo-sphere. These solutions have a period of 147 kyr ≡ /ω . In the non-smoothPP04 model these orbits are created and destroyed in border collision bifurca-tions, arising when certain fixed points intersect a discontinuity surface. In asmoothed version of the they arise instead through Hopf bifurcations from thesteady state followed by cyclic saddle-nodes.If the insolation forcing is purely periodic, of frequency ω and amplitude µ ,we find that if µ , | ω − nω /m | (with n, m = 1 , , .. ) are both small, then a2 ode-locked periodic solution exists which is a perturbation of the unforcedsolution. This solution is synchronised to the insolation forcing, but in generalhas a different phase. For small µ , the regions of existence in the ( µ, ω ) spaceare linear tongues for all values of n if m = 1, and are bounded by saddle-node bifurcations. Outside of these tongues we see quasi-periodic motion whichcombines the insolation forcing frequency and the natural frequencies of thesystem. As the amplitude µ increases, the mode locked solution persists until ittypically loses stability at a grazing bifurcation. At this point we see complexand multi-modal behaviour. As µ increases, the tongues for different values of n can overlap and expand, leading to the co-existence of different period states.Close to the boundaries of these regions we see a variety of different types ofbehaviours (with subtle domains of attraction), leading to interesting transitionsbetween the states as parameters are varied, with some behaviour having aqualitative resemblance to that at the Mid-Pleistocene Transition. When theperiodic insolation forcing is extended to being quasi-periodic we then see thestable periodic solutions perturbing to invariant tori.The layout of the reminder of this paper is as follows. In Section 2 we willbriefly review some of the observed features of the glacial cycles which we seekto reproduce in our models. In Section 3 we will review some of the existingmodels and will motivate the PP04 model for glacial dynamics. In Section 4we will explain the basic ideas, and necessity, of using non-smooth dynamicalsystems theory in the analysis of the PP04 model, showing that it takes the formof a Filippov system without sliding. In Section 5 we will study the unforcedPP04 model. We will analyse the changes in the dynamics as parameters vary,looking at both the existence of fixed points, and of periodic solutions, and thetransitions between them as a result of border collision bifurcations in the non-smooth system. In Section 6 we will give a detailed mathematical analysis ofthe existence, and stability, of the periodic solutions of the PP04 model underthe effects of both small, and large, periodic insolation forcing. In Section 7 wewill support these results through a series of numerical computations of the Ω − limit set of the solutions using a Mont´e-Carlo method, and also the transitionsbetween different states. Finally in Section 8 we will discuss the implications ofthese results to our understanding of the dynamics of the climate.
2. Observed climate dynamics and glacial cycles
It is a feature of observed climate dynamics over the last few million years,that the Earth experiences glacial cycles, which are roughly periodic variationsbetween hot and cold states. The hot (or interglacial states) tend to last forrelatively short periods compared to the longer (or glacial) states. The periodof these in the past 800,000 years has been roughly 100k years. Before then theperiod was closer to 40k years. The change between these types of behaviouris called the Mid-Pleistocene Transition (MPT), and has been studied by manyauthors [PP04, SM91, SM90, Pai98]. We see this in the following two figures.3he first shows the changes in ice volume and temperature over the last millionyears. The second shows the MPT.
Figure 1: Ices age temperatures and ice volume taken from the Vostok and EPICA ice cores.(Image from en.wikipedia.org .) The changes to the past Earth climate, shown in these figures, can be studiedthrough paleo-data sources such as coral reefs, deep sea sediments,continentaldeposits of flora and fauna and ice cores [JLP + + + +
87] particularly through the similarities between themid-June insolation forcing at 65 N and data from the δ O isotope [PJR + + δ O isotope, showed that the fluctuations of volume of ice had experiencedperiods of 23 kyr and 41 kyr which supports the contribution of orbital forc-4 igure 2: The oscillation cycles obtained from an oxygen isotope δ O . (Image from Lisieckiand Raymo 2005). ing, more precisely the variations in precession and obliquity, as the cause ofthese oscillations [Hel82], [Pai17]. The oscillations of the glacial cycles withperiods of 23 kyr and 41 kyr have been successfully reproduced by the inclusionof astronomical forcing. However, the dominant asymmetrical 100 kyr periodoscillations seen in the above figures observed for the last 450 kyrs [DSCW13]have been difficult to explain through astronomical variation theory alone. Itseems clear that a full explanation of this periodicity must involve considera-tions of internal processes occurring on the Earth, such as the time-scales forthe advance and retreat of the ice sheets. Hence in an attempt to address thisproblem, different conceptual models using the physics of ice sheets and theocean-atmosphere feedback have been developed. We now briefly review someof these.
3. Dynamical models for climate change
Whilst many sophisticated models for climate change exist, such as the GlobalClimate Models, for example the HadGem3 model [KJea18], these cannot berun to simulate the long periods associated with the glacial cycles. Hence, inorder to obtain insight into the glacial cycles it is often useful to make use ofsimpler conceptual models, which study variations of ice sheets and effects ofastronomical forcing on these ice sheets. These are usually expressed in termsof low-dimensional dynamical systems. Various such dynamical systems modelshave been proposed to explain the glacial cycles, and a good review of these isgiven in the papers [Cru12, DSCW13] and books [KE13, Dij13]. Many of thesemodels make use of some form of relaxation oscillator to explain the observed5ehaviour of the ice ages, and in particular the mechanism of repeated slowgrowth of ice sheets followed by rapid decay of ice sheets. In these models phaselocking is often observed between the astronomical variations and the internalvariability of the Earth’s climate [Dij13].
Smooth models:
Many such models use the concept of smooth excitable sys-tems [SM91, SM90, Cru12] and explain the origin of the glacial cycles throughmechanisms such as the Hopf bifurcation. In particular, Saltzman and Maasch[SM90, SM91], advocated that the glacial cycles could be viewed as limit cy-cles synchronized by the insolation forcing. The Saltzman and Maasch (SM90)[SM90] and 1991 (SM91) [SM91] models adopt the hypothesis that the increasein insolation causes a decrease in ice sheet mass and that the change in at-mospheric carbon is driven by tectonic forcing. The SM90 and SM91 modelsare smooth dynamical systems with non-linearity on the Carbon Dioxide equa-tions playing an important role in inducing the existence of limit cycle. Thesemodels interpret the Mid-Pleistocene Transition (MPT) as a bifurcation froma quasi-linear response to a nonlinear resonance [Cru12] with the SM90 modelexperiencing a Hopf bifurcation [AD15].
Non-smooth threshold models:
A second type of models, described in [Pai98,PP04, AD15, WWHM16], use the concept of thresholds to link the Carboncycle to the retreat of the ice sheets. The justification for this choice being thepresence of abrupt transitions in the paleo-climate geological records [Pai01].In 1998 Paillard (P98) suggested that the climate system can be representedas three quasi-stable states that are driven by astronomical forcing. In P98,threshold criteria are used to bring about the instability into the systems so thatit can switch between different climatic states. The presence of such thresholdsin the P98 model give it the form of a hybrid dynamical system [BBCK08].Gildor and Tziperman [GT00] proposed a (higher complexity) box model ofthe climate with thresholds to explain the MPT. This model (GT2000) coupledocean, atmosphere, sea ice and land ice behaviour, with the ocean divided intoeight boxes, the atmosphere into four vertically averaged boxes and with the seaice responding to the energy balance equations. In this model the sea ice wasconsidered to have a hysteretic response to the variations in the land ice volumewith thresholds in the land ice volume (due to growth and melting) bringingabout the switching mechanism. This mechanism suggested that the MPT couldbe as a result of climate cooling which in turn allowed sea ice cover to expand,hence activating the sea ice switch. Therefore implying that the glacial cyclesof 100 kyr timescale did not rely on the astronomical forcing. Another modelof glacial cycles is given in [AD15] (AD15). This model proposes that global icevolume relaxes to an equilibrium state depending on a climatic state and thatmelting of the ice sheet is governed by astronomical variations of the insolationforcing. In this model the climatic states are governed by a drift function whichdescribes a nonlinear relationship between the state and the ice volume, and thetransitions from the 40 kyr to the 100 kyr period states are described as a trans-critical bifurcation on the slow manifold. Non-smooth effects, and an analysis6ased on non-smooth dynamics, are also considered in the paper [WWHM16]which looks at a glacial ice-line model formulated as a Filippov system.
The model that we will study in this paper was introduced by Paillard and Par-renin in 2004 [PP04] (PP04) and describes the evolution and feedback mecha-nisms associated with the global ice volume V , the extent of Antarctic ice sheet A and atmospheric Carbon Dioxide content C . PP04 is a piece-wise smoothmodel of the glacial cycles that incorporates physical mechanisms involving theinfluence of the Antarctic ice-sheet extent on the bottom water formation. Thesein turn cause dramatic changes in the amount of atmospheric Carbon Dioxideduring the glacial-interglacial transitions when Carbon Dioxide is thought to bereleased from the deep ocean. Fuller details of the physical motivation of thismodel can be found in [PP04]In the PP04 model the equations for the change in V depend on the amount ofatmospheric Carbon Dioxide C and involve the astronomical forcing, and theextent of the Antarctic ice sheets A is then coupled with global ice volume, withfull details of the model given in [PP04]. The amount of Carbon Dioxide in theatmosphere is considered to depend on the reduction of the amount of globalice volume, the insolation forcing, and the state of the Southern Ocean. Inparticular the ocean contribution is represented by the (discontinuous) Heav-iside function H ( − F ) and is dependent on the the ’salty bottoms efficiency’parameter F ( V, A, C ) which is positive when the climate is in a glacial state,and negative when it is in an inter-glacial state. In this model the atmosphericCarbon Dioxide rapidly increases when the Southern Ocean suddenly ventilates.The ventilation process is described in [PP04], and occurs when the deep oceanstratification (which usually prevents the water mixing) ceases due to a diffi-culty in salty bottom water formation, leading to a release of Carbon Dioxide.This release of Carbon Dioxide then leads to a warming of the Earth whichdrives a rapid deglaciation process. After this event the ice sheets accumulateslowly, until the threshold value of F is again reached and another ventilationis triggered. Consequently, during the glacial periods, there is no ocean contri-bution in the system until another release of Carbon Dioxide from the ocean isinitiated.Perhaps the most important part of this model is the inclusion of the function F , which acts as a switch in the system between the glacial and inter-glacialstates. According to [PP04] F should increase when changes in V lead to globalcooling, and decrease when continental shelf areas are reduced. The function F ( V, A, C ) is then defined by F = aV − bA − cI ( t ) + d. (1)In this model F increases with global ice volume and decreases with the Antarc-tic ice sheet and Southern Hemisphere insolation forcing. The constant param-eter d controls the threshold crossing for the model from glacial to interglacial7tates, and I , is the daily insolation forcing at 60 S . The value of c is taken in[PP04] to be very small. The full PP04 model for the glacial cycle model is thendefined by the piece-wise smooth, non-autonomous dynamical system given by: dVdt = ( − xC − yI ( t ) + z − V ) τ V , (2) dAdt = ( V − A ) τ A , (3) dCdt = ( αI ( t ) − βV + γH ( − F ) + δ − C ) τ C , (4)with H being the Heavyside function defined by H ( − F ) = (cid:40) F < , F > . Here x, y, z, α, β, γ, δ are physical constants, τ V , τ A , τ C are the time-scales asso-ciated with ice formation and Carbon Dioxide growth, and I ( t ) is the insola-tion at 65 ◦ North. According to Garcia-Olivares et al [GOH13], the parameter δ can be interpreted as the Carbon Dioxide reference level and β representsthe positive feedback between the temperature and the Carbon Dioxide levelsthrough the ice volume V . Values for I ( t ) are provided in Mitsui et. al.[MA14], [DSCW13, MA14]. As a result of the introduction of the ocean contri-bution into the equation for Carbon Dioxide, the PP04 system has a derivativediscontinuity when F = 0. This is described in [PP04] as a reflection of the non-linearity of the interactions between deep stratification, bottom water formationand thermohalide circulation . In the paper [PP04], Paillard et. al. considered the values given in a table3.3 to produce their figures. These parameter values were considered from firstprinciples and were obtained experimentally. We will use the same values forour analysis, except that following a discussion with Prof. Paillard at the July2017 CliMathNet Conference, we take γ = 0 . In [DSCW13, MA14, AD15] a Fourier series representation was given for theastronomical forcing at the Northern hemisphere summer solstice at 65 latitude.The resulting expression is given by I ( t ) = 1 e (cid:88) i =1 [ s i sin( ω i t ) + c i cos( ω i t )] . (5)Here the values of e , w i , s i and c i are given in [DSCW13, MA14, MCA15] andare found through through linear regression over the past one million years to8ariables values Range a . b . c .
01 0-0.15 d .
27 0.253-0.302 x . y . z . α .
15 0-0.35 β . γ . δ . τ V kyr ) 13.2-18.1 τ C kyr ) 3.1-15 τ A kyr ) 9.5-26 Table 1: The model parameter values used in the original PP04 model. In this paper we take γ = 0 . the present. The parameter e is a scale factor used to make the function forastronomical forcing dimensionless. Different researchers take different values of e for different models. For instance Ashwin and Dietlevsen [AD15] considered e = 1 for the AD15 model. However Mitsui and Aihara [MA14] considered e = 11 . W m − for the Crucifix-De Saedeleer model and e = 18 . W m − for the SM90,SM91 and PP04 models, basing the value on the three frequencycomponents of astronomical forcing that they considered to be significant. Incontrary Mitsui et. al. [MCA15] considered the parameter e = 23 . W m − forthe PP04 model. We will take the latter value for this paper.According to Mitsui et al [MA14], the three astronomical forcing components:the precession terms at i = 1 (23.7 kyr) and i = 3 (19.1 kyr) as well as obliquityterm at i = 4 (41.0 kyr) constitute 78 percent of original insolation forcing. Itis thus reasonable to consider the simplified astronomical forcing as a quasi-periodic function comprising three harmonics that includes precession at 19or 23 kyrs and the (dominant) obliquity forcing at 41 kyr ( ω = 0 . s ≈ −
11 reported in [MA14] and setting all other coefficients tozero, and dividing by e = 23 .
58, is approximately µ = 0 . . Accordingly, forthe remainder of this paper we will use the frequency and forcing amplitude ω = 0 . , and µ = 0 .
467 (6)as parameters of a physically realistic single mode insolation forcing [MA14,DSCW13]. However, to understand the general behaviour of the model, we willexplore the dynamics which results from taking other values of these parameters.9 . The PP04 model as a Filippov system.
Although the threshold models described above, and in particular the PP04model, are non-smooth in nature, they have been studied so far as smoothsystems, See for example [MA14, MCA15]. In this analysis only the typesof dynamics peculiar to smooth dynamical systems were observed, and Mitsuiremarked that this was a limitation of the smooth analysis. Indeed, hybrid dy-namical systems are canonical examples non-smooth systems and can be studiedbest by using the theory of non-smooth dynamical systems [BBCK08] in orderto find all the dynamics present in the system. This is the motivation for theapproach used in this paper.Non-smooth dynamical systems arise in a large number of applications and asmodels of a number of phenomena. They are used in mechanical engineeringin vibro-impacting systems, or in switches in electronic circuits such as ther-mostats and also in climate models [DBH10, GST08]. Discontinuous dynamicalsystems are systems where the vector field is piece-wise smooth (discontinuous).Therefore the dynamical system is known to be non-smooth as its trajectoriesmay not be differentiable everywhere. Non-smooth dynamical systems are char-acterized by some discontinuity in their right hand side, which can be due to thediscontinuities in evolution with respect to time, or the system state reachinga discontinuous boundary [DBBC +
08, CD10]. Consequently they can be usedto represent numerous physical processes which are characterized by periodsof smooth evolution being interrupted by an instantaneous event or the sys-tems whereby the physical states switches between two or more different states[BBCK08]. Therefore, when modelling such physical states, each state is givenby a different set of differential equations [AFO05]. That is, in each region,the evolution of trajectories are defined by the smooth dynamical system whichchanges to a different defining system across the discontinuity boundary [Gle16].Thus the behaviour is that of a piece-wise smooth dynamical system [BBCK08].
A Filippov system is a general piece-wise smooth dynamical system comprisinga finite set of ordinary differential equations, which can be expressed as˙ x = N i ( x ) x ∈ S i ⊂ R n . (7)Here each subspace or region S i has a non-empty interior, and the vector field N i is smooth and defined on the disjoint open regions S i . The intersection Σ ij of S i and S j is either an R n − dimensional manifold included in the boundariesof the two regions or it is an empty set. A non empty border between anytwo or more regions S i is called a discontinuity boundary or switching manifold[BBCK08, CDBHJ12]. A piece-wise smooth system with a single discontinuityboundary (such as the PP04 model) can be defined by:˙ x = (cid:40) N ( x ) if x ∈ S N ( x ) if x ∈ S . (8)10nd we call Σ ≡ Σ According to Cort´es [Cor08] and di Bernardo et al[BBCK08], the degree of smoothness of the piece-wise smooth system depends onwhether the system exposes jumps and or switches on its state, vector field or itsJacobian. The degree of smoothness at the point x on the discontinuity bound-ary set is given by the highest order r such that the Taylor expansions of the flowseither side of Σ (assumed to be at time t = 0) agree up to terms of O ( t r − ). Thisinforms us about the behaviour of the flow as it crosses the boundary[BBCK08].Systems have degree of smoothness one if F i ( x , µ ) − F j ( x , µ ) (cid:54) = 0 for x ∈ Σ ij ∩ D [BBCK08] and are called Filippov systems. According to Piiroinen et al [PK08],an important feature of a general Filippov systems is the possibility of motionto be constrained to the discontinuity boundary where the orbit can slide. Wewill show that this does not arise in the PP04 model, which is an importantaspect of its dynamical behaviour.
We can formulate the PP04 model as a forced Filippov System. (A similarformulation of an ice-line model for the glacial dynamics is given in [WWHM16]).To do this we introduce a state vector X = ( V, A, C ) T . According to Paillard et. al. [PP04], the inclusion of the I ( t ) term in thedefinition of the function F does not affect the times when the glacial cyclesterminates or the qualitative form of the overall dynamics. The proportionalitycoefficient c = 0 .
01 in their model is very small (in comparison to a, b and d ) andtheir range of values for c includes c = 0. Setting c equal to zero significantlysimplifies the theoretical analysis of the PP04 model, without changing theobserved dynamics in any significant way. Accordingly we set c = 0 for theremainder of this paper.With this simplification, it then follows from (1) that in all regions F ( X ) = ( a, − b, T X + d ≡ c T X + d. (9)The discontinuity surface Σ is then given by the linear relationΣ = { X : c T X + d = 0 . } (10)We define the following two states corresponding to the glacial and inter-glacialstates S ≡ S + = { X : F ( X ) > , S ≡ S − = { X : F ( X ) < . } . (11)The PP04 model in S ± can then be written as:˙ X = L X + b ± + I ( t ) e (12)11ere the linear operator L and the vector e are defined by L = − /τ V − x/τ V /τ A − /τ A − β/τ C − /τ C , e = − y/τ V α/τ C . (13)It follows from a direct calculation that the linear operator L has negativeeigenvalues − λ < − λ < − λ , with corresponding eigenvectors e i , i = 1 , , b ± depend upon which region X lies in and are givenby: b + = z/τ V δ/τ C , b − = z/τ V γ + δ ) /τ C . (14)It is clear from this formulation that the PP04 model has a piece-wise linearFilippov structure. We can thus expect it to have similar dynamics to a typicalFilippov problem and to show both ’smooth’ and ’discontinuity induced’ bifur-cations (for example grazing bifurcations [ ? Sim10]) as parameters are varied.Indeed this is exactly what we will see in this paper. In the paper [WWHM16]a Filippov system of a similar form to the above was analysed for the ice-linemodel, and some of the ideas used in studying that system can be applied inthe PP04 model.We now look at the structure of the Fillipov formulation of the PP04 model.
Lemma 4.1 (i) The solutions of the PP04 system remain bounded for all time.(ii) There is an attracting region B in the X -phase space into which all trajec-tories enter. Proof As L has all negative real eigenvalues, it can be written as L = U Λ U − where Λ = diag ( − λ , − λ , − λ ) . If we set Y = U − X , p ± = U − b ± and q = U − e then ˙ Y = Λ Y + p ± + q I ( t ) . Now consider N = Y T Y / N is sufficiently largethen ˙ N = Y T Λ Y + Y T (cid:0) p ± + q I (cid:1) < − min( λ i ) N + Y T (cid:0) p ± + q I (cid:1) < . Hence N , and thus | X | , is bounded. To prove (ii) we note (from inspectionof the actual matrix) that the matrix U T U is positive definite. It follows thatbounded sets in Y correspond to bounded sets in X and vice-versa. Hence the N − ball in the Y space corresponds to a bounded set B in the X space. (cid:3) Lemma 4.2
The degree of discontinuity of the PP04 model is one. roof It is clear from the formulation that X is continuous on Σ, but that ˙ X has a jump discontinuity. The result then follows. (cid:3) The following results describe the change of F across Σ and show that we donot have sliding solutions. Lemma 4.3 (i) F and dF/dt are continuous across Σ .(ii) At any point on Σ we have ( d F/dt ) + = ( d F/dt ) − + α where α > is apositive constant. Proof (i) If F ( X ) = c T X + d. The continuity of F is immediate. It also follows immediately that dFdt = c T ddt X = c T L X + c T b ± + c T e I ( t ) . (15)Then if we define h = c T L, r ± = c T b ± , g ( t ) = c T e I ( t ) , (16)we have dFdt = h T X + r ± + g ( t ) . (17)However, it follows directly from the definition of c in (9) and of b ± in (14) that r − = r + ≡ r. So ddt F = h T X + r + g ( t ) . (18)It is clear that dF/dt is then continuous across the discontinuity surface.Similarly we have d Fdt = h T ( L X + b ± ) + ˙ g ( t ) . (19)Thus [ ¨ F ] + − = h T ( b + − b − ) ≡ α = 0 . . (cid:3) If we approach Σ from S + it follows that dF/dt ≤
0. In particular if dF/dt < S − and does not slide on Σ.It is possible for grazing to occur on Σ. This arises when F = 0 and dF/dt = 0.In the case of an unforced system this will arise when c T X + d = 0 and h T X + r = 0 .
13t follows immediately that in this case grazing on Σ occurs along a straightline, the grazing set G , which is parallel to the vector c × h . We note further that in this case we have d Fdt = h . ( L X + b ± ) . Hence the surface d F/dt = 0 is another plane in each region S ± . This canintersect G at at most one point. This rules out the possibility of sliding.Following this result, we can, without ambiguity make the following definitions:Σ + = { X ∈ Σ : dF/dt > . } Σ − = { X ∈ Σ : dF/dt < . } (20)
5. The dynamics of the unforced PP04 model
We now study the unforced PP04 model which arises when there is zero inso-lation forcing, and consequently µ = 0. In this study we show that for certainparameter values this (non-smooth) model has periodic solutions, which ariseat border collision bifurcations between the fixed points and Σ as parametersin the model change. This form of the periodic solutions are similar to thatobserved in [WWHM16]. It is easy to see that the PP04 model has two fixed points given by Z ± = − L − b ± . (21)As L has negative eigenvalues, these are both attracting nodes. We define K ± = F ( Z ± ) = − c T L − b ± + d. (22)If K + > Z + lies in S + and is a physical fixed point. Any orbit whichremains in S + for all time will evolve towards it.If in contrast K + <
0, then Z + lies in S − , and is a virtual fixed point. Ithas a stable manifold in S + and attracts trajectories in S + towards it. Suchtrajectories ultimately cross Σ and enter S − . An exactly similar situation arisesfor the fixed point Z − . A border collision bifurcation (BCB) occurs when eitherof the two fixed points crosses Σ as a parameter varies.14 .2. The dynamics of the unforced system as parameters vary. We now establish the following result which describes the changing dynamics ofthe unforced system as parameters vary.
Theorem 5.1
Let Z ± be defined as above(i) If c T Z + + d ≡ d − L − b + > then Z + is a unique globally attracting fixedpoint.(ii) If c T Z − + d ≡ d − L − b − < then Z − is a unique globally attracting fixedpoint.(iii) If c T Z − + d ≡ d − L − b − > and if c . Z + + d ≡ d − L − b + < then thesystem has a periodic solution P ( t ) and no fixed points. NOTE We see from this lemma that the unforced system has either a fixed pointor a period orbit, but not at the same time. This is in contrast to the Saltzmanand Marsh models [SM91, SM91], but it is identical to the situation describedin [WWHM16] where the periodic orbit is called a ’flip-flop’ orbit.
Proof (i) Let X ∈ S + then provided that X ( t ) ∈ S + we have X ( t ) = e Lt ( x − Z + ) + Z + . (23)Hence, if X remains in S + for all time, then (as L has negative eigenvalues) itmust asymptotically tend towards Z + .Now, suppose that X ( t ) enters S − . In this region we have X ( t ) = e Lt ( x − Z − ) + Z − . (24)Hence it is attracted towards the fixed point Z − which lies within the region S + . Thus X must reenter the region S + at some later time. We claim that thetrajectory either remains in S + for all time following this, and converges to Z + ,or has a finite number of further ’visits’ to S − before remaining S + and thenconverging to Z + To establish this result we suppose first that the trajectory enters S + at time t ,leaves at time t , renters at time t etc. so that F ( t k ) = 0 , ˙ F ( t j ) > , ˙ F ( t j +1 <
0. The function F is defined by F = c T X + d , and hence in each regions S ± ithas the general form F ( t ) = a j e − λ ( t − t j ) + b j e − λ ( t − t j ) + c j e − λ ( t − t j ) + K ± , t j < t < t j +1 , (25)where, in this case, K + > K − > . We firstly establish the following
Lemma 5.2 If c j > then the trajectory remains in S + for all t > t j . Proof . Suppose the converse. There must be a later time t j +1 for which F ( t j +1 ) = 0. Now consider the globally defined function F ∗ ( t ) = a j e − λ ( t − t j + b j e − λ ( t − t j ) + c j e − λ ( t − t j ) + K + . c j > F ∗ ( t ) tends to K + > t , but F ∗ ( t j +1 ) = 0. By considering the shape of the curve F ∗ we deduce that thereare times t j < t a < t j +1 < t b < t c so that F ∗ ( t a ) > , F ∗ ( t b ) < , F ∗ ( t c ) >K + > F ∗ ( t a ) = ˙ F ∗ ( t b ) = ˙ F ∗ ( t c ) = 0 . However, it is immediate that˙ F ∗ ( t ) is a sum of three different exponential functions. It is well known that afunction which is the sum of n different exponential functions can have at most( n −
1) zeros. Thus we have a contradiction. (cid:3)
Now consider the case of c j < S − at a time t j +1 and then back into S + at a time t j +2 . We consider the map G ( c j ) → c j +2 . Lemma 5.3 G ( z ) = α j z + β j where < α j < and < β j < . Proof
At the time t j +1 we have (from Lemma 4.3) that F ( t − j +1 ) = F ( t +2 j +1 ),˙ F ( t − j +1 ) = ˙ F ( t +2 j +1 ) and ¨ F ( t − j +1 ) = ¨ F ( t +2 j +1 ) + α . It follows, after some ma-nipulation, that the coefficients a j , a j +1 etc. obey the linear Vandermondeequation λ λ λ λ λ λ e − λ ∆ j a j − a j +1 e − λ ∆ j b j − b j +1 e − λ ∆ j c j − c j +1 = K − − K + α , (26)where ∆ j = t j +1 − t j . From the data given, K − − K − = c T ( Z + − Z − ) = 1 .
04 and α = 0 . . We deduce, on inverting the Vandermonde matrix, that e − λ ∆ j a j e − λ ∆ j b j e − λ ∆ j c j = a j +1 b j +1 c j +1 + . − . . . (27)Applying the same result at the time t j +2 with ∆ j +1 = t j +2 − t j +1 we have c j +2 = e − λ ∆ j +1 (cid:0) e − λ ∆ j c j − . (cid:1) + 1 . . (28)The form of G given in the Lemma follows immediately.Now suppose that the trajectory always re enters S − . It follows from Lemma5.2 that c j must always be negative. However, from Lemma 5.3 we have that c j +2 − c j = ( α j − c j + β j > . Thus the sequence c j bounded above (by zero) and is monotone increasing.It must therefore tend to a limit c , which in the limit satisfies c = α j c + β j .As α j and β j are always positive, this is a contradiction. We deduce that c j is eventually positive, at which point the trajectory remains in S + and hencetends to the fixed point Z + . This proves part (i).16he proof of (ii) is identical that that given above.To prove part (iii) we use the following argument, which is illustrated in Figure3. Figure 3: Schematic showing the period orbit P ( t ) and the virtual fixed points Z ± . Consider a trajectory which starts at the point X on Σ + and initially enters S + so that dF/dt >
0. As we are in case (iii), it follows that F = a + e − λ t + b + e − λ t + c + e − λ t + K + , (29)where K + < . For large t >
0, we must have
F <
0. Then there must be afirst time t = t at which F = 0 and dF/dt <
0. At this point the trajectoryintersects Σ − at the point X = X . The flow now crosses over into S − with F <
0. The resulting flow is then given by F = a − e − λ t + b − e − λ t + c − e − λ t + K − with K − >
0. By the same argument, it follows that there is a first time t > t such that F = 0 and the trajectory intersects Σ + at the point X . The conditionfor this trajectory to be a periodic solution is that X = X ≡ M ( X ) (30)This can be considered to be a fixed point condition for the nonlinear map M : Σ + → Σ + defined above.It follows immediately from Lemma 4.1 that the function M maps the finitedimensional and bounded region B ∩ Σ + into itself.17e now show that M is continuous. The trajectories in S + depend smoothlyupon the initial value X , hence F ( t ) is a smooth function of X . It thereforefollows from the implicit function theorem that, provided dF/dt < X , thenthe time t , and the point X , are continuous (indeed differentiable) functionsof X . Similarly, the point X will also be a continuous function of X . Thecontinuity of the map M then follows provided that ˙ F ( X ) < . To prove thiswe establish a contradiction. Consider the function given (29). Suppose thatat times t and t we have F = 0 and that also ˙ F ( t ) = 0 so that F ( t ) > t is close to t and t > t . As K + < t such that F ( t ) = 0 . It follows from Roll´e’s Theorem that there must be times t a and t b with t < t a < t < t b < t such that˙ F ( t a ) = ˙ F ( t ) = ˙ F ( t b ) = 0 . Now, as before, ˙ F is a sum of three exponential functions. Such a functioncannot have three zeros. Thus we have established the desired contradiction.We have thus established that the function M maps a bounded finite dimen-sional region into itself, and is continuous. The existence of a fixed point, andhence of a periodic orbit, then follows immediately from the Brouwer FixedPoint Theorem. (cid:3) NOTE A similar system was studied in [WWHM16] by using a contraction map-ping argument which could be applied directly to their problem and using whichthey could also prove uniqueness of their periodic orbit for certain parametervalues.If we take the tabulated values for the PP04 model, with γ = 0 . d = 0 . Z + = (0 . , . ,
0) and Z − = ( − . , − . ,
2) At thesepoints we have K + = c T Z + + d = − .
05 and K − = c T Z − + d = 0 .
99 so that thecondition for a periodic solution is satisfied. The time series of the componentsof the resulting periodic solution (which appears from these calculations to beunique), is then illustrated in Figure 4. These show a saw-tooth like structuressimilar to those evident in the geological reconstructed data. Such an oscillationwas observed in the original PP04 model (see [Cru12]). The relaxation oscillatorobtained for the parameters we use has a period of about 147 kyr . Note that theperiod of this unforced oscillation is higher than the observed period of 100 kyr ,but is not dissimilar. This model therefore suggests that the natural timescalesof the Earth do play a role in determining the frequency of the ice ages.
If we vary one of the parameters of the system, say d , then the periodic solutioncan lose existence at a border collision bifurcation (BCB), when either one of18 igure 4: Unforced periodic solution showing V (red), A (blue) and C (magenta), as well as F (black). the two (virtual) fixed points Z ± intersects Σ. We then see a change from aperiodic solution to a fixed point.Qualitatively, the behaviour close to the BCB is illustrated by the representativephase-plane diagram in Figure 3 given earlier. In this we show the periodicsolution when the two fixed points are virtual. The solid lines show the truedynamics in S + and S − , and the dotted lines the ’virtual’ dynamics if, forexample the dynamics in S + is extended into S − so that it approaches thevirtual fixed point. Even if the fixed point is close to Σ this periodic orbit has anon-vanishing amplitude, indeed the amplitude tends to a non-zero limit as oneof the fixed points, say Z + approaches Σ. However, as the BCB is approachedthe period of the periodic solution increases as it takes longer to approach Σ.The period rises to infinity when the fixed point Z + lies on Σ.The values of d ± at which we have a BCB occur when either K + = 0 or K − = 0.These cases arise when d ± = − c T Z ± (31)For the tabulated values we obtain d − = − . , d + = 0 . , and hence a periodic solution exists when − . < d < .
32. In Figure 5 weshow the period of the periodic solution as a function of d . In which we can seethe two BCBs at which the period tends to infinity.It is of interest, both theoretically, and also from the need to do computations,to consider how this bifurcation structure arises if we replace the non-smoothsystem by a smooth one. A convenient way to do this (see for example [Cru12])19 igure 5: The change in the period of the periodic solution as the parameter d is varied. Inthis we can see the Border Collision Bifurcations at d = − .
72 and d = 0 .
32 where the periodbecomes infinite. is to replace the (non-smooth) Heaviside function, by the regularized function H η ( z ) = 12 (1 + tanh( ηz )) . (32)For large values of η this closely approximates the Heaviside function. Usingthis approximation we integrate the system (2-4) forward in time numericallyby using the Matlab stiff ode solver ode15s . To determine the dynamics ofthe solution we then start with a random set of initial conditions and find thesolution of the dynamical system starting from these. We then take a largeenough time interval to allow the solution to converge onto its Ω − limit set. Torecord this set we then plot the maximum and minimum values of F on theasymptotic orbit. We choose to plot F as this then allows us to see how theOmega limit set of the solution interacts with the discontinuity surface. Bydoing this for a set of values of d we can determine the complete bifurcationpicture for the solutions.If η = 1600 then H η ( z ) is a very good approximation to the Heavyside function,and we expect the dynamics of the smoothed system to be very close to thatof the Filippov system describing the PP04 model. The numerically computedbifurcation picture of the asymptotic behaviour of the solution as a function of d is given in Figure 6.In this figure, as d increases, we see the value of F ( Z − ) increasing linearly with d until the BCB when F ( Z − ) = 0. The fixed point is then immediately replacedby a periodic orbit with non-zero amplitude, which is in turn destroyed at thesecond BCB when F ( Z + ) = 0.In a second figure we consider the bifurcation diagram, close to the rightmostbifurcation point, of the solutions as a function of d , when η = 400 ,
800 and20 igure 6: The bifurcations of the solutions when η = 1600 showing the two border-collisionbifurcations at d = − .
72 and d = 0 . Figure 7: The bifurcation diagram of the solutions close to d = 0 .
32 when (from left to right) η = 400 , , In this figure we see that in all cases the fixed point Z + loses stability, as d isdecreased. The stability is lost to a periodic solution in what appears to be asuper-critical Hopf bifurcation at a value of d close to, but slightly smaller than,the BCB value of d = 0 .
32. Initially the periodic orbit is close to the fixed point.However, as d is decreased further there appears to be a cyclic fold bifurcationat which point the periodic orbit expands rapidly in size, to approach the orbitto the discontinuous system. This behaviour was also observed in [Cru12]. The21udden increase in the size of the periodic orbit as d is decreased for the caseof η = 400 is shown in Figure 8 in which we plot the trajectories in the ( F, V )phase plane for d = 0 . , d = 0 . , d = 0 . , d = 0 . . Figure 8: The (
F, V ) phase plane of the solution as d is increased from 0.317 (black) to 0.318(red) showing the rapid decrease in the size of the periodic orbit. We see from these calculations that as η increases the dynamics of the smoothsystem rapidly approximates the dynamics predicted by the analysis of the non-smooth system, with the (relatively simple) border-collision bifurcation in thenon-smooth limit being replaced by a nearby, and more complex bifurcationstructure in the smooth system. In the next section we will look at how the(apparently unique) periodic solution derived above changes when a periodicinsolation forcing term is added.
6. The analytic dynamics of the periodically forced PP04 model.
We now consider the solutions of the PP04 model when the insolation forcing hasa single periodic mode. Clearly this form of forcing is unrealistic from a physicalpoint of view. However, studying such systems allows us to gain insight intothe more general case of quasi-periodic forcing, especially when one frequencyis dominant in the insolation forcing. Indeed we will give evidence at the end ofthis paper that the behaviour of the quasi-periodically forced system is a simpleperturbation of the periodically forced case.Accordingly, in this section we suppose that the insolation forcing has the form I ( t ) = µ sin( ωt ) , (33)22o that the period of the forcing is given by T = 2 πω . (34)In general we might expect to see the following types of solution behaviour forthe PP04 system:(a) Synchronised periodic solutions (both stable and unstable) of period P = nT = 2 πn/ω with n = 1 , , .. which have precisely one glacial and oneinter-glacial period (one glacial cycle) between repeats. We define theseto be (1 , n ) periodic orbits.(b) Synchronised periodic solutions with several (for example m ) differentglacial cycles between repeats, of ’average’ period P = nT /m with m =1 , , , .. . We define these to be ( m, n ) periodic orbits.(c) Quasi-periodic solutions showing at least two distinct frequencies.(d) Chaotic solutions.In practice, for appropriate choices of parameters, we see all of these types ofsolutions, possibly co-existing. Some of these solutions arise through smoothbifurcations and others (as we have seen in the previous section) from non-smooth bifurcations as we vary parameters such as µ and ω .In this section we will consider those ( m.n ) solutions which are initially smallperturbations of the periodic orbit of the unforced system constructed in thelast section, given by taking the insolation forcing amplitude µ to be small.The underlying periodic orbit has a well defined frequency ω ∗ and which (as wehave shown) intersects the surface Σ transversely. It follows from [BBCK08] thatclose to this orbit, the Poincar´e return map defined by P S : X ( t ) → X ( t +2 π/ω )is smooth. Thus we may apply the theory of Arnold Tongues [PRKK03, Sim10]to predict the existence of ’tongues’ which are curves of (say) ( µ, ω ) which definethe boundaries of the existence regions for synchronised periodic solutions forthe ( m, n ) orbits when ω ≈ nω ∗ /m . Such tongues will be expected to have | ω − nω ∗ /m | varying proportionally to µ m . (1 , n ) periodic solu-tions. It is relatively easy to construct algebraic conditions the satisfaction of which isnecessary for the existence of the (1 , n ) periodic orbits. Suppose that we have aperiodic solution X ( t ) of period P = 2 nπ/ω , and for which F ( X ( t i )) = 0. Forsuch a periodic orbit we will assume that exactly one glacial cycle exists for t in the range t ∈ [ t , t = t + P ] . For this cycle we assume that the solution is glacial if t ∈ [ t , t ], with F ( t ) > X ( t ) ≡ X + ( t ). Similarly the solutionwill be inter-glacial if t ∈ [ t , t ], with F ( t ) < X ( t ) ≡ X − ( t ). We define23he set of points X i = X ( t i ). It then follows that a periodic orbit must satisfythe conditions: X = X and F ( t i ) ≡ c T X i + d = 0 i = 0 , . (35)The differential equations satisfied by this system in the two regions are thengiven by: ˙ X ± = L X ± + b ± + µ e sin( ωt ) . (36)A particular integral of this system is given by X ± P I ( t ) = Z ± + µ p cos( ωt ) + µ q sin( ωt ) ≡ Z ± + µ r ( t ) . (37)Where (as before) Z ± = − L − b ± , p = − ( L + ω I ) − ω e , q = − ( L + ω I ) − L e . (38)We can then integrate the whole system to give X = e L ∆ (cid:0) X − X + P I ( t ) (cid:1) + X + P I ( t ) , (39)and similarly X = e L ∆ (cid:0) X − X − P I ( t ) (cid:1) + X − P I ( t ) . (40)Here we set ∆ = t − t , ∆ = t − t = P − ∆ . (41)For a given period P , the problem of existence of a periodic solution is thenreduced to finding the 5 unknowns comprising X , together with the initialphase t and the time of the transition between the glacial and inter-glacialstates at t , so that the five equations in (35) hold. (Alternatively we can take t as given and find the period P as part of the solution.) This nonlinear systemmay or may not have algebraic solutions, and we will consider this in the nextsub-section. Furthermore the algebraic solutions, if they exist may or may notlead to physically relevant climate trajectories X ( t ), defined as follows: Definition
We define a (1 , n ) periodic solution to be physical if V ( t ) > , A ( t ) > , C ( t ) > t ∈ [ t , t ] (42)and F ( X + ( t )) > , t < t < t , and F ( X − ( t )) < t < t < t . (43)Typically solutions lose algebraic existence through smooth (saddle-node orperiod-doubling) bifurcations, and lose physicality through non-smooth (graz-ing) bifurcations, where we expect to see a dramatic change in the solution asindicated in Chapter 7 of [BBCK08]. We will return to this situation later.24 .3. Small µ , synchronised, (1 , n ) periodic solutions. We consider first the question of the existence of the (1 , n ) periodic solutions.To do this we use perturbation theory, and look for explicit representationsof the periodic solutions of the forced system, which are perturbations of theperiodic solution of the unforced system when µ is small . We have shown inthe previous section that if µ = 0 (the unforced system) there is a periodicsolution with one glacial cycle of frequency ω ∗ and period P ∗ = 2 π/ω ∗ andwhich takes values X ∗ = ( V ∗ , A ∗ , C ∗ ) at the start of the glacial cycle. Extensivenumerical experiments strongly indicate that this solution is also unique (upto the arbitrary starting time t ) and attracting. Accordingly, for forcing atfrequency ω with small µ we might expect to see a synchronised (1 , n ) periodicorbit of period P provided that P ≈ P ∗ so that P = 2 πnω ≈ P ∗ = 2 πω ∗ . For the remainder of this section we will consider the periodic orbit that arisesfor the tabulated values of the parameters (so that for example d = 0 .
27 forwhich we have found that ω ∗ = 0 . . It follows that ω ≈ n ω ∗ . (44)We propose that for small µ that there is a range of ω values with | ω − n ω ∗ | = O ( µ ) such that two synchronised periodic solutions exist within this range. Overthe interval the phase t of each such solution (defined as the phase of the forcingat the start of the glacial cycle) is well defined and varies over the whole range[0 , π/ω ∗ ]. Both solutions are perturbations of the unforced solution which isgiven when µ = 0, and which has an arbitrary phase. Hence, both solutions arephysical provided that µ is sufficiently small and the parameter d is not close tothe value at which a border collision occurs for the free system. The boundariesof the regions of existence of both solutions are determined by the existence ofsaddle-node bifurcations where the solutions coalesce.This result is an immediate consequence of the following: Lemma 6.1 (i)
For each n , if µ is small then there is a set of solutions ( ω, x ) to the algebraic system, which is parameterised by t . (ii) if µ is small then the curves ( ω, V ( ω )) of the (1 , n ) orbits form ellipses whichhave the point ( nω ∗ , V ∗ ) at the centre. (iii) The size (for example the semi-major axis) of the ellipses is (for sufficientlysmall µ ) directly proportional to µ . (iv) As we go once around the small elliptical curves, the phase t increases bya factor of π/ ( nω ∗ ) . orollary 6.2 If µ (cid:28) then for each value of n , the synchronised periodicsolutions of the PP04 model exhibit saddle node (SN) bifurcations at points ( ω ,n , ω ,n ) . At which points they change to quasi-periodic orbits. Synchronisedperiodic solutions exist in the interval ω ∈ ( ω ,n , ω ,n ) . We have that | ω i,n − nω ∗ | = O ( µ ) , i = 1 , . We illustrate the conclusions of Lemma 6.1 and Corollary 6.2 in Figure 9, inwhich we plot the ellipses corresponding to the solutions ( ω, V ( t )) for the pa-rameter values µ = 0 . , . , . n = 3. In Figure 10 we plot ω as a function Figure 9: The variation of V ( t ) with ω for µ = 0 .
05 (red), µ = 0 . µ = 0 . n = 3 showing the (elliptical) curve of the solutions and the two saddle-node bifurcationpoints. It is clear that the size of the closed curve increases in proportion to µ , and that ithas a true elliptical shape for the smaller values of µ . of t for the case of the (1 ,
3) orbits which arise when µ = 0 . n = 3 as wego around the ellipse. For the problem considered we have nω ∗ = 0 .
128 Wecan see that t increases by approximately 2 π/ ( nω ∗ ) = 49 . ω ( t ) , V ( t ) , C ( t )) together for µ in-creasing from 0 .
01 to 0 .
24. We see that the closed elliptic curves exist for theseparameter values in this extended space.
Proof If µ = 0 we have a periodic solution X ∗ ( t ) of the autonomous system.This can have an arbitrary time t ∗ at the start of the glacial cycle for which F ( X ∗ ( t ∗ )) = 0. We have a well defined set of time differences ∆ ∗ and ∆ ∗ sothat the glacial period is in the interval [ t ∗ , t ∗ + ∆ ∗ ] and the inter-glacial periodin the time interval [ t ∗ + ∆ ∗ , t ∗ + ∆ ∗ + ∆ ∗ ] with ∆ ∗ + ∆ ∗ = 2 π/ω ∗ . In the forcedcase, with µ >
0, the system will exhibit phase-locking , so that we expect to seea well defined start time t in this case. Finding the form of the solution and thefrequency ω , in terms of µ and t will be part of the solution process. Suppose26 igure 10: The variation of ω with t for µ = 0 . n = 3) showing that t increases byapproximately 2 π/ (3 ω ∗ ) = 49 . V ( t ) and C ( t ) with ω when n = 3 with µ increasing from 0.01to 0.24. X ( t ) of period T = 2 πn/ω ≈ π/ω ∗ , so that ω ≈ nω ∗ , with a glacial state in the interval t ∈ ( t , t ) and an inter-glacial statefor t ∈ ( t , t ) so that t = t + 2 πnω . We define X i = X ( t i ). Thus to have a periodic solution we must satisfy thefollowing three conditions X = X and F ( X i ) = c T X i + d = 0 . (45)We now consider a solution which is a perturbation of the unforced case so thatto order O ( µ ) we have ω = nω ∗ + µα, ∆ = ∆ ∗ + µδ, X = X ∗ + µ x . To find the leading order form of the perturbed solution we then determine α, δ, x as functions of the phase t .It follows immediately that∆ = 2 πnω − ∆ = ∆ ∗ − µδ − µ πα ( ω ∗ ) . It then follows from (37) that X = e L (∆ + µδ ) ( X ∗ + µ x − Z + + µ r ( t )) + Z + + µ r ( t ) . Hence, after some manipulation X = X ∗ + µ (cid:16) e L ∆ ∗ ( δL X ∗ + x − r ( t )) + r ( t ) (cid:17) + O ( µ ) . Thus X = X ∗ + µ y + O ( µ ) , where y = e L ∆ ∗ ( δL X ∗ + x − r ( t )) + r ( t )Similarly, X = X ∗ + µ z + O ( µ ) , where z = e L ∆ ∗ (cid:0) − ( δ − πα/ ( ω ∗ ) )) L X ∗ + y − r ( t ) (cid:1) + r ( t ) . The conditions F ( X i ) = 0 are then satisfied provided that x = z and c T x = 0 , c T y = 0 .
28e define the linear operators A and A by: e L ∆ = A , e L ∆ = A . It follows that y = A ( δLX ∗ + x − r ( t )) + r ( t ) , z = A ( − ( δ + γα ) LX + y − r ( t )) + r ( t ) . Hence, as z = x we have x = A ( − ( δ + γα ) L X + A ( δLX + x − r ( t )) + r ( t ) , (46) c T x = 0 , c T y = 0 . (47)Now we look at the structure of the equations (46,47). We note that to leadingorder, as t = t + ∆ ∗ that there are vectors p , q , p , q so that r ( t ) = p cos( n ω ∗ t )+ q sin( n ω ∗ t ) , r ( t ) = p cos( n ω ∗ t )+ q sin( n ω ∗ t ) . It follows that there is a linear operator M and vectors a and b so that (46,47)can be put into the form M x δα + a cos( n ω ∗ t ) + b sin( n ω ∗ t ) = . (48)The linear operator M and the vectors a , b can all be constructed explicitly.We will make the assumption that M is invertible. Numerical evidence clearlyindicates that this is always the case. Under this assumption, for each value of t the system (48) can be solved uniquely to give the values of x , δ and α . Thesethen take the form x δα = f cos( n ω ∗ t ) + g sin( n ω ∗ t )for appropriate (constant) vectors f and g . In particular there will be uniquevalues f and g so that α = f cos( n ω ∗ t ) + g sin( n ω ∗ t ) . As t varies over the whole range of [0 , π/ ( n ω ∗ )] so α will range over theinterval. α ∈ (cid:20) − (cid:113) f + g , (cid:113) f + g (cid:21) . (49)This interval sets the limits of existence of the solutions of (46,47) and hencethe width of the tongues over which we will see synchronised periodic solutions.Clearly if W α = (cid:112) f + g then there is a phase φ α so that α = W α cos( n ω ∗ t − φ α ) . (50)29f we set [ vac ] = x T , then an identical argument implies that there are ampli-tudes W V , W A and W C , and phases φ V , φ A and φ C so that v = W V cos( n ω ∗ t − φ V ) , a = W A cos( n ω ∗ t − φ A ) , c = W C cos( ω ∗ t − φ c ) . (51)It follows immediately that the curves ( α, v ) , ( α, a ) and ( α, c ) are all ellipsescentred on the origin. (cid:3) µ solution ellipses. The values of the coefficients of the vectors f and g are determined explicitlyby the calculation above, but are hard to estimate from this. However, thebasic calculation of the (1 , n ) periodic orbits is identical for all values of n =1 , , , , .. although the precise values of the coefficients will change in each case.In particular, for small µ we expect to see small ellipses in each case, the size ofwhich is directly proportional to µ . In Figure 12 we plot the resulting ellipseswhen µ = 0 . n = 1 , , ,
4. These ellipses are computed by numericallysolving the nonlinear equations for
V, A, C, ω and t . As these solutions areparameterised by the initial time t , it is convenient in this calculation to use t as the path following variable. Each of these ellipses are centred on thevalues of ω = 0 . ω = 0 . ω = 0 .
128 and ω = 0 . n increases the size of the minoraxis of the ellipse appears to decrease, although the size of the major axis staysapproximately constant. The value of t varies over the interval [0 , π/ ( n ω ∗ )] Figure 12: The variation of V ( t ) with ω for µ = 0 . n = 2 , , as we travel around the ellipse. In particular, it follows from (49) that the valueof t changes by π/ ( n ω ∗ ) between the two saddle node bifurcation points. Thisis of interest as it demonstrates that the phase of the response ( V, A, C ) to theinsolation forcing, whilst locked to it for a particular periodic orbit, differs from30t. This phenomenon has been observed in the record of the ice ages, in whichthe Milankovitch cycles are not always seen to be in phase with the cooling andwarming periods.As an example we take the case n = 3 and µ = 0 .
1. A numerical calcula-tion in this case shows that solutions exist for ω ∈ [0 . , . π/ω ∈ [44 . , . , T = 6 π/ω = [133 . , .
57] and t ∈ [8 . , . . In Fig-ure 13 we plot two cycles of the resulting periodic orbits for the three cases t = 8 . t = 20 . t = 31 .
882 representing the left and right limitsand the middle of the range of values for which we see a solution. (1 , n ) orbits for small µ . The previous analysis has shown that ω ∗ = 0 . . Further numerical studieslead to the following approximations for small µ of the regions of existence ofthe (1 , n ) orbits. ω , = ω ∗ − . µ ω , = ω ∗ + 0 . µω , = 2 ω ∗ − . µ ω , = 2 ω ∗ + 0 . µω , = 3 ω ∗ − . µ ω , = 3 ω ∗ + 0 . µω , = 4 ω ∗ − . µ ω , = 4 ω ∗ + 0 . µ (52)In Figure 14 we give the graph of the regions of existence of the periodic solutionsfor n = 1 , , , µ > .
15. For µ > .
15 we will expect to see (as we in fact do see) the co-existence of periodicsolutions with different values of n and hence of different periods T = 2 nπ/ω . Infact, as we shall see, the original (nonlinear) problem has rather larger regionsof overlap of the existence regions. µ . The above calculation has given only a small µ analysis, showing that for small µ the width of the existence tongues and the associated ellipses of the solutionsof the algebraic system (35) increases in direct proportion to µ for all values of n . Similar results for other systems are given in [SP07].For larger values of µ nonlinear effects become important, and the ellipses de-termined above will form part of the complex surface of the solutions of (35). Inthis scenario, as we shall see, the ellipses calculated above become distorted, andthen can break up and expand as they coalesce with other curves of solutions.However, we note that (unlike the small µ case) many of the solutions of thealgebraic equations (35) for larger values of µ will not represent physical climatestates. For example this may be a trajectory starting from an initial state at t X calculated as a solution of (35) on the assumption that it remains in S + for t < t < t which may, in fact, cross Σ at a time t < t ∗ < t .In Figure 15 we plot computed regions of existence of the (1 , n ) periodic solu-tions. These regions are determined by first fixing the value of µ and solving the31 igure 13: Two cycles of the periodic solutions when n = 3. In this plot we see V (blue), A (red), C (maroon), F (black), and the insolation (green). We have Left: t = 8 .
79, Right: t = 31 .
882 and Bottom: t = 20 . igure 14: The existence regions for the periodic solutions of the linearised problem. full algebraic system numerically for a set of values of ω = nω ∗ ± δ increasing δ from 0. The calculation was done using the Matlab solver fsolve with an initialguess given by X = (5 , , , , , , , ω new against µ for which the algebraic solver breaks down. As can be seen, theregions of existence are linear (as predicted) for small values of µ . They thenexpand significantly as µ increases. This is due to a coaelescence of the small µ ellipses with other solution curves as described above.We note that the physically interesting case of ( µ, ω ) = (0 . , . ,
3) periodic solution, and we will returnto this observation later.In Figure 16 we see the set of elliptical curves for the cases of n = 2 , , µ than before. For n = 2 , µ = 0 .
25. In the case of n = 4 the coalesecenceoccurs for a larger value of µ . Indeed, we observe in general, that the coalescenceof the ellipses with other solution curves occurs for smaller values of µ as n decreases. We note further that if µ = 0 .
25 then (as expected from the linearanalysis) the regions of existence of the n = 2 and n = 3 periodic orbits overlap.As a consequence we might expect to see both n = 2 and n = 3 orbits in thiscase, with related domains of attraction for the initial data.If we take the larger, and physically relevant, value of µ = 0 .
467 then we seea more complicated curve, and the range of existence of the solutions in thiscase is more difficult to predict. In Figure 17 we show the curves of the (1 , µ increasing from µ = 0 . µ = 0 . µ < . µ = 0 .
25 and enlarge as µ increases. When µ = 0 .
467 we seethat the maximum value of ω = 0 . igure 15: A set of graphs showing the regions of existence of the (1 , n ) orbits when both µ and ω are varied. The red case n = 2 , , µ = 0 .
467 and consider the physically relevant value of ω = 0 . ω it is apparent from Figure 17 that there are (at least) two solutions, S , ≡ [ V ( t ) , A ( t ) , C ( t ) , t , ∆] to the algebraic equations, with S on theupper side of the curve of solutions and S on the lower. A careful calculationshows that these solutions are given by S = [0 . , . , . , . , . , and S = [0 . , . , . , . , . . The corresponding functions ( V ( t ) , A ( t ) , C ( t ) , F ( t )) are plotted in Figure 18,along with the insolation forcing. It is clear from this figure that only thesolution S can be physical. This is because when we consider the solution S we can see from the graph that the function F ( t ) does not keep a constant signduring either the glacial or the inter-glacial cycles.The solutions close to µ = 0 .
25 are of theoretical interest as here we see thereason for the break up of the closed elliptical curves. In Figure 19 we showthe solution existence curves for µ = 0 .
244 (left) and for µ = 0 .
245 (right). Thecurve for µ = 0 .
244 shows two separated solution branches, one of which is adistorted ellipse. As predicted earlier, these two branches then coalesce closeto µ = 0 .
25, leading to a sudden expansion of the rightmost elliptical curve. A34 igure 16: The variation of V ( t ) with ω when n = 2 , , µ = 0 . µ = 0 . µ = 0 .
25 (red) showing the break up of the ellipticalcurves when they coalesce with other solution curves as µ increases. Here n = 2 , , V ( t ) with ω for the (1 ,
3) orbit when µ increases from µ = 0 . µ = 0 .
467 showing the break up of the elliptical curve at µ = 0 . igure 18: The solutions on the upper branch ( S top) and the lower branch ( S bottom).In these graphs we plot (as functions of time) V (red), A (blue), C (purple), F (black) andthe insolation forcing (green). The solution S is not physical as F ( t ) changes sign within theglacial cycle.Figure 19: The variation of V ( t ) with ω for the (1 ,
3) orbit for µ = 0 .
244 (left) and µ = 0 . µ = 0 .
244 leading to anexpansion of the solution ellipse when µ = 0 . . t , ω ) and of ( t , V ( t )) for the case of µ = 0 .
25 is given inFigure 20. We see that unlike the case of small µ when t could take arbitraryvalues, in this case we have an upper limit of t <
78. We note, however, thatthe solutions on these curves are not necessarily physical as V Figure 20: The variation of V ( t ) (red), and of ω (blue) with t for the (1 ,
3) orbit when µ = 0 . . We can see evidence for a limit point at t = 0 . . As we have seen, not all of the orbits on the computed curves are physical, inthe sense that the function F on a solution trajectory can change sign at anintermediate point t < t ∗ < t during a glacial period, or similarly during aninter-glacial period.Also of significant interest is the stability of the resulting orbits. The rightextremes of the ( ω, V ) solution curves are in all cases marked by saddle-nodebifurcations. In general such bifurcations are associated with changes in the stability of the solutions. It is difficult to determine the stability algebraically.However a large number of numerical experiments demonstrate clearly that it isthe lower branch of the curves which is (in general) stable, and the upper branch is unstable.We will see later that as a parameter such as ω is varied, the solutions can alsolose stability at period-doubling bifurcations , where a ( m, n ) orbit is replacedby a (2 m, n ) orbit. A further loss of stability is associated with a grazingbifurcation , which is the first value of the parameter at which a solution losesphysicality with the trajectory grazing the discontinuity surface Σ. (Such eventsare known to be highly destabilising [BBCK08].) ( m, n ) periodic orbits A similar analysis can be applied to the more general ( m, n ) orbits. In suchorbits we see m glacial cycles of warming and cooling, in a period of 2 πn/ω .To construct, and analyse these, we introduce a series of m intervals ∆ i, and∆ i, with i = 1 . . . m −
1, summing in total to 2 πn/ω , being the times between37uccessive glacial and inter-glacial periods. Each such interval will start at atime t i, or t i, , with i = 0 . . . m −
1. Here each such t i, , can be computedfrom the initial time t i, of the first glacial cycle by adding up the appropriatetime periods ∆ i, , . For small µ Each ∆ i, , and ∆ i, is then a perturbation, δ i, or δ i, , of the respective times of the glacial and inter-glacial periods of theperiodic solution of the unforced problem. Similarly, we let X i, and X i, be theinitial conditions at the start of the respective glacial and inter-glacial periods.For small µ these will be perturbations x i, and x i, of the related values for theperiodic orbit of the unperturbed system. The algebraic equations for a (1 , n )orbit then extend to the following system for i = 0 . . . m − X i, = E ( t i, , ∆ i, , X i, ) , (53) X i +1 , = E ( t i, , ∆ i, , X i, ) (54) F ( X i, ) = 0 , (55) F ( X i, ) = 0 , (56) m − (cid:88) i =0 ∆ i, + ∆ i, = 2 πnω , (57) X , = X m, . (58)Here E ( t i, , ∆ i, , X i, ) is the evolutionary operator which we have constructedexplicitly. If we specify the start time t , and the amplitude µ then the system(58) constitutes 8 m + 1 equations for the 8 m + 1 unknowns X i, , X i, , ∆ i, , ∆ i, , and ω. As before, the complete system (58) can be linearised about the periodic solu-tion when µ = 0. In this case we take ω = nω ∗ /m + δω with | δω | (cid:28) . Theresulting system will be identical in form to that given in equation (48) with acorresponding linear operator M in this case. However, from the earlier discus-sion of the general rules for the asymptotic behaviour of the Arnold tongues, weexpect that δω = O ( µ m ) in this case.We present in Figure 21 an example calculation of solving this algebraic systemnumerically for the case of a periodic solution with ( m, n ) = (2 ,
5) with µ =0 . , .
02 and 0 .
05. In this figure on the left we plot V ( t ) as a function of ω ,and on the right we plot the period of the first full glacial cycle P G = ∆ , +∆ , as a function of V ( t ). It is clear from these figures that V ( t ) and P G have asingle ’cycle’ as t varies over one period, and the perturbation form the unforcedvalue scale linearly with µ , with the ( V ( t ) , P g ) curve being a perturbed ellipse.In contrast δω scales quadratically with µ and has a double cycle (in the formof a figure of eight) in this period. (A plot of the same curves for the (3 ,
5) orbitshows, as expected, similar behaviour for V ( t ) and P G and a triple cycle for δω which scales as O ( µ ) . An excellent account of the computation of Arnold tongues for general circlemaps is given in [SP07], with general surfaces for the solutions obtained forvarying parameters. The surfaces determined above (for example the ellipses38nd figure of eight can be also found in the examples computed in [SP07].
Figure 21: The computed variation of the (2 ,
5) orbit with mu = 0 . , . , .
05. On the leftwe see V ( t ) as a function of ω showing linear dependence in µ for the perturbations of V ( t )and quadratic dependence of δω . On the right the variation in the period of the first fullglacial cycle with V ( t ) showing linear dependence in µ for both.
7. More general dynamics of the PP04 model
The previous sections have allowed us to gain an analytical insight into thegeneral behaviour of the periodic solutions of the PP04 model for small periodicforcing, but give less information about the general behaviour of the system.Of course this is of most interest in a general discussion of how well the modelapplies to climate dynamics for which the periodic insolation forcing µ sin( ωt )takes larger values. We now make a systematic numerical study of this casewhich both confirms the predictions of the previous section for small µ , andalso which allows us to explore the rich dynamics of the forced PP04 system forthe case of larger values of the insolation forcing. A natural tool for analysing the PP04 climate model under periodic forcing isthe stroboscopic Poincare map P S mentioned in the last section. This map isdefined as follows Definition
Let the PP04 model be forced by the insolation function sin( ωt ),with state vector X ( t ) then P s X ( t ) ≡ X ( t + 2 π/ω ) . (59)Using this map we can construct a set of points x m defined by the iteration X m +1 = P S X m . (60)A ( m, n ) periodic orbit, as constructed above, then corresponds to an orbitwhich is an n − cycle ( X , X , . . . , X n − ) of P S for which X = P S X n − . m times. The nature of suchPoincar´e maps for Filippov flows has been studied some detail in [BBCK08]Chapter 7. In general, it follows from the theory presented in [BBCK08] thatthe map P S will be smooth if the intersection between the solution trajectoryand Σ is transversal. However it will lose smoothness if there is a grazing eventin the interval [ t, t + 2 π/ω ] leading to a non transversal intersection. As thevector field is continuous across Σ but has a derivative discontinuity, then themap will typically have a square root type behaviour close to the grazing point.We will explore the impact that this has on the dynamics of the PP04 model inmore detail in a forthcoming paper.The general dynamics of the PP04 system can now be studied by consideringthe iterations of the map P S . To do this we use a Monte-Carlo approach inwhich, for a given parameter, we take a random set of initial data (typicallyfor computations this set will have 5 members) and iterate the solution startingfrom points in this set forward. To do this calculation we take the smoothedsystem with η = 1000 in the approximation of the Heavyside function, and solveforward in time it using the Matlab code ode15s (with tolerance set to 1 e − set of random initial data, we obtaina Mont´e-Carlo plot of (hopefully) all of the possible Omega-limit sets. Thisgives significant insight into the overall dynamics of the system. It is convenientto represent the state of the whole system by plotting the values of the singlevariable F ( X i ) at the points X i . The advantage of this approach over the path-following methods used, for example, in the AUTO code [DCF + ω . Initially we take fixed small values of µ (consistent with the earlier analysis)and vary the value of ω . In Figure 22 we take µ = 0 .
05 and increase ω from0 .
08 to 0 .
13, plotting the omega-limit set of the resulting orbit in each case.It is convenient to represent these orbits by plotting the values of the function F . In this figure we can see a clear (1 ,
2) orbit for smaller values of ω and anequally clear (1 ,
3) orbit for the larger values. For ω ≈ .
107 there is a smallwindow of existence for the (2 ,
5) periodic orbit, and there is some evidence ofwindows of existence for more complex period motions. Away from these valueswe observe quasi-periodic behaviour. In Figure 23 we see (again for µ = 0 . ,
3) orbit changing to a quasi-periodic orbit when ω = 0 .
135 followed byan interval of quasi-periodic motion, which then turns into a (1 ,
4) orbit when ω = 0 . ,
7) orbit between the(1 ,
3) and (1 ,
4) orbits, and evidence of other periodic orbits.In Figure 24 we take the larger value of the forcing µ = 0 .
1. Again we seethe (1 , ,
5) and (1 ,
3) orbits with larger regions of existence, together withother types of more complex dynamics, but less evidence of a full quasi-periodicattractor. 40 igure 22: The Poincar´e section points of F on the omega limit set, as a function of ω with µ = 0 .
05 showing (as ω increases), a large window of existence for the (1 ,
2) periodicorbit, a much smaller window of existence for the (2 ,
5) periodic orbit close to ω = 0 . ,
3) periodic orbit. We can clearly seethe transition from quasi-periodic motion when ω < . ,
3) motion ata saddle-node bifurcation. All of the windows are separated by intervals of quasi periodicmotionFigure 23: The Poincar´e section of F on the Omega limit set, as a function of ω with µ = 0 . ,
3) and period (1 ,
4) motions separated by an interval of quasi periodicmotion, containing a small window with a period (2 ,
7) orbit, and evidence of other periodicorbits. igure 24: The Poincar´e section points of F on the Omega limit set, as a function of ω , with µ = 0 .
1, showing the (1 , ,
5) and (1 ,
3) periodic solutions and a variety of other types ofdynamics.
For a final calculation we take the physically relevant value of µ = 0 .
467 (seeSection 2 for a motivation of this value) and vary ω from 0 .
121 to 0 . ω , we see a (1 ,
2) periodic solutionand, for the larger values of ω , a (1 ,
3) periodic solution. There is no quasi-periodic behaviour in this case. Indeed, for a wide range of values of ω the(1 ,
2) and (1 ,
3) solutions co-exist. Two interesting transitions can be observedin this figure as ω increases. At ω = 0 .
122 the (1 ,
3) solution abruptly appears.The reason for this can be seen from studying the values of F . In particular, atthe bifurcation point, there is a value of t strictly within the glacial period, atwhich F ( t ) = 0. This is an example of a (non-smooth) grazing bifurcation (asmentioned above) at which the (1 ,
3) orbit suddenly starts to become physical.We will study this transition in more detail in a future paper. A (smooth)super-critical period-doubling bifurcation can also be seen at ω = 0 . ,
2) solution loses stability to a nearby (2 ,
4) orbit as ω increases.There is evidence of a period-doubling cascade close to this value. µ As a second calculation, we fix ω at the physically relevant value of ω = 0 . µ from zero. The resulting Mont´e-Carlo calculation ispresented in Figure 26. In this figure we see quasi-periodic behaviour for smallvalues of µ . The (1 ,
3) orbit arises at a saddle-node bifurcation at around µ =0 .
15 and persists until it is destroyed at a grazing bifurcation at µ ≈
1. For42 igure 25: A Mont´e-Carlo plot of the Omega-limit set of F for µ = 0 . ,
2) orbit and to the right a (1 ,
3) orbit. Various transitions between these orbitscan also be observed. These include a grazing bifurcation of the (1 ,
3) orbit at ω = 0 .
122 andan expansion of the (1 ,
2) orbit at ω = 0 . a short interval of values of µ there are coexisting (1 ,
2) and (1 ,
3) orbits. The(1 ,
2) orbit then persists until it too is destroyed at a grazing bifurcation when µ ≈ .
6. It co-exists with a (1 ,
1) orbit which loses stability at a period-doublingbifurcation when µ ≈ .
5. For larger values of µ we see only the (1 ,
1) periodicorbits, completely locked to the forcing. At the physically interesting value of µ = 0 .
467 (see Section 2) we see only a (1 ,
3) periodic orbit.
The co-existence, for example, of the (1 ,
2) and (1 ,
3) solutions when ω = 0 . µ = 0 . ,
4) and (1 ,
3) orbits when ω = 0 . V, A, C ) suchthat the omega-limit set of the iterations of the map P S is either the (1 , , ,
3) orbit. It is problematic to find the full three dimensional sets, so forconvenience we find a two-dimensional projection by fixing t initial = 0 , A = 0 . ,
3) orbit as ω increases from 0.124 to 0.128.43 igure 26: The Poincar´e section points on the Omega limit set of F , as a function of µ withfixed ω = 0 . µ increases from quasi-periodic,to (1 , ,
2) periodic orbits and then a (1 ,
1) periodic solution, with regions of co-existence.The (1 ,
3) orbit starts at a saddle-node bifurcation and terminates at a grazing bifurcation.The (1 ,
1) periodic solution shows evidence of a period-doubling bifurcation at µ ≈ . t initial = 0 , A = 0 .
55, when ω = 0 .
124 (left) and ω = 0 . µ = 0 . , ,
2) (left) or (2 , ω = 0 .
124 and plot ( t, F ) fora solution in which we take initial conditions in the green region but close to thered boundary with ( V (0) , A (0) , C (0)) = (0 . , . , . ,
3) periodic orbit, which thenultimately evolves to a (1 ,
2) orbit. We note that there is a dramatic change inthe behaviour of the system when t ≈ kyr . This occurs when there is a localminimum at which F < grazing transition [BBCK08].
Figure 28: The time evolution of F ( t ) for the system started close to the boundary ofthe domain of attraction. This figure shows the slow evolution from a (1 ,
3) periodic or-bit to a (1 ,
2) periodic orbit. Here ω = 0 .
124 and µ = 0 .
467 and the initial conditions are( V (0) , A (0) , C (0)) = (0 . , . , . As a separate calculation we take µ = 0 .
467 and ω = 0 . V (0) , A (0) , C (0)) = (0 . , . , . . In the resulting intermittent dynamics wesee a (2 ,
4) orbit evolve into a larger amplitude (1 ,
3) orbit in a manner whichqualitatively resembles that at the mid-Pleistocene transition. The sudden ex-pansion in the solution amplitude (and the consequent change in period) againseems to occur just after the function F grazes zero. We will return to this inthe forthcoming paper on grazing transitions in the PP04 model.
8. The implications of these results for climate modelling.
From the results that we have obtained for the PP04 model, we have shown thatif there is no insolation forcing on the system and − . < d < .
32, then thereis a periodic orbit of period of about 140 kyr. Numerically this orbit appears to45 igure 29: The time evolution of F ( t ) when ω = 0 .
128 and µ = 0 .
467 with initial conditions( V (0) , A (0) , C (0)) = (0 . , . , . . In this case we see a slightly unstable, low amplitude,(2 ,
4) orbit evolve into a larger amplitude (1 ,
3) orbit. be both stable and unique. The existence of this orbit suggests that the Earth’sclimate, if left alone without the contribution of insolation forcing, will haveperiodic glacial cycles. In these it will spend most of its time, say about 120kyr, in the glacial state and less time in the inter-glacial state. On the otherhand, if the d is greater than 0 .
32 or less than − .
72, we have stable equilibriaand the climate can get locked into either a glacial state or an inter-glacial state. (1 , orbit under changes to the inso-lation forcing. When purely periodic insolation forcing is introduced, and we consider the phys-ically relevant values of ( µ, ω ) = (0 . , . ,
3) pe-riodic orbit. This orbit has period 6 π/ω = 123 kyr, which is slightly longerthan the observed period of 100 kyr. (We note that 100 kyr is very close tothe period of the (2 ,
5) periodic orbit. However, we have not seen any evidenceof this orbit existing close to the realistic parameter values). From extensivenumerical experiments, for these parameter values, the (1 ,
3) orbit appears tobe unique, globally stable, and indeed strongly attracting, for all physical initialstates. The resulting orbit and a short transient is shown in Figure 30Of course this analysis has only been made for the case of periodic forcing . Inpractice the Milankovich cycles lead to quasi-periodic insolation forcing. Thestructural stability of the (1 ,
3) orbit constructed above means that this orbitpersists, appropriately perturbed to an invariant torus, when quasi-periodicforcing is introduced, with a small additional forcing. We demonstrate this byconsidering an insolation forcing of the form I ( t ) = µ sin( ω t ) + µ sin( ω t ) . igure 30: The evolution of the solution when ( µ, ω ) = (0 . , . ,
3) periodic solution. In this figure V is shown in red, A in green and C in blue. Provided that µ is not too large, the (1 ,
3) orbit in this case is replaced by aquasi-periodic orbit on a torus in the phase space close to the original periodiccurve. This is illustrated in Figure 31 which we compare with the above figureFigure 30. The study in more detail of the quasi-periodic forced PP04 modelwill be given later in a later paper, where we consider the break up of the tori forlarger forcing µ . Similar results for quasi-periodic forcing of the PP04 model(and other similar reduced climate models) are described in the paper by Ashwinet. al [ADCvdH18] (see also [Cru13]) in which apparently chaotic behaviour ofthe solutions was observed for certain types of quasi-periodic forcing. When µ = 0 .
467 and ω = 0 . , ω close to 0 .
128 we also see stable (1 , ,
4) solutions. Note that if ω = 0 .
128 then theperiod of the (1 ,
2) orbit is 98.17 kyr and of the (1 ,
3) orbit is 147 kyr. If ω is fixed and the initial data is taken close to the boundary of the domains ofattraction of these orbits, the we see transitions, for values of ω close to 0.128,both from (1 ,
3) orbits to (1 ,
2) orbits and from (period-doubled) (2 ,
4) orbits to(1 ,
3) orbits. The latter transition, in particular, has some resemblance to thequalitative changes in the behaviour of the climate at the MPT. During suchtransitions there is a long transient motion close to one form of periodic orbit,before the solution converges on the other. Such examples of transitions raisethe hope of understanding the MPT through a bifurcation type of analysis.47 igure 31: The time solution of the system showing the quasi-periodic forced solution withthe (1 ,
3) periodic solution evident for µ = 0 . ω = 0 . µ = 0 . ω = 0 . ,
3) periodic solution is perturbed to a more complex orbit ona torus. However, the basic form of the (1 ,
3) orbit remains.
However, much more work needs to be done on this to explore the varioustransitions possible given the large number of parameters that can be varied inthe PP04 model. We will be discussing, in particular, sudden transitions due tograzing bifurcations in a forthcoming paper.
9. Conclusions
In this paper we have made a first mathematical study using the theory of non-smooth dynamical systems of the (periodically forced) PP04 model for climatechange. This has revealed the existence of stable and unstable periodic orbits,with subtle domains of attraction and transitions between them. The stableorbits calculated for physically realistic values of the parameters persist undersmall additional quasi-periodic forcing and have a similar form to those of theobserved glacial cycles. The results make an interesting comparison to thoseof descriptions of the glacial cycle using smooth dynamics systems models, forexample [EKKV17].Much more work needs to be done on the PP04 model to understand fully boththe transitions in the whole of the parameter space and also the effect of ad-ditional larger terms in the quasi periodic forcing. Both of these will be thesubject of further work, which will look in more detail at the effect of grazingbifurcations on the stability of the orbits in the PP04 model and how these48grazing) transitions change when the insolation forcing is quasi-periodic. Fur-thermore additional work is needed to understand better the effect of includingadditional climatic terms into the PP04 model. However, we conclude that thePP04 model both has a rich structure as a discontinuous dynamical system, andis a plausible explanation of the glacial cycles. As such it deserves much furtherstudy.
Acknowledgement
This research was funded in part by an award from the Botswana Interna-tional University of Science and Technology (BIUST). We would like to thankProf. Rachel Kuske (Georgia Tech) and Prof. Paul Glendinning (University ofManchester) for many stimulating conversations related to this work, and theanonymous referees for their very insightful comments on an earlier version ofthis work.
BibliographyReferences [AD15] Peter Ashwin and Peter Ditlevsen,
The middle pleistocene transi-tion as a generic bifurcation on a slow manifold , Climate dynam-ics (2015), no. 9-10, 2683–2695.[ADCvdH18] Peter Ashwin, Charles David Camp, and Anna S von der Heydt, Chaotic and non-chaotic response to quasiperiodic forcing: Limitsto predictability of ice ages paced by milankovitch forcing , Dynam-ics and Statistics of the Climate System (2018), no. 1, 1–20.[AFO05] Jan Awrejcewicz, Michal Feˇckan, and Pawel Olejnik, On continu-ous approximation of discontinuous systems , Nonlinear Analysis:Theory, Methods & Applications (2005), no. 7, 1317–1331.[BBCK08] Mario Bernardo, Chris Budd, Alan Richard Champneys, and Pi-otr Kowalczyk, Piecewise-smooth dynamical systems: theory andapplications , vol. 163, Springer Science & Business Media, 2008.[CD10] Alessandro Colombo and Fabio Dercole,
Discontinuity induced bi-furcations of nonhyperbolic cycles in nonsmooth systems , SIAMJournal on Applied Dynamical Systems (2010), no. 1, 62–83.[CDBHJ12] Alessandro Colombo, M Di Bernardo, SJ Hogan, and MR Jeffrey, Bifurcations of piecewise smooth flows: Perspectives, methodolo-gies and open problems , Physica D: Nonlinear Phenomena (2012), no. 22, 1845–1860.[Cor08] Jorge Cortes,
Discontinuous dynamical systems , IEEE Controlsystems magazine (2008), no. 3, 36–73.49Cru12] Michel Crucifix, Oscillators and relaxation phenomena in pleis-tocene climate theory , Philosophical Transactions of the RoyalSociety A: Mathematical, Physical and Engineering Sciences (2012), no. 1962, 1140–1165.[Cru13] ,
Why could the ice ages be unpredictable , Clim. Past (2013), 2253–2267.[DBBC +
08] Mario Di Bernardo, Chris J Budd, Alan R Champneys, PiotrKowalczyk, Arne B Nordmark, Gerard Olivar Tost, and Petri TPiiroinen,
Bifurcations in nonsmooth dynamical systems , SIAMreview (2008), no. 4, 629–701.[DBH10] M Di Bernardo and SJ Hogan, Discontinuity-induced bifurcationsof piecewise smooth dynamical systems , Philosophical Transac-tions of the Royal Society A: Mathematical, Physical and Engi-neering Sciences (2010), no. 1930, 4915–4935.[DCF +
98] Eusebius Doedel, Alan Champneys, Thomas Fairgrieve, BjornSandstedte, and Xainjun Wang,
Auto97:c (continuation and bifur-cation software for ordinary differential equation, with homcont) ,Concordia University, Technical Report (1998).[Dij13] Henk A Dijkstra,
Nonlinear climate dynamics , Cambridge Uni-versity Press, 2013.[DSCW13] Bernard De Saedeleer, Michel Crucifix, and Sebastian Wieczorek,
Is the astronomical forcing a reliable and unique pacemaker forclimate? a conceptual model study , Climate Dynamics (2013),no. 1-2, 273–294.[EKKV17] Hans Engler, Hans Kaper, Tasso Kaper, and Theodore Vo, Dy-namical systems analysis of the maasch-saltzman model for glacialcycles , Physica D: Nonlinear Phenomena (2017), 1–20.[Gle16] Paul Glendinning,
Classification of boundary equilibrium bifur-cations in planar filippov systems , Chaos: An InterdisciplinaryJournal of Nonlinear Science (2016), no. 1, 013108.[GOH13] Antonio Garc´ıa-Olivares and Carmen Herrero, Simulation ofglacial-interglacial cycles by simple relaxation models: consistencywith observational results , Climate dynamics (2013), no. 5-6,1307–1331.[GST08] Marcel Guardia, TM Seara, and MA Teixeira, Topological equiv-alences for planar filippov systems , Talk during Problems in Non-smooth Dynamical Systems, University of Bristol (2008), 28–29.50GT00] Hezi Gildor and Eli Tziperman,
Sea ice as the glacial cycles cli-mate switch: Role of seasonal and orbital forcing , Paleoceanogra-phy and Paleoclimatology (2000), no. 6, 605–615.[Hel82] Isaac M Held, Climate models and the astronomical theory of theice ages , Icarus (1982), no. 2-3, 449–461.[HIS +
76] James D Hays, John Imbrie, Nicholas J Shackleton, et al.,
Vari-ations in the earths orbit: pacemaker of the ice ages , Science (1976), no. 4270, 1121–1132.[IBB +
93] John Imbrie, Andr´e Berger, EA Boyle, SC Clemens, A Duffy,WR Howard, G Kukla, J Kutzbach, DG Martinson, A McIntyre,et al.,
On the structure and origin of major glaciation cycles 2. the100,000-year cycle , Paleoceanography (1993), no. 6, 699–735.[JLP +
87] Jean Jouzel, Cl Lorius, JR Petit, C Genthon, NI Barkov,VM Kotlyakov, and VM Petrov,
Vostok ice core: a continuousisotope temperature record over the last climatic cycle (160,000years) , Nature (1987), no. 6138, 403.[KE13] Hans Kaper and Hans Engler,
Mathematics and climate , vol. 131,Siam, 2013.[KJea18] Till Kuhlbrodt, Colin Jones, and et. al.,
The lowresolution versionof hadgem3 gc3.1: Development and evaluation for global climate ,Journal of advances in modelling earth systems (2018), 2865–2888.[MA14] Takahito Mitsui and Kazuyuki Aihara, Dynamics between orderand chaos in conceptual models of glacial cycles , Climate dynam-ics (2014), no. 11-12, 3087–3099.[MCA15] Takahito Mitsui, Michel Crucifix, and Kazuyuki Aihara, Bifur-cations and strange nonchaotic attractors in a phase oscillatormodel of glacial–interglacial cycles , Physica D: Nonlinear Phe-nomena (2015), 25–33.[Pai98] Didier Paillard,
The timing of pleistocene glaciations from a sim-ple multiple-state climate model , Nature (1998), no. 6665,378.[Pai01] ,
Glacial cycles: toward a new paradigm , Reviews of Geo-physics (2001), no. 3, 325–346.[Pai17] , Climate science: Predictable ice ages on a chaotic planet ,Nature (2017), no. 7642, 419.51PJR +
99] Jean-Robert Petit, Jean Jouzel, Dominique Raynaud, Narcisse IBarkov, J-M Barnola, Isabelle Basile, Michael Bender, J Chappel-laz, M Davis, G Delaygue, et al.,
Climate and atmospheric historyof the past 420,000 years from the vostok ice core, antarctica , Na-ture (1999), no. 6735, 429.[PK08] Petri T Piiroinen and Yuri A Kuznetsov,
An event-driven methodto simulate filippov systems with accurate computing of slidingmotions , ACM Transactions on Mathematical Software (TOMS) (2008), no. 3, 13.[PP04] Didier Paillard and Fr´ed´eric Parrenin, The antarctic ice sheet andthe triggering of deglaciations , Earth and Planetary Science Let-ters (2004), no. 3-4, 263–271.[PRKK03] Arkady Pikovsky, Michael Rosenblum, Jurgen Kurths, and J¨urgenKurths,
Synchronization: a universal concept in nonlinear sci-ences , vol. 12, Cambridge university press, 2003.[Sim10] David John Warwick Simpson,
Bifurcations in piecewise-smoothcontinuous systems , vol. 70, World Scientific, 2010.[SM90] Barry Saltzman and Kirk A Maasch,
A first-order global model oflate cenozoic climatic change , Earth and Environmental ScienceTransactions of the Royal Society of Edinburgh (1990), no. 4,315–325.[SM91] , A first-order global model of late cenozoic climatic changeii. further analysis based on a simplification of co 2 dynamics ,Climate Dynamics (1991), no. 4, 201–210.[SP07] Frank Schilder and Bruce Peckham, Computing arnold tongue sce-narios , J. Comp. Phys. (2007), 932–951.[WWHM16] James Walsh, Esther Widiasih, Jonathan Hahn, and RichardMcGehee,
Periodic orbits for a discontinuous vector field arisingfrom a conceptual model of glacial cycles , Nonlinearity29