A Discontinuous ODE Model of the Glacial Cycles with Diffusive Heat Transport

We present a new discontinuous ordinary differential equation (ODE) model of the glacial cycles. Model trajectories flip from a glacial to an interglacial state, and vice versa, via a switching mechanism motivated by ice sheet mass balance principles. Filippov’s theory of differential inclusions is used to analyze the system, which can be viewed as a nonsmooth geometric singular perturbation problem. We prove the existence of a unique limit cycle, corresponding to the Earth’s glacial cycles. The diffusive heat transport component of the model is ideally suited for investigating the competing temperature gradient and transport efficiency feedbacks, each associated with ice-albedo feedback. It is the interplay of these feedbacks that determines the maximal extent of the ice sheet. In the nonautonomous setting, model glacial cycles persist when subjected to external forcing brought on by changes in Earth’s orbital parameters over geologic time. The system also exhibits various bifurcation scenarios as key parameters vary.


Introduction
Earth's climate has varied drastically over long time scales, cycling between warmer interglacial periods and very cold glacial states (ice ages). The beginning of a glacial cycle is characterized by a slow descent into a much colder world, in which massive ice sheets advanced into North America and Eurasia, fed by the moisture resulting from a corresponding drop in sea level of about 350 feet [1]. The ice ages are followed by a relatively rapid retreat into interglacial periods, leading to a sawtooth pattern evident in the paleoclimate data ( Figure 1). While large-scale climate models are able to provide insight into glacial cycle dynamics [2], the use of a comprehensive Earth System Model to accurately simulate even the last glacial cycle remains a challenge [3]. As such, approaches such as the analysis of data from deep ocean sediment cores [4] and Greenland and Antarctic ice cores [5], and the use of conceptual (or low order, reduced) climate models, also play important roles in the investigation of climate phenomena associated with glacial cycles [6][7][8][9][10][11][12][13][14][15]. As the time scales of related climate phenomena can differ to a great degree, conceptual models often involve a switching mechanism and nonsmooth model equations [13,[16][17][18][19][20]. The model introduced here will contain a climate switching mechanism, leading to an associated nonsmooth geometric singular perturbation problem. We prove the existence of a unique attracting (nonsmooth) limit cycle that corresponds to the glacial-interglacial oscillations.
Numerous questions regarding the glacial cycles remain unanswered, beginning with those related to the frequency and amplitude of the oscillations. Many believe the pacing of the ice ages is related to changes in Earth's orbital parameters, known as Milankovitch cycles [21]. For example, the obliquity (or tilt) of the Earth's spin axis has varied over the past 5 million years [22] with a period of about 41 kyr. It is known that, between 3 and 1.2 million years ago, the glacial cycles oscillated with a dominant period of 41 kyr as well [4].
Over the past 800 kyr, ice ages were of longer duration and the period of the glacial cycles grew to about 100 kyr, with the amplitude of the cycles increasing as well ( Figure 1). Interestingly, the eccentricity of the Earth's orbit has varied over geologic time with a period of roughly 100 kyr. Understanding this change to larger amplitude, smaller frequency glacial cycles between 1200 and 800 thousand years ago-known as the Mid-Pleistocene Transition (MPT)-and the nature of the climate response to the external forcing provided by the Milankovitch cycles, continues as a major open problem in paleoclimate. From a dynamical systems perspective, a phenomenon such as the MPT is testament to the complexity of the climate system and the often nonlinear response of the climate system to external forcing.
Additional questions concern internal climate feedbacks and the possibility of self-sustained oscillations, which represent the main focus of this paper. Foremost among these from a conceptual modeling perspective is positive ice-albedo feedback. An expanding ice sheet leads to a higher albedo and reduced absorption of solar radiation, thereby reducing temperature and further enhancing growth of the ice sheet. Similarly, a retreating ice sheet leads to an increase in absorption of solar radiation and higher temperatures, contributing to further ice sheet loss.
The role of ice-albedo feedback as it relates to zonally averaged (by latitude) surface temperature was studied in the seminal work of Budyko [7] and Sellers [14] in 1969. Budyko's surface temperature model was coupled with an ODE for the evolution of the edge of the ice sheet (known as the ice line) in [15], where it was shown the ice line either converges to a small stable ice cap equilibrium position or descends to the equator (with parameters aligned with the modern climate). The latter outcome is known as a runaway snowball Earth episode.
The potential for self-sustained climate oscillations might arise via consideration of competing feedbacks associated with ice-albedo feedback [8,9]. For example, at the onset of a glacial age, a descending ice sheet leads to an increase in the surface temperature gradient, due to the increased albedo in the polar region. This leads to a positive feedback in that the poleward transport of heat and moisture by the atmosphere, with the moisture precipitating out as snow on the glacier, is enhanced, all things being equal [23]. In this way, the ice sheet continues to grow.
A competing, negative feedback is that the efficiency of the meridional heat and moisture transport is decreased as the ice sheet advances. In the conceptual model setting, North [24] proposed that the efficiency of the meridional transport might have two values, the smaller of which corresponds to the portion of the planet covered by ice. That is, a planet with extensive ice cover will be less efficient in the transport of heat and moisture poleward (also see [25]). The thinking is that an increase in sea ice during a glacial advance serves to reduce heat advection in the oceans. Due to atmospheric-oceanic coupling, together with extensive sea ice, evaporation and the transport of moisture will each be impaired.
Similarly, as albedo decreases during a glacial retreat, the meridional temperature gradient weakens, impacting the poleward transport of moisture. However, the efficiency of the meridional transport might be assumed to increase as the sea ice retreats, with ocean advection and the coupling between the atmosphere and the ocean enhanced. Thus, the efficiency of the diffusive transport across latitudes can be viewed as a negative feedback with respect to ice sheet growth.
We investigate the interaction between this negative feedback and the positive temperature gradient feedback described above in the conceptual model, ODE setting. We do this in part by assuming the diffusive meridional transport efficiency parameter is smaller during a glacial advance, relative to that during a glacial retreat, due to the ice cover.

