Nonsmooth frameworks for an extended Budyko model

In latitude-dependent energy balance models, ice-free and ice-covered conditions form physical boundaries of the system. With carbon dioxide treated as a bifurcation parameter, the resulting bifurcation diagram is nonsmooth with curves of equilibria and boundaries forming corners at points of intersection. Over long time scales, atmospheric carbon dioxide varies dynamically and the nonsmooth diagram becomes a set of quasi-equilibria. {However, when introducing carbon dynamics, care must be taken with the physical boundaries and appropriate boundary motion specified. In this article, we extend an energy balance model to include slowly varying carbon dioxide and develop a nonsmooth framework based on physically relevant boundary dynamics. Within this framework, we prove existence and uniqueness of solutions, as well as invariance of the region of phase space bounded by ice-free and ice-covered states.


1.
Introduction. There is significant evidence that large glaciation events took place during the Proterozoic era (2500-540 million years ago). In particular, this evidence points to the existence of glacial formations at low latitudes, see the review articles [18,28] and the references therein. One theory on the exodus from such an extreme climate was put forward by Joseph Kirschvink [22], who advocated that there was accumulation of greenhouse gases in the atmosphere, e.g. CO 2 . His theory purports that during a large glaciation chemical weathering processes would be shut down, thus eliminating a CO 2 sink. Moreover, volcanic activity would continue during the glaciated state. The combination of these effects would slowly lead to enough build-up of atmospheric carbon dioxide to warm the planet and start the melting of the glaciers. Once a melt began, a deglaciation would follow rapidly due to ice-albedo feedback.
There has been a wealth of modelling work on "snowball" events ranging from computationally intensive global circulation models (GCMs) to low dimensional conceptual climate models (CCMs). In 1969, Mikhail Budyko and William Sellers independently proposed energy balance models (EBMs); these were CCMs capturing the evolution of the temperature profile of an idealized Earth [5,30]. Many others, for example [7,1,27], have followed in the footsteps of Budyko and Sellers and used similar conceptual models capable of exhibiting snowball events. The low dimensionality of CCMs allows for a dynamical systems analysis, and hence a deeper investigation into some of the key feedbacks such as greenhouse gas and the ice-albedo effect.
From the point of view of dynamical systems, many of these early works share a similar theme by focusing on a bifurcation analysis with respect to the radiative forcing parameter, one that depends on changes in atmospheric CO 2 and other greenhouse gas levels. The reader may find figures similar to those displayed in Figure 1 in earlier works [1,18,29]. These figures illustrate the glaciation state of an Earth with symmetric ice caps in terms of atmospheric CO 2 i.e. ice line latitude at 90 o means Earth is ice free while 0 o means Earth is fully glaciated. In both figures, the effect of the radiative forcing due to CO 2 and other greenhouse gases is treated statically as a parameter in a simple bifurcation analysis. Mathematically, the bifurcation analysis already extends beyond the realm of smooth dynamical systems. On one hand, the bifurcation diagrams are obtained conventionally, with the bold curves indicating the stable branches and the dotted curves showing the unstable ones. On the other hand, the extreme states (ice-free and ice-covered) are not true equilibria of the system, but treated as such in the literature as evidenced by the labeling of "ice-free branch" and "ice-covered branch". Indeed, the extreme states of the ice line are special because they serve as physical boundaries for the dynamics: the ice line cannot extend beyond the pole nor, because of the north-south symmetry assumed in the models, can it move downward of the equator. Therefore, mathematical models of Snowball Earth must reflect this physical imposition and be treated with a nonsmooth systems perspective.
Since glacial extent varies over time and dynamic processes such as chemical weathering affect the level of atmospheric CO 2 , the bifurcation diagrams in Figure  1 can be viewed as phase planes with dynamic variables consisting of the glacier extent (ice line) and radiative forcing due to greenhouse gases. In particular, if a global glaciation did occur, and Kirschvink's argument about the accumulation of atmospheric CO 2 held, then we should expect orbits of this dynamical system to traverse the extreme ice latitudes, i.e. the equator and the pole. Following Kirschvink's idea further, we treat greenhouse gas effects on the energy balance by incorporating a slow CO 2 variable and treat the ice latitude as a relatively fast variable. The goal of the present article is to analyze such an interplay using a nonsmooth systems framework, thus providing mathematical support for Kirschvink's hypothesis. In this work, we treat the bifurcation diagram (similar to those in Figure 1) as a set of quasi-steady states of a slow-fast system and specify boundary motion based on Kirschvink's hypothesis and the long-term carbon cycle.
We present the topics as follows. In the next section, we motivate the model of interest. In Section 3, we propose an extension of the model that is amenable to a nonsmooth dynamical systems treatment. In particular, we show that a nonsmooth value, DA = A 0 − A, which we can think of as the radiative forcing due to increased optical thickness of the atmosphere (due to increased greenhouse gases, for example) relative to the present.
[23] Roughly following Roe and Baker [2010], we can calculate the stability of the equilibria by first noting from equation (7) that so that which leads to If xs DA ð Þ > 0 on a steady state solution then reducing the radiative forcing causes the ice latitude to move equatorward, and vice versa, so the state is stable. Similarly when xs DA ð Þ < 0 the solution is unstable and when xs DA ð Þ becomes infinite there is a bifurcation. This is an example of the "slope-stability" theorem, a more general consideration of which can be found in Cahalan and North [1979]. Using equation (10) it is easy to show that the solution line is stable for xs greater than the bifurcation at 0 < xs < 1 and unstable for x s less than this bifurcation (use @p @xs = S(x s )(a 1 − a2(xs)) from equation (5) and note that s2 is negative). We show a bifurcation diagram of the standard Budyko-Sellers model in Figure 10, generated using the solution and stability criteria outlined in this section.

Modifying the Budyko-Sellers Model to Produce the Jormungand State
[24] The Jormungand global climate state depends on the difference between the albedo of snow covered and bare sea ice, which is not incorporated into the standard Budyko-Sellers model. To include this effect we let a2 take a snow covered value (a2 s ) poleward of a transition latitude and a bare sea ice value (a 2 i ) equatorward of this transition latitude. Mathematically, we take where xi is the sine of the latitude of the transition from bare to snow covered sea ice, and Dx i is the width (in sine latitude) of the transition region. In a global climate model xi and Dxi are set by atmospheric dynamics, but we must impose xi and Dx i in the Budyko-Sellers model. We can still use equation (7) to solve the model as long as we recalculate as and ap(xs) using equation (11). [25] To calculate stability of the states, we calculate xs DA ð Þ using the method described in section 3.1. We find a relation similar to equation (10) Note that in the standard Budyko-Sellers model @s @xs = 0, but in the modified Budyko-Sellers model @s @xs is nonzero near the ice-snow transition.
[26] Our parameter choices for the modified Budyko-Sellers model are motivated by results from CAM, but ultimately selected so that the modified Budyko-Sellers model bifurcation diagram is similar to that of CAM. The fact that it is possible to obtain a similar bifurcation diagram at all with our simple modification to the model (equation (11)) suggests that our explanation for the Jormungand state is sound. Furthermore, we will subject the model to a sensitivity study (section 3.3) that shows that the behavior we describe is fairly robust, and suggests the qualitative control on the Jormungand state of the physical processes that model parameters represent. [27] We take xi = 0.35 based on the CAM simulations (section 2) and Dx i = 0.04. We choose albedo and ice-snow transition parameters of the modified Budyko-Sellers model to correspond roughly with those produced by CAM (section 2). We also take C B = 1.5, as opposed to C B = 2.4 in the work by Budyko [1969] and C B = 2.5 in Figure 10, because meridional heat transport in CAM is relatively inefficient in the Jormungand state (poleward heat transport across 45°of ∼2.25 PW (10 15 W) in the Jormungand state as opposed to ∼4.75 PW in the warm state). This is the result of both a strong inversion in the mid and high latitudes, which reduces eddy activity, and extremely cold temperatures in the mid and high latitudes, which reduces moist energy transport. Finally, we let Ts = 0°C since this value is more appropriate if the ice Assuming an albedo runaway did occur, the climate would be dominated by the dry atmosphere and the low heat capacity of the solid surface (Walker, 2001). It would be more like Mars (Leovy, 2001) than Earth as we know it, except that the greater atmospheric pressure would allow surface meltwater to exist. Diurnal and seasonal temperature oscillations would be strongly amplified at all latitudes because of weak lateral heat transfer and extreme 'continentality' (Walker, 2001). Despite mean annual temperatures well below freezing everywhere, afternoon temperatures in the summer hemisphere would reach the melting point (Walker, 2001). Evaporation of transient melt water would contribute, along with sublimation, to maintain low levels of atmospheric water vapour, and glaciers would feed on daily updrafts of this moisture (Walker, 2001). The global mean thickness of sea ice depends strongly on sea-ice albedo ( 1.4 km for albedo 0.6) and meridional variability is a complex function of solar incidence, greenhouse forcing (see below), zonal albedo, ablation or precipitation, and equatorward flowage of warm basal ice (Warren et al., 2002).
Climate physicists originally assumed that no ice-albedo catastrophe ever actually occurred because it would be permanent: the high planetary albedo would be irreversible. A saviour exists, however, and Kirschvink (1992) identified it as the buildup of an intense atmospheric CO 2 greenhouse through the action of plate tectonics in driving the long-term carbon cycle (Walker et al., 1981; Caldeira and Kasting, 1992; Kirschvink, 1992). On a snowball Earth, volcanoes would continue to pump CO2 into the atmosphere (and ocean), but the sinks for CO2 -silicate weathering and photosynthesis -would be largely eliminated (Kirschvink, 1992). Even if CO2 ice precipitated at the poles in winter, it would likely sublimate away again in summer (Walker, 2001). CO2 levels would inexorably rise and surface temperatures would follow, most rapidly at first and more slowly later on (Fig. 7) due to the nonlinear relation between CO2 concentration and the resultant greenhouse forcing (Caldeira and Kasting, 1992). With rising surface temperatures, sea ice thins but ground ice sheets expand in some areas due to a stronger hydrological cycle. If CO2 outgassing rates were broadly similar to today (a reasonable assumption for 600-700 Ma), then the time needed to build up the estimated 0.12 bar CO2 required to begin permanent melting at the equator, assuming a planetary albedo of 0.6, would be a few million years ( Effect of a 30% reduction in meridional heat transport is shown, as is the estimated solar flux at 600 Ma. Of three possible stable points for Es ¼ 1.0, the Earth actually lies on the partially ice-covered branch at point 1. An instability due to icealbedo feedback drives any ice-line latitude < 30°onto the ice-covered branch. A pCO2 » 0.12 bar is required for deglaciation of an ice-covered Earth, assuming the planetary albedo is 0.6 and Es ¼ 1.0 (Caldeira and Kasting, 1992). The snowball Earth hypothesis is qualitatively predicated on these findings and infers a hysteresis in pCO2 (and consequently surface temperature) following the circuit labelled 1-7. Starting from point 1, lowering of pCO2 causes ice lines to migrate stably to point 2, whereupon runaway ice-albedo feedback drives ice lines to the ice-covered branch at point 3. Normal volcanic outgassing over millions of years increases pCO2 to point 4, initiating deglaciation. Reverse ice-albedo feedback then drives ice lines rapidly to the ice-free branch at point 5, where high pCO2 combined with low planetary albedo creates a transient ultra-greenhouse. Enhanced silicate weathering causes lowering of pCO2 to point 6, whereupon polar ice caps reform and ice lines return to the partially ice-covered branch at point 7. In the 1960s, Budyko was concerned with the small icecap instability, which predicts a possible switch to the ice-free branch (e.g. disappearance of Arctic sea ice) due to anthropogenic global warming.   ice line model based on a latitude-dependent EBM coupled with a simple equation for greenhouse gas evolution has unique forward-time solutions. We also show that a similar result can be obtained by embedding the system in the plane and utilizing Filippov's theory for differential equations with discontinuous right-hand sides in Section 3.2. We end Section 3 with an analysis of the system dynamics and conclude the paper with a discussion in Section 4.
2. The equations of motion.
The main idea is that the change in temperature is proportional to the imbalance in the energy received by the planet. The amount of short wave radiation entering the atmosphere is given by Qs(y)(1 − α(η, y)), where Q is the total solar radiation (treated as constant), s(y) is the distribution of the solar radiation, and α(η, y) is the albedo at latitude y given that the ice line is at latitude η. The outgoing longwave radiation is the term A + BT (y, t); it turns out that the highly complex nature of greenhouse gas effects on the Earth's atmosphere can be better approximated by a linear function of surface temperature than by the Stephan-Boltzmann law for blackbody radiation (σT 4 ), see the discussion in [17]. The parameter A is particularly interesting here, because it is related to greenhouse gas effects on the climate system, which we will describe further in Section 2.2. The transport term C( 1 0 T (y, t)dy − T ) redistributes heat by a relaxation process to the global average temperature. Note that Budyko's equation assumes symmetry of the hemispheres and so one typically works only with the northern hemisphere, i.e. y ∈ [0, 1]. A more detailed discussion of this model can be found in [33]. Table 2 lists the polynomial expressions of some key functions in this model using standard parameter values.
The evolution of the ice-water boundary, or what is often called the ice line, η, provides a positive feedback to the system through the albedo function α(η, y). In [33] and [27], a critical ice line annual average temperature, T c , is specified. Above this temperature ice melts, causing the ice line to retreat toward the pole, and below it ice forms, allowing glaciers to advance toward the equator. One way to model this is described in [36], where the augmented ice line equation governing η is written as In [26], McGehee and Widiasih imposed equatorial symmetry on the set of solutions to (1)-(2) by expanding the temperature profile T (y, t) in even Legendre polynomials in the spatial variable y. They found that the space of even Legendre polynomials contains all equilibrium solutions of the model and is invariant for the system when using the quadratic approximation to the insolation distribution s(y) that was proposed by North [27]. By using only the first two terms in the expansion, the authors reduced the infinite-dimensional system (1)-(2) to one consisting of five ordinary differential equations, four representing the temperature profile and one for the ice line. They also found that the five-dimensional temperature profile-ice line system could be further reduced to a single equation, due to the existence of a one-dimensional invariant manifold parametrized by the ice line. The equation with where α(η) = 1 0 s(y)α(η, y)dy. The parameter A appearing in (3) is indeed the greenhouse gas parameter in [26], and is playing a similar role as that in the bifurcation diagram shown in Figure 1a.
In the same work, they examined the consistency of the coupled temperature profile-ice line system (1)-(2) with physical interpretation. The orbital (i.e. Milankovitch) forcings affect the amount of sunlight that the planet receives, hence, its temperature. However, the planet's response to these forcings is delayed. The Budyko equation has been shown to be a good approximation to Earth's temperature with the assumption of a 2500 year delay between the orbital forcings and the climate response [25]. Building on this assumption, the authors in [26] show that the parameter ρ must be small, and approximate it to be about 0.001.
In what follows, we couple an equation for atmospheric greenhouse gas evolution through the (former) parameter A. We couple this equation to the reduced model of McGehee and Widiasih [26] and use this as a proof of concept for the nonsmooth questions of interest. Numerical simulations produce similar dynamics in higherdimensional versions of the model, but making a rigorous reduction of the infinitedimensional system to the nonsmooth system studied in this article would require new center manifold results for discontinuous vector fields. We leave this for future work. From here on, we highlight the explicit dependence of the equation for η on A by writing where h(A, η) = h 0 (η; A) and h(A, η) is thought of as a real valued function over the plane R × R.  when y = η α 2 when y > η Table 2. Functions as in [26].

2.2.
Incorporating greenhouse gases. We now make an argument for a very simple form of an equation for greenhouse gas evolution. In the Budyko and Sellers models [5,30], the parameter A plays the important role of reradiation constant. The Earth absorbs shortwave radiation from the sun, and some of this is reradiated in the form of longwave radiation. The current value of A is measured using satellite data to be approximately 202 W/m 2 [17,33].
However, throughout the span of millions of years the longwave radiation parameter A is not constant, as the amount of heat reradiated to space depends crucially on greenhouse gases, especially carbon dioxide. This relationship is highly complex, but there have been some attempts at creating simplified models. One example comes from [7], where A was expressed as a polynomial expansion in the logarithm of atmospheric carbon dioxide for CO 2 near modern values measured in parts per million. Such an expansion is unlikely to directly apply over the time period associated with entrance into and exit from Snowball Earth, where carbon dioxide would have varied over many orders of magnitude. Following [1] (see also Figure 1a), we choose to vary A as a proxy for CO 2 rather than attempting to build a model for their direct relationship. Intuitively, adding carbon dioxide to the atmosphere decreases its emissivity, allowing less energy to escape into space. Therefore, outgoing longwave radiation should vary inversely with CO 2 .
Since the land masses were concentrated in middle and low latitudes prior to the global glaciation period, Kirschvink postulated that the stage was set for an ice-covered Earth. Then 'On a snowball Earth, volcanoes would continue to pump CO 2 into the atmosphere (and ocean), but the sinks for CO2 silicate weathering and photosynthesis would be largely eliminated' [18,22]. In more recent work, Hogg [19] put forth an elementary model for the evolution of greenhouse gases consistent with Kirschvink's theory. In short, he argued that the main sources for atmospheric CO 2 were due to an averaged volcanism rate and ocean outgassing. The main carbon dioxide sink was said to be due to the weathering of silicate rocks. For our purposes, it will be enough to work with volcanism and weathering. Let V denote the rate of volcanism and W the weathering rate, where the latter is assumed to depend on the location of the iceline, or more specifically the amount of available land or rock to be weathered. Putting this together with the assumption that the reradiation variable A varies inversely with CO 2 gives where δ > 0 and η c = V W , the ratio of volcanism to weathering. Let g(A, η) := δ(η − η c ).
Next, allowing A to vary in the ice line equation (4), we now have a system of equations for A and η: To understand the timescale of A, we refer to the elucidation of the snowball scenario by Hoffman and Schrag in which they argue that weathering took place at a much slower rate than the ice-albedo feedback (see Fig. 7 on page 137 [18]).
The idea of packaging the essence of the long term carbon cycle into one simple equation is not novel, and is consistent with earlier findings [13,19,24]. The novelty comes from connecting such an equation to an energy balance model and analyzing the dynamics of the coupled system. Indeed, the parameter A ex on page 644 [24] is the variable η in equation (6) and it represents the effective area of exposure of fresh minerals. The parameter η c can therefore be thought of as a critical area. The coupling of the ice line η and the greenhouse gas variable A follows naturally.
In its current form, the model does not restrict dynamics to the physical region 0 ≤ η ≤ 1. There are orbits that exit this interval on either end, and we must therefore create reasonable assumptions for projection of this motion onto the boundaries. Furthermore, η should evolve along the physical boundaries, as suggested by the bifurcation diagrams in Figure 1. In the next section, we develop a nonsmooth version of the model that respects these boundaries and prove existence and uniqueness of solutions. We use this as a proof of concept for how to treat physical boundaries beyond which a vector field is undefined.
3. A nonsmooth system for a glaciated planet. In this section, we introduce suitable constraints on the ice dynamics so that trajectories with initial conditions in the physical region remain in this region for all time. Moreover, we obtain physically reasonable motion on the boundaries.
As mentioned previously, the pole and the equator define physical constraints of the ice line; it must stay in the unit interval [0, 1]. However, the model in its current form does not take this into account. A natural choice is to setη = 0 when η reaches a physical boundary. In addition to this, one needs to account for the stability/instability of these states, thus allowing orbits to exit the boundary when it becomes unstable. More specifically, it would be natural to define the evolution of (A, η) via Ȧ = g(A, η) where f (A, η) = 0 {η = 0 and h(A, η) < 0} or {η = 1 and h(A, η) > 0} h(A, η) otherwise (8) and g, h are as in Table 2. The second equation forces η to stop at the physical boundary exactly when it is about to cross that boundary. When this happens, η should remain constant while A continues to evolve. As A evolves, so does the value of h and hence η should be allowed to re-enter the interval (0, 1) as soon as h becomes positive on the lower boundary, η = 0, or negative on the upper boundary, η = 1. This is analogy with Kirschvink's hypothesis, as we expect orbits that enter the snowball state to recover when enough carbon dioxide has built up in the atmosphere. However, the right-hand side of (7) is discontinuous at the boundaries and so solutions will not necessarily be differentiable there. Thus, we define a solution to an initial value problem for (7) with initial condition (A(0), η(0)) = (A 0 , η 0 ) to be an absolutely continuous function X(t) = (A(t), η(t)) that satisfies the integral equation X(t) = X 0 + t 0 (g(A(s), η(s)), f (A(s), η(s))) ds (9) with f and g as above. Note that since h is smooth, the integral equation is equivalent to the system (Ȧ,η) = (g, h) inside the strip 0 < η < 1.
Remark 3.1. The discontinuous vector field arising from (7) defines a projection rule for orbits that reach the boundary: when the vector field points out of the strip R × [0, 1], the forward evolution is defined by projection of the vector field onto the boundary, and when the vector field points into the strip, the system is left to evolve without interruption.
Because there is no physically reasonable way to define the vector field outside of the interval 0 ≤ η ≤ 1, a system with such a rule does not fit readily into an existing analytical framework. Indeed, common frameworks for nonsmooth systems assume that the vector field is defined in a neighborhood of the curve along which it is discontinuous. Thus it is not clear that classical results from dynamical systems apply in this setting. Fundamentally, the existence and uniqueness of solutions needs to be verified. In what follows, we prove that the system with the rule specified above has unique forward-time solutions and that the physical region is forward invariant. We also outline an alternative approach that involves defining a nonphysical "confinement" vector field outside the physical region and serves to handle and correct the behavior of numerically generated solutions that may overshoot the boundary. This second approach also guarantees existence, uniqueness, and forward invariance. Further, its solutions trace out the same curves in phase space as the previous approach, though they differ in their smoothness.
3.1. First framework: the projection rule. In this section, we show that the projection rule results in existence and uniqueness of solutions.
Theorem 3.2. Let 0 < η c < 1 and t > 0. Then an initial value problem for (7) with initial condition (A 0 , η 0 ) and 0 ≤ η 0 ≤ 1 has a unique forward-time solution of the form (9). The time derivative of the solution has at most a finite number of discontinuities and satisfies (7) in between them. Furthermore, the strip R × [0, 1] is forward invariant.
Proof. First, a solution is created in an ad hoc manner by concatenating pieces of solutions to the associated smooth systems on each subregion of the phase space. Since f and g are smooth on the open set {0 < η < 1}, we focus only on cases where entrance into or exit from η = 0, 1 occurs. Without loss of generality, consider an initial condition (A 0 , η 0 ) with 0 < η 0 < 1 such that the smooth solution (A(t), η(t)) reaches the line η = 1 at some time t 1 < t, i.e. lim t→t − 1 η(t) = 1. From Table 2, we see that h(A, η) = ρ 1.5 (f (η) − A) where f is a cubic function of η. Thus, as can be seen in Figure 2, the set {h = 0} divides the strip into two regions; h > 0 on the left-most region and h < 0 on the right-most region. Define 1) > 0, then h being Lipschitz continuous implies that it must be positive in a neighborhood of (A 1 , 1). Therefore, we extend the solution according to (8) and (9) via . Since 0 < η c < 1, note that δ(1 − η c ) > 0 and hence there is a single orbit satisfying the smooth system of differential equations that passes the curve h = 0 from left to right and touches the boundary at the single point (A 2 , 1) forming a tangency. Moreover, this orbit is strictly contained in the interval 0 < η < 1 on either side of the tangency and also satisfies the integral equation (9) where f is replaced by h. Thus, taking an initial condition (A(t 2 ), η(t 2 )) = (A 2 , 1) and concatenating the forward evolution of this trajectory with the previous piece of the solution to the nonsmooth system provides a further extension via A(s) = A 2 + s t2 g(A(τ ), η(τ )) dτ η(s) = 1 + s t2 h(A(τ ), η(τ )) dτ, for all t 2 ≤ s ≤ t such that 0 < η(s) < 1. Suppose t 3 is the next time the created solution reaches the boundary. In the case that the orbit again reaches the ice-free state η = 1, we extend the solution exactly as above and this definition leads to another exit from this boundary at the point (A 2 , 1) as before. In this case we have created a closed orbit. Analogous extensions can be made when the ice-covered boundary η = 0 is reached. Any solution created by this method of concatenation that reaches the same boundary twice lies on a closed orbit. This in turn implies that the created solution can reach the boundary at most finitely many times on the interval [0, t], since it will spend the same amount of time there for each subsequent passage. This solution is absolutely continuous and differentiable everywhere except at the points where it enters one of the lines η = 0, 1. By definition, it satisfies (7) everywhere except at these points. Now suppose there is another absolutely continuous solution (Ā(t),η(t)) of (7) with (Ā(0),η(0)) = (A 0 , η 0 ). Due to smoothness of h and g, both trajectories must reach the boundary at the same time and hence they agree for all 0 ≤ s < t 1 . If they reach the boundary at the tangency point, then the integrands of (9) continue to be h and g and hence classical uniqueness arguments can be used until such time that the solutions reach another boundary. If they reach the boundary at (A 1 , 1) such that h(A 1 , 1) > 0, then both satisfy (10)- (11) and uniqueness persists until t 2 such that h(A(t 2 ), 1) = h(Ā(t 2 ), 1) = 0.
Indeed, since g and h are Lipschitz continuous, classical arguments invoking Grönwall's Inequality can be used to complete the uniqueness proof. In particular, for t 2 ≤ s ≤ t such that 0 < η(s) ≤ 1 one finds that where K 1 and K 2 are the Lipschitz constants of g and h. Similarly, uniqueness persists across any entrance into or exit from η = 0, 1, which is enough to guarantee (A(t), η(t)) = (Ā(t),η(t)). Forward invariance of the strip is also immediate from the rules specified in (8).
Note that Theorem 3.2 only guarantees forward-uniqueness of solutions of (7), and in general solutions are not unique in backward time; any solution that passes through a point on an attracting portion of the boundary (i.e. where h(A, η) > 0 on η = 1 or where h(A, η) < 0 on η = 0) could have have reached that point directly from the interval 0 < η < 1 or from further back along the boundary itself.
Remark 3.3. The solutions guaranteed by Theorem 3.2 are of the same type as those guaranteed by Carathéodory's Existence Theorem [9] (see also [16]). However, the existence theorem does not apply here because it requires that the right-hand sides of the differential equations be continuous in (A, η). In the above theorem, we have shown that the projection rule specified in (7) and Remark 3.1 has enough continuity to guarantee existence and uniqueness of solutions. More specifically, note that this rule guarantees upper semi-continuity of the vector field, which is also a key requirement for existence and uniqueness of Filippov's notion of solution [16], and which we explore in the next section. We will see that the solution given by the theorem above agrees with Filippov's solution, but the latter may have multi-valued derivative.
3.2. Alternative approach: Filippov framework. In this section, we introduce a Filippov systems perspective by way of defining a confinement vector field that points into the unit interval 0 < η < 1 from outside. This perspective is useful for numerical simulation because it handles the case when orbits overshoot the physical boundary. We wish to choose a confinement field that preserves the speed of motion along the boundary guaranteed by the projection rule (7) and obtain solutions which have the same forward evolution. We will see that Filippov's convention [16] preserves this speed if we work with the following piecewise-Lipschitz vector field on R × R:Ȧ having an initial value (A(0), η(0)) ∈ R × R.
In the above scenario, the boundaries η = 0, 1 become curves along which the vector field is discontinuous (so-called switching manifolds) and reasonable motion along these curves must be specified. Because we have chosen to work with Filippov's convention, we take all convex combinations of the limiting vectors from either side, i.e. linear combinations of the Lipschitz vector fields. A more general convention allowing for nonlinear combinations was discussed in [20], and taking such a perspective may require a different choice of external vector field to preserve the speed of boundary motion. We do not discuss the nonlinear case further here. Instead, we consider the differential inclusioṅ and let f (A, η) be the vector field defined by right-hand sides of (13)- (15). A Filippov solution to (13)-(15) will be defined as an absolutely continuous function X(t) = (A(t), η(t)) such thatẊ ∈ f (X) for almost all time t. We observe that the vector field defined by (13)- (14) points "into" the interval 0 < η < 1 except along the curve h(A, η) = 0 where it is parallel to η = 0, 1. Thus, according to Filippov's theory, there can only be sliding or crossing behavior along the boundary away from the two transverse intersections with h = 0. All crossing orbits must enter the interval 0 < η < 1. As we will see, the points where h = 0 signal the transition from (attracting) sliding to crossing. We now show the existence and uniqueness of Filippov solutions and the forward invariance of the strip. Proof. This is a direct application of Filippov's existence and uniqueness result, Theorem 14, p. 228, in [15]. Using similar notation, let f (z) denote the vectorvalued function [g(z), h(z)], the normal direction N is the vector [0, 1], and anything with a subscript N denotes the projection onto N . By design, the confinement vector fields above and below the strip R × [0, 1] point into the strip. We consider a small neighborhood G around a point z 0 away from the nullcline h = 0 on the line η = 1. Letting f + and f − denote the vector field above and below η = 1, we see that f + − f − is continuously differentiable, and the condition that either f + N < 0 or f − N > 0 is satisfied for all points in G. Then Filippov's theorem asserts the existence and uniqueness of a forward-time solution to system (13)- (15) as well as continuous dependence on initial conditions. Lastly, Remark 1, p. 229 [15] asserts that the result holds at the intersection of η = 1 with the h-nullcline, since at that point, f + = 0 = f − . A similar argument applies for initial conditions near η = 0.
Lastly, since any solution can only cross into or slide along the boundary, the strip is forward invariant.
Remark 3.5. There are many possible ways of extending system (6) beyond the strip R × [0, 1]. The choice of extension as done in system (13)-(15) assures the following points: 1. The strip R × [0, 1] is attracting. 2. The switching from sliding to crossing is entirely determined by zeroes of h(A, η). This point is especially important in the modeling of the physical system, since h signals the ice line's advance or retreat and h should be "off" for some positive amount of time at the extreme ice line locations (pole or equator) and should "turn back on" when the system has sufficiently reduced or increased greenhouse gases. 3. Filippov's theory asserts that at the attracting sliding segment on η = 1, the time derivative must be dη dt = −|h| + h 2 almost everywhere, while at η = 0, it is dη dt = |h| + h 2 almost everywhere. Hence, the effective speed at which orbits slide along the boundaries is zero and the motion is thus entirely determined by the governing function of the greenhouse gas effect, g(A, η). Again, this aspect is especially important in the modeling of the physical system, since at the extreme ice line location (pole or equator), the ice-albedo feedback shuts off and the greenhouse gas effect is the only driver of the system.

3.3.
Dynamics of the model. We now explore the dynamics of the system (7). We show that its boundary dynamics guarantee exit from the snowball scenario, and use the Filippov perspective to run simulations. We find that the system possesses a stable small ice cap regime as well as large amplitude periodic orbits that oscillate between ice-covered and ice-free states.
The η-nullcline, given by h(A, η) = 0, is cubic in η and linear in A (see Table  2). It has a single fold in the region 0 ≤ η ≤ 1 at η = η f ≈ 0.77, as can be seen in Figure 2. The A-nullcline is the horizontal line η = η c . If 0 < η c < 1, the system has a single fixed point. The eigenvalues associated with this fixed point are From this and Figure 2, we see that if the fixed point lies above the fold and in the physical region, i.e. η f < η c < 1, then it is stable, and it is unstable when 0 < η c < η f . This is reminiscent of the "slope-stability theorem" from the climate literature [6]. In that work, the authors related changes in slopes of equilibrium curves to changes in stability of equilibria for a globally averaged energy balance model. They found that a small ice cap or ice-free solution could be stable, a large ice cap was unstable, and a snowball state was again stable.
Remark 3.6. In the case that η c = 0, 1, one of the physical boundaries is entirely composed of equilibria, and stability is determined by the direction of the vector field near it. Physically, this degenerate scenario comes about when there is an absence of volcanism in the ice-covered state or a perfect balance between weathering and volcanism in the ice-free state. We do not discuss these cases further. Figure 2. The physical region of the phase space and possible fixed points of the system given by the η-nullcline, h = 0. The location of the equilibrium is determined by the critical effective area of exposed land 0 < η c < 1. Solid black portions of the curve represent stable equilibria while the dashed lines denote unstable equilibria. Solid black portions of the boundary are attractive sliding regions and dashed boundaries are crossing regions.
3.3.1. Sliding and escape from Snowball Earth. Recall that η c was defined to be a ratio of volcanism to weathering. We now prove that the system always recovers from the ice-covered regime provided this parameter is positive, i.e. there is some contribution from volcanic outgassing to atmospheric greenhouse gas content. The following result can be easily achieved following the proof of Theorem 3.2: Corollary 3.7. Suppose η c > 0. Then any orbit of the system (7) that enters the line η = 0 must exit in finite time.
Proof. Since dA dt η=0 < 0, any orbit that reaches the line η = 0 must eventually reach the unique point (A * , 0) such that h(A * , 0) = 0. Since h is negative to the right of this point and positive to the left, there is a single orbit that reaches the tangency point from the interval 0 < η < 1, passes through it from right to left, and then reenters the open interval. Due to the forward-time uniqueness guaranteed by Theorem 3.2, the sliding solution must follow this orbit and reenter the physical region.
Remark 3.8. Similarly, any orbit that enters an ice-free state must exit in finite time. The required condition is that there is not a perfect balance between volcanism and weathering, i.e. η c = 1.
3.3.2. Numerical simulations. As mentioned above, the Filippov framework allows for more robust numerical simulation of the model. Let (A c , η c ) be the location of the single fixed point. We find that small ice cap states corresponding to fixed points with η f < η c < 1 are possible attractors of the system, as are periodic orbits born from the Hopf bifurcation at η c = η f .
The fold of the η-nullcline is a canard point, and the mechanism that produces the large amplitude periodic orbit in Figure 3 is a canard explosion [11,2,8,4,12,23]. While mathematically interesting, canard orbits are unlikely observables for planar systems because they occur within a parameter range that is exponentially small in δ. However, they can be robust in higher dimensions which makes more complex oscillatory behavior possible, see e.g. [3,10,32,35]. With this in mind, we remark that an interesting avenue for future research involves studying an extended version of the system (13)- (14) with an additional variable (w in [26]) and foregoing their invariant manifold reduction.
In the simulations, trajectories that leave the upper boundary and re-enter the physical region are rapidly attracted to a slow manifold that is within O(δ) of the ηnullcline (the so-called critical manifold ) as guaranteed by Fenichel theory [14,21]. This is why the trajectories in the figures appear to have sharp corners at the exits from ice-covered and ice-free states, when in fact they are differentiable there and tangent to the relevant boundary.  η c = 0.6. The + symbol marks the initial condition and the horizontal long-dashed line is the A-nullcline. In (a), the orbit reaches the ice-free state and slides until it reaches the intersection of the folded curve with this boundary. It then enters the physical region and approaches the small ice cap equilibrium. In (b), the fixed point is unstable and the orbit oscillates between the ice-free and ice-covered boundaries. Parameters are as in Table 1 and δ = 0.01. Simulations were performed using Mathematica 9.

3.4.
Application to the Jormungand model. The previous occurrence of a complete snowball event is still a highly contested topic in geology. In the models discussed here, the only possible stable fixed point is a small ice cap. In this section, we present a variant of the system, as introduced by Abbot, Voigt, and Koll [1], and employ the nonsmooth framework as above. By modifying the albedo function α(η, y) so that it differentiates between the albedo of bare ice and snow-covered ice, the new system has an additional stable state that corresponds to a large ice cap. In [1], this was called the Jormungand state because it allows for a snake-like band of open ocean at the equator. Moreover, there are again attracting periodic orbits. These can be seen in Figure 4.
In this version of the model, the albedo is defined as follows: where α 2 (y) = 1 2 (α s + α i ) + 1 2 (α s − α i ) tanh M (y − 0.35). Here α w is the albedo of open water, α i is the albedo of bare sea ice, and α s is the albedo of snow-covered ice. The model assumes sea ice aquires a snow cover only for latitudes above y = 0.35. We modify h(A, η) by replacing α with α J : where we note that although α J (η, y) has a discontinuity across y = η, both α J (η, η) and α J (η) = 1 0 α(η, y)s(y)dy are smooth functions of η. For this system, we can immediately apply Theorem 3.2 and Corollary 3.7 to obtain the following result.
Corollary 3.9. The physical region is forward invariant with respect to the system (7) where h is replaced by h J . Solutions exist and are unique in forward time, and any trajectory that enters an ice-covered state must exit in finite time provided η c = 0.   Table 3 and δ = 0.01.

Parameters Value Units
0.45 dimensionless α s 0.8 dimensionless Table 3. Parameter values as in Table 1 unless specified above. Additional values taken from [1].
4. Discussion. In this article we have extended a class of energy balance models to include a greenhouse gas component. The resulting system was nonsmooth due to physical constraints on the extent of Earth's glaciation and did not immediately fit into an existing analytical framework. However, we defined two frameworks that produced the (same) expected motion, in agreement with physical arguments from the climate literature. Both approaches prevented the dynamics from exiting the physically relevant region; the first was based on a rule that projected the vector field onto the boundary when it was reached and the second defined a "confinement vector field" that pointed into the region. In both cases, we showed the existence and uniqueness of forward-time solutions. We further proved that the system always escapes from the ice-covered scenario, in agreement with Kirschvink's hypothesis about carbon dioxide accumulation in a snowball scenario due to a shut down of chemical weathering processes [22]. Using our model, we found that small ice-cap and large amplitude ice-cover oscillations were possible attractors of the system. We then applied our results to the case of a "Jormungand" world and showed that it is possible to obtain further oscillatory dynamics between no ice cover, large ice-cover, and full ice-cover states.
General mathematical questions brought about by this work have to do with building a general framework for such models. A first step toward dealing with physical boundaries might be a general theory for semiflows generated by smooth vector fields on manifolds with boundary. In addition, the stable portions of the physical boundaries in our model should be thought of as equilibria of the fast subsystem, i.e. as part of the critical manifold from geometric singular perturbation theory [21]. In fact, they are much more normally hyperbolic than the curve h = 0; they are reachable in finite time! However, the current theory cannot be directly applied to the full discontinuous system and we remark that developing an analogous Fenichel theory [14] for singularly perturbed discontinuous systems with sliding is an interesting avenue for future research.
Another interesting future direction is to study the explicit effect of the temperature (here we have considered it only through the equation for η based on the reduction by McGehee and Widiasih [26]) on the dynamics of the system. In a nonsmooth system, the addition of a stable dimension may destroy an existing stable periodic orbit [31]. Moreover, an additional dimension could result in more interesting glacial dynamics, such as mixed mode oscillations.
Finally, it is natural to ask how the shape of the η-nullcline changes with physical parameters and how this affects system dynamics, e.g. could the lower fold in Figure  (