ON THE BUDYKO-SELLERS ENERGY BALANCE CLIMATE MODEL WITH ICE LINE COUPLING

. Over 40 years ago, M. Budyko and W. Sellers independently introduced low-order climate models that continue to play an important role in the mathematical modeling of climate. Each model has one spatial variable, and each was introduced to investigate the role ice-albedo feedback plays in inﬂuencing surface temperature. This paper serves in part as a tutorial on the Budyko-Sellers model, with particular focus placed on the coupling of this model with an ice sheet that is allowed to respond to changes in temperature, as introduced in recent work by E. Widiasih. We review known results regarding the dynamics of this coupled model, with both continuous (“Sellers-type”) and discontinuous (“Budyko-type”) equations. We also introduce two new Budyko-type models that are highly eﬀective in modeling the extreme glacial events of the Neoproterozoic Era. We prove in each case the existence of a stable equilibrium solution for which the ice sheet edge rests in tropical latitudes. Mathematical tools used in the analysis include geometric singular perturbation theory and Filippov’s theory of diﬀerential inclusions.


1.
Introduction. The hierarchy of climate models begins on one end with conceptual energy balance and low-order models, progresses through intermediate complexity models, and culminates with ever-evolving Earth system models. Energy balance models (EBM) focus on major climate components, such as incoming solar radiation and planetary albedo, and their interactions, thereby providing a broadbrushstroke view of the dynamics of climate. Earth system models (ESM), on the other hand, require enormous computational resources as they seek to simulate the fully coupled atmosphere-ocean-land-cryosphere system. A typical ESM currently uses a horizontal resolution of roughly 100 km with 95 vertical levels [35], computating estimates of variables such as temperature, pressure, humidity and wind The most advanced comprehensive climate models effectively represent the current ability to simulate the climate system, and it is natural and appropriate to take the output of those models as the basis for predicting the future climate. However, it is the understanding of these responses-an understanding that depends on the presence of an emergent simplicity in the forced responsewhich provides a level of confidence that justifies advising policy-makers and the public to pay heed to these predictions.
We also note the analysis of EBM lends itself to a mathematical treatment, with techniques from dynamical systems theory finding increasing use.
Many low-order climate models arise via finite dimensional (ODE) approximations of the governing system of PDEs, the latter derived from the corresponding physics. The most well known example, perhaps, is that provided by E. Lorenz in [20], in which chaotic behavior was found in a general atmospheric circulation model consisting of three ODEs; other examples include [21], [34], [36], [38]. Low-order models are often used to investigate low-frequency variability in various climate components, such as in ocean-atmosphere interactions ( [3], [4], [32], [37]). The time scale for models such as these is typically months or years.
Another approach to low-order models arises from the direct analysis of reactiondiffusion PDE defined on the 2-sphere arising from the balance equation for the heat fluxes involved. One is typically interested in the evolution of the 10-year mean surface temperature in this class of models. Existence of global solutions and a trajectory attractor are often the focus of the analyses. An introduction to this PDE approach to EBM is provided in Part II of [6]; other illustrative works include [7], [8], [9], [16].
We focus here on the seminal zonal average EBM introduced independently by M. Budyko [5] and W. Sellers [33] in 1969 to investigate ice-albedo feedback in the climate system. This work serves in part as a tutorial centered on these models, with a particular emphasis placed on their coupling with an ice sheet edge capable of responding to changes in temperature, introduced by E. Widiasih [40]. In addition, we present new albedo functions appropriate for modeling the climate of the extensive glacial periods of the Neoproterozoic Era (roughly 700 million years ago). We prove the existence of a stable equilibrium solution for finite approximations of the infinite-dimensional coupled model having an ice edge in tropical latitudes. This so-called Jormungand climate state was proposed in [1] as an alternative to the Snowball Earth hypothesis, in which the Earth is posited to have been completely covered with ice as a consequence of ice-albedo feedback, as discussed below. The Jormungand state was found in simulations using an appropriately parametrized general circulation model (GCM) in [1].
In the following section we present the models of Budyko and Sellers, commenting where appropriate on their shortcomings. The introduction of a dynamic ice line, coupled with the one-dimensional EBM, is presented in Section 3. We survey recent results for this system of equations with both continuous ("Sellers-type") and discontinuous ("Budyko-type") albedo functions.
Section 4 introduces the Jormungand climate model, for which we review known results for the Sellers-type model. We then introduce two new Budyko-type albedo functions, proving in each case the existence of a stable equilibrium solution with ice line near the equator. The analysis requires tools from invariant manifold theory and the theory of non-smooth dynamical systems.
We discuss additional applications of the techniques introduced in Section 4 and we summarize this work in the concluding section.
2. One-dimensional energy balance models. The dimension of a low-order climate model typically refers to the number of spatial variables. The models of Budyko and Sellers, concerned with the mean annual surface temperature in latitude zones, are each one-dimensional in this regard. (From a dynamical systems perspective, however, the state space of the model PDE below is infinite-dimensional.) At a given latitude, EBM of this type were historically expressed as a true energy balance, written in the general form The "energy in" E in -term corresponds to the incoming solar radiation (or insolation), modulo the surface albedo, a measure of the extent to which insolation is reflected back into space. The E out -term models the outgoing longwave radiation (OLR) emitted by the Earth. Budyko chose an empirical approach, setting E out = A + BT , where T ( • C) is the mean annual surface temperature as a function of latitude, and A and B are constants determined by satellite data. Sellers chose a variant of the Stephan-Boltzmann Law for blackbody radiation, setting T in Kelvin. The lefthand side of (1) represents the zonal radiation balance. The models of Budyko and Sellers differ significantly in the E tran -term, used to model the transport of heat energy across latitudes. Sellers' model includes the net meridional fluxes of both water vapor (due to atmospheric currents) and sensible heat (due to atmospheric and ocean currents). Budyko, on the other hand, set E tran = C(T − T ), where T is the mean annual global surface temperature and C > 0 is an empirical constant. We note others have modeled meridional energy transport as a diffusive process, setting E tran = D∇ 2 T , where D is a constant typically chosen to fit the present climate. Finite Legendre polynomial expansions are often used in the analysis in this case, due to the form the diffusion term takes in spherical coordinates ( [15], [19], [26], [27], [28], Chapter 10 in [13], for example).
We focus here on Budyko's model due the simplicity of the OLR and energy transport terms. In addition, Budyko considered interactions between the atmosphere and hydrosphere, so one assumes the planetary surface consists of water and ice. While many early works assumed the energy was in balance as in (1), we consider the corresponding dynamical equation first introduced in [15]: The variable y is the sine of the latitude, chosen for convenience (the area of an infinitesimal latitudinal strip at y is proportional to dy). As the model assumes a symmetry across the equator we take y ∈ [0, 1], with y = 0 at the equator and y = 1 at the north pole. We refer to y as the "latitude" in all that follows, assuming no confusion will arise on the reader's part. The function T = T (y, t) ( • C) is the annual mean surface temperature on the circle of latitude at y, and R is the heat capacity of the Earth's surface (with units J/(m 2 • C)). The lefthand side of (2) represents the change in energy stored in the Earth's surface at y; units on each side of (2) are Watts per meter squared (W/m 2 ).
The mean annual global insolation is denoted by Q, with Q = 343 W/m 2 currently (generally, Q depends on the eccentricity of the Earth's orbit [24]). The distribution of insolation across latitude is provided by s(y), a function depending on the obliquity, or tilt, of the Earth's spin axis, and satisfying 1 0 s(y)dy = 1.
We note s(y) can be computed from astronomical principles [24], having explicit formulation where β denotes obliquity (see Figure 1). The function α(y) represents the albedo, so that Qs(y)(1 − α(y)) is the E in -term in (1), while the A + BT -and C(T − T )-terms, introduced above, represent the OLR and meridional heat transport, respectively. We also note T = T (t) = 1 0 T (y, t)dy.

Equilibrium solutions of Budyko's equation. The models of Budyko and
Sellers were introduced to investigate how solutions of energy balance (1) vary as parameters, such as Q, change. The corresponding time frame is thus hundreds of millennia. Working with (2), first note an equilibrium solution T * = T * (y) satisfies Integrating (4) from y = 0 to y = 1 yields from which it follows the mean global temperature at equilibrium is Substituting (6) into (4) and solving for T * yields Note T * (y) is completely determined once the albedo function is specified. Many early models (including that of Sellers) expressed the albedo as a function of y and T , typically in the form Here, T c is a critical temperature above which ice is assumed to melt and below which ice forms. Thus if T (y) > T c there is no ice and the albedo assumes a smaller value α w corresponding to open water, which is less reflective than snow and ice. For T (y) < T c there is ice at y and the albedo has the larger value α s of snow-covered ice.
With this albedo formulation, however, equilibrium solution (7) can take a variety of nonsensical forms, an example of which is plotted in Figure 2 (after [23]). If the explicit dependence on T is removed and α = α(y), then R. McGehee has shown equilibrium solution (7) is globally attracting under mild hypotheses (T (y) ∈ L 1 [0, 1], for example; see [23]).  One focus of Budyko and Sellers was the role played by positive ice-albedo feedback in determining climate. A cooling of the planet would lead to the expansion of an existing ice sheet, thereby increasing the albedo and reducing the temperature, leading to further ice sheet growth. The ice sheet would retreat were the planet to warm, reducing the albedo and increasing the temperature, so that the ice sheet would melt further. Might there exist a tipping point, the existence of an ice sheet sufficiently large, for example, so as to trigger a runaway Snowball Earth event? (Sellers was also interested in the potential consequences of a vanishing Arctic ice cap, as well as anthropogenic effects on climate, in [33].) Following Budyko [5] we assume the planet has an ice cap, with ice existing everywhere above a certain latitude y = η, and with no ice south of y = η. The edge of the ice sheet η is called the ice line. The albedo function, now parametrized by η, is defined by EBM with an inherent discontinuity at η (as in (9)) have been labeled "Budykotype" [7]. Sellers, on the other hand, assumed a continuous transition in the albedo function across the ice line, and EBM with this property are of "Sellers-type." We use this terminology in all that follows. Equilibrium solutions (7) with albedo function (9) are plotted in Figure 3 for a variety of η-values. Note that α = α(η) = α 1 η 0 s(y)dy + α 2 1 η s(y)dy is parametrized by η, so that an equilibrium solution exists with ice sheet edge at η for all η ∈ [0, 1], that is, equilibrium solutions (7) are parametrized by η as well. We write T * η (y) in place of T * (y) when we wish to emphasize the dependence of (7) on η.    Each equilibrium solution has a discontinuity at η; typically one sets the equilibrium temperature at the ice line to be the average One then requires the equilibrium temperature at η to equal the critical temperature. With parameter values appropriate for today's climate there are two such equilibrium solutions, one having a small ice cap and the other a much larger ice cap (see Figure 3). Budyko was particularly interested in how these two preferred solutions vary with changes in Q.
Setting T * (η) in (10) equal to T c gives a relationship between the position of the ice line at equilibrium and the parameter Q, which we plot in Figure 4. Note there are no equilibrium solutions with T * (η) = T c once Q decreases by roughly 5% from its current value. The implication is then that the planet plummets into a Snowball Earth episode for these Q-values (due to ice-albedo feedback), a possibility relevant to the Neoproterozoic Era, for example, when Q was 94% of its current value [29]. Notably, however, neither Budyko nor Sellers considered the stabilitytype (with respect to initial conditions) of these preferred equilibrium solutions with T * (η) = T c .
The analysis above provides one example of the utility of simple one-dimensional EBM in exploring the role played by major climate system components (such as ice-albedo feedback). A significant drawback of the models of Budyko and Sellers, nonetheless, is that the ice line remains fixed for all time. One would expect the ice line to descend to lower latitudes if the absorbed insolation were to decrease, for example, and to move poleward if the planet were to warm. Similarly, if T * (η) < T c (as with T * η1 (y) in Figure 3) ice should form, and if T * (η) > T c (as with T * ηi (y), i = 3, 4, in Figure 3) the ice line should retreat. This shortcoming was remedied in [40], in which a dynamic ice line was coupled to (2). We survey recent approaches to the analysis of these coupled equations, for both Budyko-and Sellers-type systems, in the following section.

3.
A coupled temperature-ice line model. An equation for movement of the ice line in response to changes in temperature was introduced in [40] in the simplest fashion. One considers the integro-differential system where ε > 0 is a parameter governing the relaxation time of the ice sheet (y and η are as in the preceding section). The temperature distribution T (y, t) evolves according to Budyko's equation (11a), while the dynamics of η are determined by the temperature at the ice line, relative to the critical temperature. The ice sheet retreats toward the pole if T (η, t) > T c , and moves equatorward if T (η, t) < T c .
In the analyses of (11) carried out to date expression (3) is replaced by the approximation is the second Legendre polynomial [27] (see Figure 1). We assume s(y) = 1+s 2 p 2 (y) in all that follows. We turn first to results concerning the Sellers-type version of (11) presented in [40].
a continuous approximation of (9). We note the computation of T * (y) follows precisely as in Section 2.1, with α(η) (and hence T * η (y) in (7)) again parametrized by η, given albedo function (13). The graphs of T * η (y) become continuous approximations of those plotted in Figure 3.
In this setting, E. Widiasih proved the existence of an equilibrium with a stable, small ice cap for a discrete-time version of (11), namely where T n (y) and η n are the temperature at latitude y and the position of the ice line, respectively, at year n. (Note T n = T n (y) is solely a function of y in this discrete model, with time evolution indicated by the subscript n.) The function F in (14) is given by with K the number of seconds in one year, and System (14) is in spirit an Euler's method approach to (11) with time step 1 year.
For analytical convenience, y (and so η) is embedded in the real line in [40], so that the state space becomes B × R (the reader is referred to [40] for details). The norm on B × R is given by Widiasih studied iteration of the map The intuition for her result stems from the following. Note the position of the ice line remains fixed when ε = 0 in (14). For ε = 0, fixed points of (14) thus arise when F (T, y) = 0, that is, for T * (y) = T * η (y) as in (7). Hence system (14) admits an invariant manifold of fixed points T * = (T * η (y), η) : η ∈ R when ε = 0. If the temperature profiles were to converge quickly to expression (7) relative to the movement of η (e.g., if ε was sufficiently small), then perhaps T * η (η) would serve as a good approximation to T (η) in (15). Widiasih proved this is indeed the case: Using Hadamard's graph transform method [31], she proved the existence of an attracting invariant manifold P * for map (16) for ε positive and sufficiently small, using tools such as the Contraction Mapping Theorem in certain Banach spaces. In addition, she was able to use T * η (η) to determine the dynamics on the invariant manifold. Theorem 3.1. [40] For sufficiently small ε, there exists an attracting one dimensional invariant manifold for the map m. That is, 1. There exists a Lipschitz continuous map Φ * : R → B such that There exists a closed set B ⊂ B for which the invariant manifold is attracting.
That is, given any η ∈ R, if T ∈ B then the distance between m k (T, η) and P * in the B × R-norm decreases to 0 exponentially fast, where m k denotes the k th iterate of m.
It is also worth noting that ice forms if the planet were ever ice free (h(1) < 0), and that a runaway Snowball Earth only occurs if η is ever positioned south of η 1 . More generally, Theorem 3.1 implies the behavior of solutions to the infinitedimensional dynamical system (11) is well approximated by the dynamics of the single ODE h = εh(η), with h as in (17).
dynamical systems foundation. Decreasing Q from its present value shifts the solid curve in Figure 5 downward, so that the tipping point illustrated in Figure 4 and discussed in Section 2.1 can be seen to be the result of a saddle-node bifurcation (see Figure 6). The increased analytical difficulty arising from the addition of equation (11b) must be noted, however. In the following section we highlight an alternative approach in which finite dimensional approximations to (11) are used.  Figure 6. The saddle-node bifurcation for the Sellers-type coupled model with albedo function (13). Parameters as in Figure  5.
3.2. The temperature-ice line model: Budyko-type. We begin this section by recalling Fenichel's Theorem ( [10], [11]; [17] is a nice introduction to geometric singular perturbation theory). As the case in which the invariant manifold is the graph of a function is pertinent here, we present the theorem within that framework.
Consider a system of the formẋ where x ∈ R n , y ∈ R m and ε is a real parameter. Assume f and g are each C ∞ on R n+m × I, where I is an open interval containing 0. Assume the set of rest points when ε = 0, is the graph of a C ∞ function y = h 0 (x). Assume further that Λ 0 is normally hyperbolic with regard to system (18) when ε = 0.
[10] If ε is positive and sufficiently small, there exists a manifold Λ ε within O(ε) of Λ 0 and diffeomorphic to Λ 0 . In addition, Λ ε is locally invariant under the flow of (18) and C r (including in ε) for any r < ∞.
in both x and ε, for any r < ∞.
(ii) Given that x parametrizes the manifold Λ ε , the equation describing the dynam- (iii) Changing the time scale (τ = εt), (19) can be rewritten We then haveẋ (iv) The function y = h 0 (x) (and y = h ε (x) as well) is defined on a compact subset K of R n . In our applications of Theorem 3.2 below, K can be any compact subset of R.
An analysis of temperature-ice line model (11) via finite Legendre polynomial expansions was presented in [25]. As mentioned in the introduction, the reduction of infinite-dimensional dynamical systems such as (11) to a finite number of approximating ODEs is common in a mathematical approach to low-order climate modeling. McGehee and Widiasih were the first to use such an approach to analyze an approximation of energy balance model (11) from a dynamics perspective, that is, to investigate the behavior of solutions as time varies. We conclude this section with a presentation of their approach. (We note we are now again in the setting T = T (y, t), that is, T is a function of both y and t.) In this analysis we return to the use of Budyko's discontinuous albedo function (9). Recall the model assumption of symmetry across the equator y = 0. This implies in particular that T (y, t) is an even function of y. Also recall we are using s(y) = s 0 p 0 (y) + s 2 p 2 (y), with s 0 = 1, s 2 = −0.482, and with p 0 (y) = 1 and p 2 (y) = 1 2 (3y 2 − 1) the first two even Legendre polynomials, to approximate the insolation distribution (3). Note this implies equilibrium solution (7) of Budyko's equation (2) is piecewise quadratic, with a discontinuity at η due to albedo function (9). This collectively serves to motivate the use of the approximation in the analysis of equation (2) and system (11). The problem then becomes to determine the coefficient functions u i (t) and v i (t), i = 0, 2.
To that end, substituting (20) into equation (2) yields Equating coefficients of p 0 and p 2 , respectively, in each equation in (21) (and recalling s(y) = p 0 (y) + s 2 p 2 (y)), one arrives at the four-dimensional system with T given by (22). With temperature T at the ice line η defined as in (20) we note, for future reference, T (η) = 1 2 (u 0 + v 0 ) + 1 2 (u 2 + v 2 )p 2 (η) (for ease of exposition, we suppress the dependence on t here and in the following).
The change of variables w = 1 where and in the new variables. We see equations (24b)-(24d) decouple and are each linear. As −(B +C)/R < 0, we have as t → ∞ for all initial conditions, where we have set Thus (24) admits a globally attracting invariant line on which the dynamics are governed by equation (24a) which, using (25) and simplifying, becomes One can write (27) in the form where We see that for each fixed η, equation (28) has an attracting equilibrium solution given by w * = F (η), similar in spirit to the role played by solutions (7) in the analysis of PDE (2). Indeed, plugging u * 2 and v * 2 into (25) and using the formulation of F (η) above, one can show Hence for fixed η, the variable w is simply a translate of the global mean temperature [25]. Preferred equilibrium solutions again satisfy T (η) − T c = 0. Plugging u * 2 and v * 2 into (26) allows one to write T (η) − T c = 0 equivalently as Note (30) can be written The intersection of the curves w = F (η) and w = G(η) in the (η, w)-plane then provide equilibria w * = F (η) satisfying T (η) = T c . With parameters appropriate for today's climate there are again two such solutions, illustrated in Figure 7.
To investigate the stability of these equilibria for the coupled model we reintroduce Widiasih's equation (11b) which, in the context of the approximation presented here, becomesη = ε(T (η) − T c ) = ε(w − G(η)), with G(η) as in (31). We thus arrive at the two-dimensional systeṁ When ε = 0 the ice line η remains fixed, and we have a globally attracting curve of rest points w * = F (η). By Fenichel's Theorem 3.2, system (32) has an attracting invariant manifold for ε sufficiently small, on which the dynamics are well-approximated by the single ODĖ η = ε(F (η) − G(η)).
we have returned, via ODE approximations and invariant manifold theory, to the discussion of the dynamics presented in Section 3.1 subsequent to the statement of Theorem 3.1. Once again the small ice cap equilibrium solution is stable, while the large ice cap solution is unstable, for ε sufficiently small. Alternatively, one can analyze system (32) for any ε > 0, as in [25]. Label the two equilibrium solutions (η i , w i ), i = 1, 2, as in Figure 7. Via linearization of (32) at the equilibria, McGehee and Widiasih showed (η 1 , w 1 ) is a saddle and (η 2 , w 2 ) is a sink for any ε > 0 (see Section 4.3). The stable and unstable manifolds of (η 1 , w 1 ) are also plotted in Figure 7 for a particular choice of ε. Note the saddle (η 1 , w 1 ) correlates with the large ice cap equilibrium solution for the infinitedimensional Sellers-type coupled model (Section 3.1), the latter solution admitting a one-dimensional unstable subspace in B × [0, 1]. The sink (η 2 , w 2 ) corresponds to the stable small ice cap equilibrium for the Sellers-type coupled model.
There is then a globally attracting invariant curve for system (32) for any ε > 0, comprised in part by the unstable manifold of (η 1 , w 1 ) for η ∈ (0, η 2 ). (Any extension of this part of the invariant curve to η ∈ (η 2 , 1) is not unique.) That the analysis holds for all ε > 0 is an advantage of the ODE approximation approach, relative to the function space setting of the Sellers-type coupled model. We see that the coupling of Budyko's equation (2) with an ice line equation allows for an analysis of the EBM from a dynamical systems perspective, in the settings of both infinite-dimensional function spaces and finite-dimensional ODE approximations. In the following section we show that (11) is also a highly effective EBM when modeling the extensive glacial events of the Neoproterozoic Era. In particular, and as first found using a GCM in [1], we find a stable large ice cap equilibrium solution for both the Sellers-and Budyko-type coupled models. The existence of this stable solution prevents (and serves as an alternative theory to) a runaway Snowball Earth event (recall Figure 6). The analysis of a modification of (11) appropriate to the Neoproterozoic Era requires the theory of non-smooth dynamical systems for the Budyko-type model.

Neoproterozoic glacial episodes and the Jormungand climate model.
There is ample and varied scientific evidence indicating that ice sheets flowed into the ocean in tropical latitudes during two glacial episodes in the Neoproterozoic Era, roughly 630 and 715 million years ago ( [29] and references therein). This has led many to conjecture the planet was in fact in a Snowball Earth state during this era, due to positive ice-albedo feedback. Recall the EBM presented in Section 3 also indicated the Earth would become completely ice covered if the ice line η were to advance sufficiently far south.
Evidence also exists, however, indicating that photosynthetic organisms and more biologically complex sponges flourished both before and immediately after these glacial episodes [1]. This has led others to posit that an open strip of ocean water remained ice free during these glacial maxima, implying that the runaway advance of the ice sheet toward the equator was impeded in some fashion.
A mechanism by which the ice sheet might halt shy of the equator was presented in [1]. Using a GCM parametrized for the Neoproterozoic climate, Abbot et al found that evaporation exceeded precipitation in a latitudinal belt centered around 15-20 • , due in part to the extreme cold. Ice forming at these latitudes would not be covered with snow, on average. This "bare" sea ice would then have an albedo α i , with α w < α i < α s . If the ice sheet advanced sufficiently, this decrease in albedo from α s to α i would serve to increase the absorbed insolation and raise the temperature near the ice line, potentially resulting in a stable ice cap with ice edge near the equator. Such a stable, large ice cap equilibrium solution, designated a Jormungand climate state by Abbot et al, was indeed produced using a GCM in [1]. We are interested in finding a Jormungand solution in a coupled temperature-ice line EBM of the type discussed in Section 3.
In the following sections we return to model (11), albeit with an appropriately adjusted albedo function. We first review recent work for the Sellers-type model, in which a stable Jormungand state was shown to exist [39]. We then turn to new results concerning ODE approximations to the Budyko-type coupled model which, as we'll see, leads to a non-smooth analog of h(η) in (33). As such, results from the theory of non-smooth dynamical systems will be brought to bear on the problem. 4.1. The Jormungand model: Sellers-type. Consider the infinite-dimensional dynamical model (11), with state space B × R as in Section 3.1. The primary idea introduced in [1] with regard to albedo and the Neoproterozoic Era is the following. Fix a latitude y = ρ, and assume any ice forming south of ρ has no snow cover due to the precipitation shortfall in those latitudes. A Snowball Earth event in this scenario, for example, would have bare sea ice from y = 0 to y = ρ and snow-covered ice from y = ρ to y = 1. We note the ice line η is still taken as the boundary between open water and ice, be it bare (η < ρ) or snow-covered (ρ < η).
(34) Intuitively, (34) is a continuous approximation of a function that is 3-step, piecewise constant for η ∈ (0, ρ) and equal to (13) for η ∈ (ρ, 1) (see Figure 8). Equilibrium solutions T * (y) are once again computed as in Section 2.1, using (34) in both the computation of α(η) in (5) and in equation (7). As previously, there are infinitely many equilibrium solutions T * (y) = T * η (y), one for each ηvalue. A few such solutions are plotted in Figure 9, with parameters adjusted to the Neoproterozoic Era, as discussed in [1] and [39] (in particular, T c has been set to 0 • C). Consider Widiasih's map (16), m : B × R → B × R, with albedo function (34). Precisely as before, a manifold of fixed points T * = (T * η (y), η) : η ∈ R of map m exists when ε = 0. In [39] it is shown that Theorem 3.1 applies to this Jormungand Sellers-type model, implying the existence of a one-dimensional attracting manifold P * = (Φ * η (y), η) : η ∈ R within O(ε) of T * and on which the dynamics are well approximated by the single ODĖ for ε sufficiently small. The plot of (35) is shown in Figure 10. Mimicking the analogous discussion in Section 3.1, there are three fixed points (Φ * ηi (η i ), η i ) on the m-invariant manifold, with n i close to then i illustrated in Figure 10, where n 1 < n 2 < n 3 . Note, however, the large ice cap solution corresponding to η 1 is stable in this case. For η < η 1 , h(η) > 0 and the ice line retreats toward η 1 , while for η ∈ (η 1 , η 2 ), h(η) < 0 and the ice line descends toward η 1 . Similarly, the equilibrium solution corresponding to η 3 is stable, while that with ice line at η 2 is unstable. An adjustment to the albedo function as in (34) leads to an attracting Jormungand state for the infinite-dimensional temperature-ice line coupled model (11).  Figure 10. The Sellers-type h(η) (35) for the Jormungand model. Parameters as in Figure 9.
We remark that the Snowball Earth state is unstable in the sense that h(0) > 0 in Figure 10. This is due to the increased absorption of insolation at latitudes south of y = ρ, so that bare sea ice melts (η approaches η 1 ). On the other hand, an ice free Earth remains unstable, as in Section 3.1, since we again have h(1) < 0.
Before continuing on to the ODE approach, we pause to present the bifurcation diagram shown in Figure 11. We plot the behavior of the ice line for the Jormungand model as a function of the change in insolation relative to Q Neo = 321 W/m 2 , the insolation value roughly 700 million years ago [29]. A saddle-node bifurcation occurs as Q decreases by roughly 1%. The climate then heads to the stable Jormungand state, however, rather than to a Snowball Earth, as Q decreases further. The stable Jormungand equilibrium disappears as Q decreases by about 5.5% in a second saddle-node bifurcation, at which point the systems heads to a Snowball Earth.  Figure 11. Bifurcations for the Sellers-type Jormungand model with albedo function (34) as a function of changes in the solar constant, relative to Q Neo = 321 W/m 2 . Parameters as in Figure  9.
In the following sections we examine finite ODE approximations of two new Budyko-type coupled models designed for the Neoproterozoic Era. Each model makes use of discontinuous albedo functions analogous in spirit to (34). In each case we show the equations admit an attracting Jormungand state, although the analysis leads to the theory of non-smooth dynamical systems.

4.2.
Two theorems from non-smooth dynamical systems theory. We take a moment to state two theorems useful to the analyses to follow. We start with Filippov's Theorem, which we present in the context pertinent to our work (see [12] for more general considerations).
Suppose R n can be separated into two subspaces S − and S + by a hypersurface Σ. Assume Σ is the level set of level 0 of a smooth function g, so that x ∈ Σ if and only if g(x) = 0. Let v ± : S ± → R n be smooth vector fields defined on S ± , and assume v ± each extend smoothly to Σ. For x ∈ R n , consider the differential inclusionẋ Solutions, while in S − , have flow described byẋ = v − (x), and solutions evolve according toẋ = v + (x) while in S + . For x ∈ Σ,ẋ must lie in the closed convex hull of the two vectors v − (x) and v + (x) (Σ is called the switching boundary). A solution to (36) in the sense of Filippov is an absolutely continuous function x(t) satisfyingẋ ∈ V(x) for almost all t. (Noteẋ(t) is not defined at times for which x(t) arrives at or leaves Σ.) We recall the following existence theorem [2]. (ii) If one assumes v ± are linearly bounded on S ± then, given that V is assumed to be bounded when set-valued (i.e., for x ∈ Σ), there are constants c 1 and c 2 such that V(x) ≤ c 1 x + c 2 for all x ∈ R n . In this case solutions to (36) exist on [0, ∞).
We also recall the following result for ODEs with continuous, but non-smooth, right-hand side (Theorem 3.1 in [18]).

Theorem 4.2.
Suppose v : R n → R n is continuous, and let x 0 = x(0) ∈ R n . Then there is a solution ofẋ = v(x) on the interval (−δ, δ) for some δ > 0.

Remark 3. If one assumes in addition that v(x) is linearly bounded and locally Lipschitz, then the solution guaranteed by Theorem 4.2 is both unique and defined on (−∞, ∞).
We now turn to the Budyko-type Jormungand model.
The albedo has the appropriate average value at any point of discontinuity, akin to (9). We use the term "instantaneous" since the albedo jumps from α s instantly to α i as η decreases through ρ.
Let a ± = G ± (ρ) and b ± = F ± (ρ) denote the w-values of the intersections of the nullclines shown in Figure 12 with Σ. Note a − < b + < b − < a + given our parameters. A solution of (54) starting in the region η < ρ flows according to the smooth vector field v − . If this solution hits Σ above (ρ, a + ) or below (ρ, a − ) it crosses Σ transversally. Such solutions are continuous, but not smooth, at points (ρ, w).
We see the (non-smooth) stable manifold W s (η 1 , w 1 ) of the saddle point (η 1 , w 1 ) separates basins of attraction of the two stable equilibria (ρ, w * ) and (η 2 , w 2 ). A solution crossing Σ from left to right above W s (η 1 , w 1 ) arrives at η = ρ with the "temperature" w large enough so that the ice continues to recede, with η heading toward the small ice cap solution having η = η 2 .
A solution crossing Σ just below W s (η 1 , w 1 ) as η increases also has a w-value large enough so that the ice continues to melt. In this case, however, ice starts to form once the trajectory passes through the η-nullcline. The ice line then descends and the trajectory hits Σ in finite time, with η remaining fixed at ρ while the temperature adjusts (w → w * ). In this model, (ρ, w * ) represents the stable Jormungand state.
This behavior is in contrast to that of the model with albedo function (9) (Section 3.2), in which case the trajectory starting just below the upper branch of W s (η 1 , w 1 ) heads toward the Snowball Earth scenario (η → 0) once the ice line "turns around" (recall Figure 7).
The instinctive adjustment in the albedo function from (9) to (37) to model the glacial episodes of the Neoproterozoic Era yields a stable Jormungand state when approximating (11) via ODEs. This approach, however, leads to a discontinuous vector field, with the large ice cap equilibrium solution sitting at η = ρ. An alternative to this instantaneous change in albedo as η passes through ρ could encompass a transitional phase. In the following section we argue that such a model is indeed appropriate for these periods of vast glacial extent by taking the role of Hadley cells into consideration. We present a new Budyko-type EBM that admits a stable, large ice cap equilibrium solution in the following section.
4.4. The Jormungand model: Transitional Budyko-type. The overturning circulations in the atmosphere, in which warm air rises at the equator and sinks at more northerly latitudes, are known as Hadley cells. As the warmer, moist air rises near the equator the water vapor condenses and precipitates out. The descending column of air in the northern boundary region of the Hadley cell is thus drier so that, as argued in [1], ice forming below y = ρ would not be covered in snow, on average. Abbot et al state, "... it is atmospheric dynamics that set the boundary between high-albedo snow and low-albedo bare sea ice." In related work, Poulsen and Jacobs found the Hadley circulation intensified as the sea ice advanced into low latitudes in simulations with a coupled GCM-sea ice model [30]. We use this idea in temperature-ice line EBM (11) by incorporating a transitional phase in the albedo, from α s to α i , as η decreases through ρ. Specifically, we set α w , y < η α i (η), η < y < ρ α s , ρ < y, α(y, η) = α w , y < η α s , y > η, and the temperature at the ice line (61) becomes One can rewrite the (w, v 2 )-system on Π as (note by (62) T − w is a function of η and v 2 ), and H − (η) = Ls 2 (1 − α i (η)). For fixed η, we have a globally attracting curve of rest points w * = F − (η, H − (η)) (which we note is a 4th degree polynomial). Preferred equilibrium solutions again satisfy T (η) = T c = 0 • C, which can be written w = G − (η, H(η)) with ) is a cubic polynomial in η). This leads to consideration of the function admits an attracting invariant manifold for sufficiently small ε, on which the dynamics are prescribed by a single ODE of the forṁ With the parameters aligned with the Neoproterozoic Era,η = εh − (η) has an equilibrium point atη s ≈ 0.069, so that (63) has an equilibrium solution at which η = η s ≈η s (see Figure 14). When ε = 0 the eigenvalues at rest points (η, F − (η, H − (η)), H − (η)) of (63) are λ 1 = 0, λ 2 = −B and λ 3 = −(B + C)/R. Thus for small ε, eigenvalues λ 2 and λ 3 of the Jacobian at the rest point P s = (η s , F − (η s , H − (η s )), H − (η s )) will still be negative. The remaining eigenvalue can be written in the form λ 1 = ε(h − ) (η s ) + O(ε 2 ). As (h − ) (η s ) < 0, we have λ 1 < 0, and hence P s is a sink for system (63). In this scenario P s represents the stable Jormungand state.
The approximations represented by h ± (η) are only within O(ε) of the (η, w, v 2 )system on the invariant manifolds Λ ε ± , respectively. Hence the functions h ε − (η) and h ε + (η) whose graphs provide Λ ε ± for small ε need not satisfy h ε − (ρ) = h ε + (ρ). In this case we would again have a one-dimensional Filippov system, as in Section 4.3. We note and (h ε + ) (ρ) have opposite signs. A solution with η(0) ∈ (ρ, n 1 ) would reach η = ρ in finite forward time, while a solution with η(0) ∈ (η s , ρ) would reach η = ρ in finite backward time. This Filippov equilibrium solution attracts from the right and repels from the left in this case, in contrast to that depicted in Figure 13.
If h ε − (ρ) = h ε + (ρ) then (0, η 1 ) is in the stable set of η = η s . The function is continuous, linearly bounded and locally Lipschitz on (0, 1) in this case (following from the fact h − (η) and h + (η) are quartic and cubic polynomials, respectively, and Fenichel's Theorem 3.2). By Theorem 4.2, given η 0 ∈ (0, η 1 ) there is a unique solution ofη = εh ε (η) defined for all t ∈ R and satisfying η(0) = η 0 . Given h ε (η) < 0 on (η s , η 1 ) and h ε (η) > 0 on (0, η s ), the ice line monotonically approaches η = η s from either side. We see the transitional albedo function (57) leads to a stable ice cap with ice edge close to the equator. This stands in contrast to the result arrived at when using the instantaneous albedo function (37), in which case the (Fillipov) stable equilibrium solution satisfied η s = ρ. The much larger stable ice cap in the transitional albedo setting is a consequence of the progression from α s to α i as η decreases through ρ. This leads to an overall increase in planetary albedo relative to that when using (37), thereby lowering the temperature to a significant extent.

5.
Conclusion. The models of Budyko and Sellers are remarkable in that they allow for a far-reaching investigation into the role of major climate components such as ice-albedo feedback, despite being conceptual in nature. With today's climate parameters, for example, Budyko's model exhibits an equilibrium temperature profile for which the ice line is in the vicinity of the southern edge of today's Arctic ice cap. E. Widiasih introduced a dynamic ice line that allowed her to prove the small ice cap equilibrium is attracting in an appropriately defined function space setting. She proved that temperature-ice line pairs approach equilibrium, with the ice line in the Arctic, assuming the relaxation time of the ice sheet was sufficiently large (that is, ε in (11) is sufficiently small) relative to changes in temperature.
Historically, many have analyzed the Budyko-Sellers temperature model by introducing approximations of model functions via finite Legendre polynomial expansions, thereby converting the PDE to a finite system of ODEs. McGehee and Widiasih were the first to use such an approach to analyze the coupled temperatureice line model (11). Among their findings was the existence of a stable, small ice cap equilibrium solution for the approximating system of ODEs for all ε > 0.
Surprisingly, the Budyko-Sellers model is also effective when modeling the extreme glacial events of the Neoproterozoic Era, as first introduced by Abbot, Voigt and Koll. Due to the relative paucity of precipitation in a latitudinal band centered around 15-20 • during these very cold times, ice forming in lower latitudes would generally be bare, and would therefore have a lower albedo. Via GCM simulations, the corresponding increase in absorbed insolation near the ice line was shown by Abbot et al to be sufficient to halt the advance of the ice sheet shy of the equator.
A Sellers-type albedo function incorporating bare sea ice was introduced in [39], where it was shown that Widiasih's Theorem 3.1 could be applied. The result was a proof of the existence of an attracting temperature-ice line pair with the ice line in tropical latitudes. Thus in the function space setting, the coupled model (11) was able to recreate what had first been realized via simulations with a GCM.
We have also introduced two new Budyko-type albedo functions appropriate for the Neoproterozoic glacial epochs. The instantaneous Budyko-type model is the natural generalization of the standard albedo function typically used with Budyko's equation. In this setting Filippov's theory of differential inclusions was invoked to prove the ODE approximation to the coupled equations (11) had a stable (Filippov) equilibrium solution. In this case the ice line at equilibrium rests precisely at the latitude at which the albedo switches from that of bare sea ice to that of snowcovered ice.
Again working with ODE approximations, and motivated by consideration of the behavior of Hadley cells as the glacier descends to lower latitudes, we introduced a transitional Budyko-type model. Here we relied on the theory of geometric singular perturbations to prove the existence of an attracting rest point with ice line near the equator. Thus for both Budyko-type albedo functions, we proved the existence of an attracting "Jormungand" equilibrium state for the ODE approximation to system (11). These results correlate well with the function space treatment of the Jormungand model in [39].
It is reasonable to wonder how the inclusion of land masses might affect the behavior of solutions to (11). Starting with the simplest geometric configuration, an annulus of land was incorporated into the coupled model in the function space setting in [39]. The idea was to make an appropriate adjustment to the albedo function, one that still renders the model accessible to analysis. The ODE approximations introduced by McGehee and Widiasih and advanced in this article would provide another fruitful line of inquiry into the study of the model in the case in which land was included.