Results
We present a new conceptual model of the glacial cycles in which the latitudinal transport of heat is modeled as a diffusive process. The zonal surface temperature PDE is coupled to the ice line, as in [15]. The spectral method is used to approximate the resulting infinite dimensional system with a finite system of ODEs. Standard invariant manifold theory [26] is used to reduce this finite system to a single governing equation on a slow manifold, corresponding to the evolution of the edge of the ice sheet.
Consideration of ablation and accumulation zones on the ice sheet leads to the introduction of a second model variable, as in [20]. The incorporation of a switching mechanism based on an ice sheet mass balance equation into this two-dimensional system, as described below, leads to a differential inclusion. Using Filippov's approach to differential inclusions [27], a unique attracting (nonsmooth) limit cycle is shown to exist, with the Contraction Mapping Theorem the major tool used in the proof. This periodic orbit, a consequence of the competing feedbacks discussed in the Introduction, represents internal, self-sustained climate oscillations.
Similar in spirit to Welander's mixed ocean layer heat-salt oscillator [28], two virtual equilibria play a fundamental role in the limit cycle produced here. Varying a model parameter leads to a type of nonsmooth Hopf bifurcation, arising as a border collision bifurcation, as is discussed along with additional bifurcation scenarios in the following.
It is of interest to note that forcing our model with variations in eccentricity and obliquity does not significantly affect the glacial-interglacial oscillations in a qualitative sense. The dynamics of this nonautonomous, nonsmooth model are presented below. This paper is organized as follows. The model equations are derived in Section 3. The switching mechanism and the set up of the Filippov flow are presented in Section 4. The nonsmooth model flow is analyzed in Section 5, in which the existence of a unique attracting limit cycle is proved. We revisit, and place within the context of the model, the competing feedbacks discussed above by considering the flux at the albedo line throughout the glacial cycle in Section 6. We discuss bifurcations exhibited by the model in Section 7, while the effect of forcing the model with Milankovitch cycles is presented in Section 8. The concluding section contains comments pertaining to open questions, both of a mathematical nature and as they relate to Earth's glacial cycles, that might be investigated in a conceptual fashion with our model.

The Zonal Average Surface Temperature Equation
The model development begins with the equation for the annual average surface temperature T = T(y, t) ( • C) at y = sin θ, where θ denotes latitude. Hence, T is a zonal average, with a single spatial variable y. Due to an assumption of symmetry of the temperature distribution about the equator as in [7], one takes y ∈ [0, 1]. Note that y = 0 and y = 1 correspond to the equator and North Pole, respectively.
The equation of a type first introduced by Budyko [7] and Sellers [14], is derived from energy balance principles. The behavior of equilibrium solutions of Equation (1) in response to changes in a given parameter was well studied by North and colleagues in the 1970s and early 1980s [24,[29][30][31].
The units on each side of (1) are watts per square meter (W/m 2 ). The constant R, with units joules per square meter per degree Celsius, denotes the global average heat capacity of the Earth's surface; the value R = 0.5 × 10 9 J/(m 2 • C) used here lies within the range of estimates provided in [32].
Q (W/m 2 ) is the global average annual incoming solar radiation (or insolation). While Q = Q(e) depends upon the eccentricity e of the Earth's orbit over long time scales [11], varying between roughly 342.95 and 343.52 over the past 5 million years [22], we set Q = 343 W/m 2 in this section.
The function is the insolation distribution function [11], accounting for the fact that, on Earth, tropical latitudes receive greater annual mean insolation than do polar regions. The function s(y) satisfies the normalizing condition 1 0 s(y)dy = 1.
In addition, s(y) depends on the obliquity β of the Earth's spin axis. The obliquity has varied between roughly 22 • and 24.5 • over geologic time, with the current value of β approximately equal to 23.4 • .
The function where α 1 < α 2 , represents the reflectivity or albedo of the Earth's surface. The ice line is denoted by η; following Budyko in [7], one assumes η separates latitudes of high albedo α 2 poleward of η from latitudes of lower albedo α 1 equatorward of η. Taken together, the expression Qs(y)(1 − α(y, η)) in (1) models the energy into the climate system. The (A + BT)-term represents the outgoing longwave radiation (OLR), or energy out of the climate system. The parameters A > 0 (W/m 2 ) and B > 0 (W/(m 2• C) have been estimated empirically [33] (see Table 1). That the Earth's OLR is essentially linear over a wide range of surface temperature values is presented in [34].
The final term in (1) models the meridional transport of heat by the atmosphere and ocean as a diffusive process, with D > 0 (W/(m 2 • C) the diffusion coefficient. Recalling y = sin θ and that T(y, t) has no longitudinal dependence, the spherical diffusion operator takes the form (see, for example, [24]). A mean annual temperature model with symmetry about the equator must satisfy boundary conditions, namely, that the gradient of the temperature vanishes at the North Pole and at the equator. That this is the case in our approach is shown in Section 3.2. While Budyko modeled meridional heat transport in [7] via the expression where C > 0 and T ave is the annual global mean surface temperature, Equation (4) has the advantage that the Legendre polynomials P n (y) are eigenfunctions of the spherical diffusion operator d dy (1 − y 2 ) d dy P n (y) = −n(n + 1)P n (y), n = 0, 1, 2, ....
We make use of this observation in the following section.

The Spectral Approach
We approximate Equation (1) via a finite system of ODEs by taking advantage of (6). For a fixed ice line η, this is the approach used by North in works such as [24,29]. (In the following section, the ice line will be allowed to evolve over time, coupled to the temperature equation as in [15].) Full details of the following derivation can be found in [35]. We assume the obliquity β = 23.4 • in the following.
Nadeau and McGehee provide a mechanism for the computation of the coefficients s 2i (as functions of the obliquity β) in [36]. We also note the dependence of α 2i on the ice line η, arising from the substitution of albedo function (3) into the integrand in (11). In particular, α 2i is a polynomial in η of degree 4N + 1. Substituting (7)-(9) into PDE (1) and equating coefficients of P 2i , one arrives at a system of N + 1 ODEs that can be placed in the form where In this fashion, solutions to the N + 1 equations in system (12) serve to approximate solutions to Equation (1). For future reference, note that, for a fixed ice line η, T 2i (t) → f 2i (η) as t → ∞ for each i, leading to a globally attracting curve of rest points and an associated equilibrium temperature distribution each parametrized by η. E. Widiasih introduced an equation for the evolution of the ice line η that was coupled to the zonal average surface temperature in a simple fashion in [15]. This additional equation will aid consideration of the feedbacks discussed in the Introduction. We incorporate this model enhancement in the following section.

A Dynamic Ice Line
It is natural to allow the ice line to respond to changes in the zonally averaged surface temperature, a concept introduced and formalized in [15]. The equation considers the temperature at the ice line η, relative to a critical temperature T c . If T(η, t) > T c , η increases as the ice line moves poleward, while η decreases and the ice line moves equatorward if T(η, t) < T c . Coupling Equation (15) with system (12), and evaluating (7) at η, leads to the system of N + 2 ODEs Note that when ρ = 0 we return to the fixed ice line case and a globally attracting curve of rest points on which T 2i = f 2i (η), as described above. Via standard invariant manifold theory [26], this curve of rest points perturbs to a nearby globally attracting invariant manifold, on which system (16) can be approximated by Equation (16b) That is, Equation (17) is a good approximation to the coupled temperature-ice line system (16) if it is assumed that the ice line changes slowly relative to changes in the temperature. Evidence that this latter assumption is reasonable is provided in [12]. Graphs of h(η), a polynomial of degree 4N + 3, are plotted in Figure 2 for various D-values (the use of subscripts G (glacial) and I (interglacial) in Figure 2 is fully explained in Section 4). Note that, for D = 0.3 and N = 1, Equation (17) has an unstable equilibrium with the ice line at roughly η = 0.2 corresponding to a very large ice cap, and a stable, small ice cap equilibrium at about η = 0.8 (the blue curves in Figure 2a,b). To each of these equilibrium η-vales, there is an associated equilibrium temperature profile T * η (y) given by Equation (14).
A diffusion coefficient value of D = 0.394 is required to place the stable equilibrium near η = 0.94 when N = 1, roughly positioned to match the Earth's current polar ice cap extent and climate (the red curves in Figure 2a,b). That is, if at equilibrium η = 0.94 when N = 1, the choice of D = 0.394 yields a global mean surface temperature of 15 • C, approximately equal to today's value. We note the existence of an unstable equilibrium η-value between η = 0.94 and the North Pole when η = 0.394 ( Figure 2b). While working solely with equilibrium solutions of temperature Equation (1) and a fixed ice line, and with "stability" considered with respect to changes in the parameter Q rather than with changes in η as considered here, North analyzed the existence of this small unstable ice cap in [30]. In particular, North discussed the connection between the use of the diffusive heat transport mechanism (4) and this small unstable ice cap, while also commenting on the lack of such an equilibrium solution when using Budyko's transport term (5).
We note the stable-unstable small ice cap pair is lost through a saddle node bifurcation as D increases through roughly D = 0.4 (see Figure 3). Hence, if the meridional heat transport is sufficiently efficient, the climate tends to either an ice-free Earth or a snowball Earth, depending on the initial position of the ice line relative to the large unstable ice cap η-value. Finally, the use of higher order approximations does not change the qualitative nature of solutions of Equation (17), as illustrated in Figure 2c, for which all expansions were taken out to degree 6.

Incorporating Accumulation and Ablation Zones
Note the role of η as defined in albedo function (3) is more accurately viewed as an "albedo line", or more germanely to our model, as a snow line. At locations above η, the surface has a higher albedo than at latitudes south of η. We now recast η as the snow line, as in [20], considering the region poleward of latitude η to be snow-covered. This leads to the introduction of the new variable ξ to designate the ice edge. The region between ξ and η is interpreted as the ablation zone, where wasting of the glacier occurs. The surface between η and the North Pole is an accumulation zone, comprised of latitudes at which snowfall exceeds ablation processes (see Figure 4).  Motivation for this modeling choice comes from numerous sources, the first being that accumulation and ablation of ice play an important role in the advance, retreat, and size of a glacier [37]. In addition, study [2], in which an ice sheet model is driven by climate parametrizations obtained by a GCM in an effort to simulate the last four glacial cycles, found that the larger the ice sheet was at equilibrium, the larger was the ablation zone. Moreover, Abe-Ouchi et al. also found that an increased ablation rate was necessary during the retreating phase in order to faithfully reproduce the qualitative features of the paleoclimate data and, in particular, the rapid interglacial retreats that are known to have occurred.
Additional motivation was provided by the fact that areas near the ice sheet margins can become stagnant, allowing for the growth of superglacial forests [38] which reduce the local surface albedo. Other factors contributing to a lower surface albedo in the ablation zone include the presence of dust and debris [39], and the existence of "aging snow" [40]. Indeed, Gallée et al. also found that an increased ablation rate during glacial retreat-here due to the lower albedo of old snow and ice-was crucial in simulating the last glacial cycle. More currently, increased ablation due to the reduced albedo of darker bare ice is playing a role in the present day mass loss of the Greenland Ice Sheet [41].
We thus think of the surface poleward of the snow line η as snow-covered and having the high albedo α 2 , while regions equatorward of η have a lower albedo α 1 . The recasting of η as the snow line continues to accurately reflect our choice of albedo function (3).
The evolution of the ice edge ξ is, conceptually, driven by a mass balance principle. Although we do not explicitly consider ice volume and mass here, these quantities could be deduced from ξ were one to assume a specific ice sheet shape with longitudinal symmetry, as in [9,42].
Assume snow is accumulating between η and the pole at a (dimensionless) rate a, while ablation occurs between ξ and η at a (dimensionless) rate b. We then set Note that, when the accumulation a(1 − η) is greater than the ablation b(η − ξ), dξ/dt < 0 and the ice edge moves equatorward. If ablation supersedes accumulation, dξ/dt > 0 and the ice edge retreats poleward. We thus consider the system of equations Note the model allows for the possibility that η < ξ, as would occur during a glacial advance when new snow falls south of ξ. In this case, Equation (19b) dictates that ξ follows the snow line equatorward as the snow turns into ice over time.
A final strongly motivating work was that of Welander in [28], in which a simple mixed ocean layer "flip-flop" heat-salt oscillator was presented. As in [28], our model incorporates a mechanism, based on the differing ablation rates for glacial advance and retreat as discussed above, that enables switching from one regime to the other. The corresponding Filippov flow is described in the following section.

The Flip-Flop and Associated Filippov Flow
The theory of ODEs with discontinuous vector fields continues to advance, often based on the work of Filippov [27] and motivated in part by mathematical challenges arising in the modeling of systems that include phenomena such as switches and collisions [43,44]. In terms of our model, we preface the definition of the discontinuous vector field by first focusing on the diffusion coefficient D, which appears in each of the functions f 2i (η), i > 0, in (13). It follows that the function h(η) in (17) can be thought to be parametrized by D. Recall that an increase in D increases the efficiency of the meridional heat transport. We saw previously (Figure 2a) that an increase from D = 0.3 to D = 0.394 leads to a warmer world with a stable albedo line closer to the North Pole.
As mentioned in the Introduction, North put forth the possibility of using two D-values, with the smaller value corresponding to the ice-covered portion of the surface. While North only considered the fixed albedo line setting, system (19) has the advantage that the albedo-and ice-lines evolve with time. One might then imagine the value of D decreasing with an advancing ice sheet and increasing during a glacial retreat, playing its role in the negative feedback discussed in the Introduction. Here, we move to the discontinuous limit as it were, assuming a smaller constant value D = D G throughout the glacial advance, and a larger constant value D = D I throughout the retreat to an interglacial.
We will thus consider two climate regimes, one an "interglacial period" with D = D I = 0.394, and a "glacial period" for which D = D G = 0.3. We let h I (η) and h G (η) denote the function h(η) in Equation (17) when D = D I and D = D G , respectively. This represents step one in building the flip-flop characteristic of our model. We note the analysis to follow holds for any values D G < D I for which there are distinct small stable ice cap equilibria.
We pause before developing the model switching mechanism to revisit the positive ice albedo-temperature gradient feedback discussed in the Introduction. Let η = η G and η = η I denote the stable equilibrium solutions of Equation (17) when D = D G and D = D I , respectively, as in Figure 2. Let T * G (y) and T * I (y) denote the corresponding equilibrium temperature profiles (14), which for N = 1 are plotted in Figure 5. Note the temperature gradient during the glacial period (η = η G ) is larger than that during the interglacial period (η = η I ). We thus see that, at equilibrium, the model exhibits the desired effect associated with the positive ice albedo-temperature gradient feedback, namely, larger ice sheets due to enhanced poleward moisture transport brought on by an increased temperature gradient. As for the switching mechanism from one climate state to the other, we recall the numerous studies mentioned above that indicate an increased ablation rate during glacial retreats is needed to simulate the glacial cycles in a qualitatively accurate manner. Our second step will be to select (dimensionless) ablation rates b G < b I to incorporate this assumption. The two climate regimes are then modeled by the glacial advance system and the interglacial system We set B = {(η, ξ) : η ∈ [0, 1], ξ ∈ [0, 1]}, and we define the vector fields and Let φ G ((η, ξ), t) and φ I ((η, ξ), t) denote the flows associated with systems (20) and (21), respectively. To switch climate regimes, we select a critical ablation parameter b ∈ (b G , b I ). When b(η − ξ) − a(1 − η) > 0 (i.e., the mass balance is negative), we assume the ice sheet is retreating and the system evolves according to (21). The system is in glacial advance mode when b(η − ξ) − a(1 − η) < 0 (i.e., the mass balance is positive). In this fashion, the set on which b(η − ξ) − a(1 − η) = 0, namely becomes a switching manifold (or discontinuity boundary) [45]. A trajectory for either system (20) or (21) that intersects Σ in a crossing region switches to the alternate regime, as discussed in detail in Section 5.
In spirit, this switching from one flow to the other is entirely analogous to that found in Welander's mixed ocean layer model [28]. More formally, Σ separates B into the domains Note that each of V G and V I extend smoothly to Σ. For x = (η, ξ) ∈ B, consider the differential inclusion [45] dx dt While in S 1 solutions have flow φ G , while trajectories evolve under the flow φ I when in S 2 . For x ∈ Σ, one assumes dx/dt lies in the convex hull of the two vectors V G (x) and V I (x). A solution to (23) in the sense of Filippov is an absolutely continuous function x(t) satisfying dx/dt ∈ V(x) for almost all t (note x(t) is not differentiable at times t for which x(t) enters or leaves Σ). As V G and V I are continuous on Σ ∪ S 1 and Σ ∪ S 2 , respectively, the set-valued map V(x) is upper semicontinuous, as well as closed, convex, and bounded for all x ∈ B and all t ∈ R. These properties ensure that, for each x 0 ∈ Int(B), there is a solution x(t) to differential inclusion (23) defined on an interval [0, t f ] with x(0) = x 0 [45].

Stable Virtual Equilibria
Let ψ G (η, t) denote the flow associated with Equation (20a), and let ψ I (η, t) denote the flow associated with Equation (21a). Recall that h G and h I are polynomials of degree 4N + 3, differing only in the choice of diffusion coefficients D G < D I .
The analysis of discontinuous system (23) is aided by the observation that the evolution of η over time is independent of that of ξ. In particular, the Jacobian matrix for the vector field V G is (and similarly for the vector field V I ). As − b G /R < 0 and − b I /R < 0, the stability of any equilibrium solutions (η, ξ) of systems (20) and (21) will depend entirely on the signs of h G (η) and h I (η). With N = 1 and with parameter values as in Table 1, the degree seven polynomial h G (η) admits two equilibrium solutions η G,u ≈ 0.197 and η G,s ≈ 0.789 in the interval [0, 1], with h G (η G,u ) > 0 and h G (η G,s ) < 0 (see Figure 2). Thus, system (20) has a saddle Q G,u = (η G,u , ξ G,u ) and a stable node Q G,s = (η G,s , ξ G,s ). While ξ G,u < 0, the equilibrium Q G,s lies in B as ξ G,s ≈ 0.641.
As only the stable node for system (20) plays a role in the analysis, we let Q G and η G denote Q G,s and η G,s , respectively, in all that follows.
Note the line is invariant under the glacial flow φ G , with solutions on L 1 converging monotonically to Q G over time. In addition, every φ G -trajectory having initial condition in the region above Σ and to the right of L 1 , namely converges to Q G under the flow φ G . This follows from the fact that a ψ G -trajectory with initial condition η 0 ∈ [η G , 1] decreases monotonically to η G (recall Figure 2), and hence for such an η 0 the corresponding φ G -trajectory of (η 0 , ξ 0 ) satisfies ξ(t) → ξ G due to the nature of Equation (20b). Due to the ablation rate assumption b G < b, note for η < 1. This implies that the ξ-nullcline for vector field V G lies below the switching manifold Σ. At points above the ξ-nullcline system, (20) satisfies dξ/dt < 0 (and dη/dt < 0 for η > η G ); hence, φ G -trajectories starting in S 1 move monotonically to the left and downward as they converge to Q G (see Figure 6).  Figure 6. The Filippov flow set-up. Above Σ and to the right of L 1 , φ G -trajectories move left and downward (blue arrows). Below Σ and to the left of L 2 , φ I -trajectories move right and upward (red arrows). The blue and red dashed lines are the ξ-nullcines for systems (20) and (21), respectively. Vector fields V G (blue) and V I (red) are schematically indicated on Σ. We refer the reader to the text for further details.
Every φ G -trajectory starting in S 1 will thus cross Σ before approaching Q G and hence, following Filippov's convention, switch to the interglacial flow φ I as the trajectory enters S 2 . The fact the ξ-nullcline for system (20) lies below Σ implies Q G is a virtual equilibrium point [43]. (ii) x is a virtual equilibrium point of (23) if either V G (x) = 0 and x ∈ S 2 , or if V I (x) = 0 and x ∈ S 1 . (iii) x is a boundary equilibrium point of (23) if V G (x) = V I (x) = 0 and x ∈ Σ.
We now consider the interglacial flow. With N = 1 and with parameter values as in Table 1, the polynomial h I has three zeros η 1 ≈ 0.22, η 2 ≈ 0.936 and η 3 ≈ 0.955, with h I (η 1 ) > 0, h I (η 2 ) < 0 and h I (η 3 ) > 0 (recall Figure 2). Thus, system (21) has two saddles (at which η equals η 1 and η 3 , respectively), with a stable node Q I = (η I , ξ I ) in-between. The line is invariant under the flow φ I , with φ I -trajectories on L 2 converging monotonically to Q I over time.
The ξ-nullcline for vector field V I satisfies for η < 1 due to the assumption b I > b, and hence lies above Σ. We thus have that Q I is a virtual equilibrium for the flow φ I since Q I ∈ S 1 . Let Λ denote the strip and set Λ i = Λ ∩ S i , i = 1, 2. Since dη/dt > 0 and dξ/dt > 0 for system (21) in Λ 2 , φ I -trajectories starting in Λ 2 move monotonically to the right and upward as they approach Q I , and hence must intersect and cross Σ. In this fashion, the systems "flip-flops" between glacial and interglacial periods in this region of phase space.

Tangencies on the Switching Manifold
To gain further insight into the Filippov flow, it is important to note where on Σ the vector fields V G and V I are each tangent, given that such points may bound a sliding region [43]. To that end, a vector normal to Σ is Note that the line s G has positive slope (b G − b < 0), with s G (0) < 0 and s G (1) = 0. Recalling the plot of h G in Figure 2, the equation h G (η) = s G (η) thus has a unique solution η = µ G with η G < µ G . Moreover, as decreases to 0, µ G will decrease monotonically toward η G as the line s G approaches the η-axis. One may check V G · N > 0 at x 1 = (η, ξ) ∈ Σ for which η ∈ [η G , µ G ), so that V G points into S 2 at such x 1 . Additionally, V G points into S 1 at x 2 = (η, ξ) ∈ Σ for which η ∈ (µ G , 1], since V G · N < 0 at such x 2 (see Figure 6).
Similarly, V I · N = 0 if and only if h I (η) = s I (η), where a line with negative slope (b I − b > 0) satisfying s I (0) > 0 and s I (1) = 0. Hence, h I (η) = s I (η) has a unique solution η = µ I ∈ (η G , η I ), with µ I increasing to η I as decreases to 0. Moreover, V I points into S 1 at x 1 = (η, ξ) ∈ Σ for which η ∈ (µ I , η I ] (as (V I · N)(x 1 ) < 0), while V I points into S 2 at In the analysis to follow, we choose small enough to ensure µ G < µ I , which is possible given η G < η I . The points u G = (µ G , γ(µ G )) and u I = (µ I , γ(µ I )) on Σ are known as invisible fold singularities [44]. This means the φ G -trajectory passing through u G at t = t 0 is both tangent to Σ at u G and contained in S 2 for t near t 0 . Similarly, the φ I -trajectory passing through u I is both tangent to Σ at u I and contained in S 1 on some time interval.
Label the points z G = (η G , γ(η G )) and z I = (η I , γ(η I )) on Σ. The above discussion implies that φ I -trajectories starting in Λ 2 cross Σ moving to the right and up between u I and z I , while φ G -trajectories starting in Λ 1 cross Σ between z G and u G moving to the left and down. The portions of Σ between z G and u G and between u I and z I are known as crossing regions for the Filippov flow (23). As vector fields V G and V I are each smooth, Filippov solutions to differential inclusion (23) that pass through the crossing regions are unique.
The portion of Σ between u G and u I , known as a repelling sliding region [43], plays no role in the following analysis since Filippov trajectories beginning in Λ 1 , Λ 2 or in one of the two crossing regions never intersect this sliding region.

The Filippov Return Map and a Unique Limit Cycle
As shown in the preceding section, each φ G -trajectory starting in I 1 crosses Σ in I 2 . For each x ∈ I 1 , there is a unique time t (x) such that φ G (x, t (x)) ∈ I 2 . For future reference, note t (x) → ∞ as → 0. Define the mapping r G, : Similarly, for each x ∈ I 2 , there is a unique t (x) with φ I (x, t (x)) ∈ I 1 , since the φ I -trajectory of any such x crosses Σ in I 1 . This defines the map r I, : The composition r = r I, • r G, : I 1 → I 1 is then a continuous return map for the Filippov flow. Via our choice of the compact interval J, r has a fixed point in I 1 by Brouwer's Fixed Point Theorem. This implies the existence of a (nonsmooth) periodic orbit Γ for the Filippov system for any > 0 for which µ G < µ I .
That the periodic orbit Γ is unique and attracting for sufficiently small is easily motivated. Note that, when = 0, the φ G -trajectory of z I = (η I , γ(η I )) ∈ I 1 limits on the regular equilibrium point (η G , γ(η I )) ∈ L 1 . Similarly, the φ I -trajectory of z G = (η G , γ(η G )) ∈ I 2 limits on the regular equilibrium point (η I , γ(η G )) ∈ L 2 . For sufficiently small, the periodic orbit Γ is approximated by the rectangular "orbit" comprised of the φ G -trajectory of z I , the portion of L 1 from (η G , γ(η I )) down to z G , the φ I -trajectory of z G , and the portion of L 2 back up to z I (akin to Figure 5 in [20]).
With this insight, the return map r : I 1 → I 1 can be shown to be a contraction map as diam(r (I 1 )) can be made as small as we please by decreasing sufficiently, implying the existence of a unique fixed point to which all r -orbits converge. Table 1, there existsˆ > 0 so that, for all ∈ (0,ˆ ], the Filippov system (23) admits a unique attracting limit cycle.

Theorem 1. With parameters as in
Proof. Let c 1 ∈ (0, 1). Let x = (η 0 , ξ 0 ) be any point in I 1 . Recall ψ G (η 0 ) → η G as t → ∞, where ψ G is the flow associated with ODE (20a), here with initial condition ψ G (0) = η 0 . As the interval J used to define I 1 is compact, there exists a time T 1 so that for all t ≥ T 1 , for all η 1 , . By the continuity of φ G with respect to initial conditions and time, there exists δ = δ(x) > 0 so that, for all y ∈ B δ(x) ∩ I 1 , t (x) (y) > T 1 , where r G, (x) (y) ∈ I 2 . Note that, for any ≤ (x), t (y) > T 1 for y ∈ B δ(x) ∩ I 1 .
With parameter values as in Table 1, recall ODE (21a) has an unstable equilibrium η I,u lying between the stable equilibrium η I and the North Pole. The line is invariant under the interglacial flow φ I , with each φ I -trajectory on L 3 converging to Q I,u = (η I,u , ξ I,u ) in forward time. Thus, any point x = (η, ξ) ∈ S 2 with η G ≤ η ≤ η I,u will be in the stable set of the unique attracting limit cycle Γ, since the trajectory φ I (x, t) must intersect I 1 . At points x = (η, ξ) in S 2 and to the right of L 3 , however, vector field V I points to the right. Hence, a φ I -trajectory with initial condition x may hit the boundary η = 1 prior to intersecting Σ, and thus would not converge to Γ. We comment further on trajectories intersecting the boundary of B in the concluding section. The small unstable ice cap exists due to the choice D = D I = 0.394. As D decreases from D I = 0.394, the small unstable ice cap disappears, as can be gleaned from the bifurcation plot in Figure 3, with the ice-free Earth moving from "stable" to "unstable" as D decreases. If D I were chosen so that the small unstable ice cap equilibrium at η I,u did not exist, the stable set of the unique attracting limit cycle Γ would include all x = (η, ξ) ∈ B with η greater than the η-value corresponding to the large unstable interglacial ice cap.

Feedbacks and Meridional Flux at the Albedo Line
Placing the limit cycle depicted in Figure 7 in the context of the climate feedbacks discussed in the Introduction, we consider first the glacial advance. At the onset of the glacial age, the albedo in the polar region increases, leading to a steeper temperature gradient. This in turn enhances the meridional transport of moisture to northern latitudes where it precipitates out as snow, leading to further ice sheet growth. The albedo line η descends and, over time, the ice edge ξ follows, with the ice sheet mass balance positive in this regime.
As the ice edge descends, the assumption of reduced meridional transport efficiency comes into play, inhibiting the positive ice albedo-temperature gradient feedback. Accumulation on the ice sheet is reduced and eventually comes to a halt, while inertia keeps the ice edge ξ from responding quickly [9]. Thus, the ablation zone between ξ and η expands prior to the retreat, as is the case with the large model presented in [2].
The expansion of the ablation zone triggers the switch to the interglacial state as the mass balance becomes negative. In this mode, transport efficiency is increased, whereas the meridional temperature gradient is reduced. The albedo line moves poleward with the ice edge slowly following. The effect of the increase in the diffusion coefficient can be seen in the fact the albedo line stops its retreat shy of the North Pole as, eventually, snow again begins to precipitate out and the cycle begins anew.
The behavior described above can be summarized by considering the total meridional heat flux across the albedo line η over time. The flux at η is computed by integrating the diffusive heat transport term in (1) with respect to y, while making use of expression (7) for the temperature. The total heat flux moving poleward across latitude y is proportional to Using the assumption T 2i = f 2i (η) for each i and evaluating (26) at η yields Expression (27) is then proportional to the total meridional heat flux across η. The plot of (η) is shown in Figure 8. Note the increase in the poleward flux at η coincides with the advance of the albedo line. Due to the associated extensive ice cover, the transport efficiency is reduced, and, while the ice sheet edge continues to descend, the flux and the albedo line remain steady. When the system flips to the interglacial state, the poleward flux at η decreases as the albedo line retreats. Although the temperature gradient now weakens, the increased efficiency ensures moisture reaches northern latitudes, and eventually the albedo line stabilizes with the ice edge retreating. This simple flip-flop model serves to illuminate the interplay between the two competing feedbacks.

Bifurcation Scenarios
Recall that the ablation rates were chosen to satisfy b G < b < b I . As b G increases to b, the V G ξ-nullcline approaches the switching boundary Σ (see Figure 6 and Equation (24)). When b G = b, the equilibrium point Q G lies on Σ; when b G > b, Q G lies in S 1 and thus becomes a visible equilibrium. Similarly, Q I can be made to pass through Σ by letting b I decrease through b.
As b G increases through b, the limit cycle Γ is lost through a border collision bifurcation [44], a type of nonsmooth Hopf bifurcation. For b G > b, a trajectory starting in Λ 1 will converge to Q G , never intersecting Σ and so never triggering the flip-flop. In Figure 9a, we plot the limiting maximum and minimum values of the norm of trajectories beginning in Λ 1 as b G passes through b (assuming is sufficiently small).
with s replaced by the negative of the mass balance term, b(η − ξ) − a(1 − η). Similarly, replacing b G and b I in (28) with D G and D I , respectively, leads to a smooth transition between the glacial advance and retreat diffusion constants. In this way, system (19) becomes a smooth dynamical system that approximates our flip-flop model for large M-values.
We again vary b G , but now in the smooth (η, ξ)-system. The plot in Figure 9b displays the existence of a stable (spiral) equilibrium for b G larger than roughly 1.535. For b G < 1.535, there is a stable limit cycle created through a supercritical Hopf bifurcation, with the maximum (red) and minimum (blue) values of the norm on the limit cycle plotted as a function of b G . The existence of a Hopf bifurcation has played an important role in other models of the glacial cycles; see, for example, [10,46].
Returning to the discontinuous model, one can also vary the diffusion constant parameters. For example, the choice of D G determines the extent of the glacial state ice cap: reducing D G leads to larger ice caps (recall Figure 1). Varying D G then leads to limit cycles with larger amplitude and period ( Figure 10), akin in spirit to the behavior associated through the MPT. Welander's mixed ocean layer flip-flop model [28] was placed in the nonsmooth setting and analyzed in [47]. The periodic orbit in the Filippov analysis of Welander's equations arises in a nonsmooth Hopf bifurcation through a fused focus [27] as a parameter is varied. We note that the corresponding bifurcation in our model is entirely analogous.
Recall that there is a repelling sliding region between the invisible fold singular points u G and u I on Σ, assuming the time constant for the evolution of the ice edge ξ is sufficiently small. As is increased, u G and u I merge in a double invisible tangency point (what Filippov termed a fused focus) and then pass each other, with an attracting sliding region now on Σ between u I and u G and no periodic orbit. This is the same type of bifurcation exhibited by Welander's model and presented in [47].
While simple in nature, the model presented here exhibits relatively rich dynamical behavior.

Forcing with Milankovitch Cycles
The eccentricity e of the Earth's orbit, which varies with a period of roughly 100,000 years, has an effect on the amount of insolation reaching the planet and thus on climate on geologic time scales. The global annual average insolation as a function of e is where Q 0 is the insolation value corresponding to e = 0 [11]. The magnitude of the change in Q due to eccentricity is roughly 0.6 W/m 2 . This change in incoming solar radiation by itself is too weak to account for the glacial-interglacial global temperature changes [5].
The obliquity β of the Earth's spin axis plays an important role in our climate over geologic time. The obliquity varies with a period of about 41,000 years, ranging between roughly 22 • and 24.5 • . For smaller β-values, less insolation is received in the polar regions, leading to an increased temperature gradient. This scenario is then favorable for the growth of a polar ice cap. Conversely, greater insolation is received in the polar regions for larger β-values, leading to a decreased temperature gradient. This scenario is thus favorable for the potential retreat of a large ice sheet. The roles of eccentricity and obliquity-and precession, if seasons are considered-in influencing the Earth's glacial-interglacial history have been well-studied (see, for example, [17,23,48]).
Recall that the insolation latitudinal distribution function (2) depends upon β. In [36], a mechanism to compute s(y, β) as a polynomial in both y and β, to any desired degree of accuracy, is presented. As discussed in [36], the second degree expansion s(y, β) = 1 − 5 8 P 2 (cos β)P 2 (y) (30) well approximates equation (2) for planets with obliquity values similar to Earth's. We incorporate (29) and (30) into the Filippov system (23), using values of e and β computed by Laskar [22]. A difficulty that arises immediately is the proximity of the (interglacial) small, saddle ice cap equilibrium to the stable node (η I , ξ I ) when D = 0.394 (Figure 2b). In every simulation run with Milankovitch forcing, the φ I -trajectory shoots poleward of the stable manifold (25) of the small ice cap saddle and heads to the boundary η = 1. We thus use the reduced value D G = 0.38, for which there is no small ice cap saddle equilibrium (recall Figure 3), for this simulation.
Results of a representative solution are shown in Figures 11 and 12. Note the trajectory continues to pass through the discontinuity boundary when in both glacial and interglacial modes, thereby producing the glacial cycles. The top image in Figure 12 shows the behavior of η and ξ over the past 1 million years, that is, through the past 10 model glacial cycles. The effect of Milankovitch forcing is clearly evident in the plot of the albedo line η. The lower plot in Figure 12 displays the meridional flux at the albedo line (η) (27) (red, scaled and translated), as well as the obliquity β (blue). We focus on the obliquity as the magnitude of the change in insolation due to time-varying eccentricity is too small to account for the glacial-interglacial global climate changes, as mentioned above. We also note the emphasis here is placed on η rather than ξ as it is the albedo line (or snow line) that principally drives the model dynamics, due to its role in the competing feedbacks discussed in the introduction.
Note that maximum and minimum β-values correspond to minimum and maximum values, respectively, of (η). This is consistent with the fact that smaller β-values lead to larger temperature gradients and increased meridional transport, while larger β-values lead to reduced temperature gradients and diminished meridional transport.
Huybers argued in [17] that high obliquity is essential for the initiation of a glacial retreat. Indeed, the dominant signal seen in the climate data prior to the MPT is obliquity, as mentioned previously. Huybers also noted that after the MPT the climate data indicate that at times two maximum β-values are skipped before initiation of a deglaciation. He further posited that perhaps the 100 kyr cycles after the MPT are the average of 80 kyr (two β-cycles per glacial cycle) and 120 kyr (three β-cycles per glacial cycle).
Huyber's proposition is supported by the model output. As can be seen in the lower plot in Figure 12, one half of the model glacial cycles skip two maximum β-values. Correspondingly, the flux at the albedo line passes through two minimum values during one half of the ten glacial cycles in Figure 12, prior to decreasing to the the flux value at η during an interglacial. The conceptual model introduced here illustrates Huyber's concept quite well, while bringing the poleward flux at η in response to obliquity forcing into play in a novel way.

Discussion and Future Work
We have introduced a new conceptual model of the glacial cycles, with an energy balance annual and zonal mean surface temperature equation with diffusive heat transport (1) the starting point. We recast the coupled dynamic ice line introduced in [15] as the albedo line, entirely consistent with Budyko's piecewise-constant albedo function (3). As in [20], the new ice edge variable ξ is incorporated to allow for consideration of ablation and accumulation zones on the ice sheet.
The present work differs from [20] in several important regards. The meridional transport term used in [20] is Budyko's relaxation to the mean term (5). The associated zonal average surface temperature equation again presents an infinite dimensional dynamical system, one which is approximated in a manner far different from that used here in Section 3 [12]. Lacking a diffusive process, it is the critical temperature T c that is varied between glacial and interglacial states in [20]. One motivation for this choice was provided by [8], in which the value of T c was assumed to depend upon changes in the deep ocean temperature.
The model presented here was strongly motivated by work of Raymo and Nisancioglu [23], in which competing feedbacks associated with meridional transport and related to ice-albedo feedback are discussed in depth. It is natural to connect these feedbacks with poleward flux, as we have done here, which then leads to consideration of a diffusive transport term.
We note the use of diffusive heat transport smooths the albedo function discontinuity in the sense the temperature profile is continuous ( Figure 5; this is not the case when using transport term (5)). The discontinuity in our model vector field arises from the use of a larger ablation rate during interglacial retreats, and a switching mechanism motivated by mass balance principles as in [20].
The model focuses on the interplay of feedbacks associated with ice-albedo feedback, namely, the temperature gradient (positive) and transport efficiency (negative) feedbacks, similar in spirit to considerations in [23]. We assume a large ice sheet results in less efficient heat transport, and hence we use a smaller "glacial" parameter D G than "interglacial" parameter D I .
An ad hoc method is used to analyze the resulting nonsmooth geometric singular perturbation problem. A nonsmooth attracting limit cycle is shown to exist, provided the ice edge moves slowly enough. The model oscillations continue when forced by time-varying obliquity and eccentricity values. We note the competing feedbacks mentioned above serve to determine the maximal ice sheet size.
Many questions of a mathematical nature are raised by this flip-flop glacial cycle model. It would be of interest to extend the analysis presented here to the boundary of B, particularly the segments η = 1 and ξ = 1. The Filippov flow would have to be extended to include this portion of state space, perhaps along the lines of work by Barry et al. in [16].
Interesting questions concerning bifurcations would then arise for the extended flow and a state space that includes the boundary of B. What becomes of the saddle-node bifurcation indicated in Figure 3 with dynamics defined on the line η = 1? In addition, an investigation of model dynamics in the case of the small unstable ice cap scenario (Figure 2b) would then be of interest when subjected to Milankovitch forcing, which can push trajectories poleward toward η = 1 as discussed previously.
Another intriguing problem concerns making, say, the stable equilibrium Q G visible and placing it near or on the discontinuity boundary Σ, and then determining the effect of Milankovitch forcing in this situation. Might Milankovitch forcing push the system to tip over to the interglacial mode, and how might any such tipping depend upon the time constant for the ice sheet edge?
We also note the work presented here is the first part of a larger project aimed at investigating the Pliocene-Pleistocene Transition (PPT) via our mathematical model. The PPT refers to the development of perennial sea and land ice in the Northern Hemisphere, of an eventual extent similar to that which we have today. The study [49] used ocean cores drilled at the present Arctic Ocean summer ice margin to reconstruct the local climate during this transition. The following findings come from this work.
Roughly 5 Mya (million years ago) the Arctic Ocean was ice free at the drilling site, or covered by first-year winter ice. Sea ice expanded from the central Arctic Ocean for the first time about 4 Mya, at a time when atmospheric CO 2 levels were similar to today's and the global average temperature was 1-2 • C warmer than at present. This early Pliocene climate is often viewed as an analog of a future warmer Earth.
Not until about 2.6 Mya, during the PPT, did Arctic Sea ice expand to its modern winter limits. The glacial cycle model presented here concerns the Earth's climate from this point to the present.
A long-term decrease in atmospheric CO 2 concentration accompanied the amplification of the polar cooling described above through the PPT. The parameter A in (1) is often used as a proxy for the concentration of greenhouses gases such as atmospheric CO 2 : larger CO 2 values lead to reduced OLR (a decrease in A), and smaller CO 2 concentrations lead to increased OLR (an increase in A).
Our model can be enhanced by taking A to be a dynamic variable, similar to the work of Barry et al. in [16]. Consideration of an ice-free state would entail an analysis of system (23) on the boundary η = 1 or ξ = 1. Either scenario might be associated with a decrease in CO 2 , that is, an increase in A, as atmospheric CO 2 is drawn down by a silicate weathering process [50] enhanced by the lack of ice cover. In future work, the idealized mathematical model presented here will be extended in an effort to gain insight into the interactions of the large scale mechanisms and feedbacks in play through the PPT.