Diffusive Heat Transport in Budyko's Energy Balance Climate Model with a Dynamic Ice Line

M. Budyko and W. Sellers independently introduced seminal energy balance climate models in 1969, each with a goal of investigating the role played by positive ice albedo feedback in climate dynamics. In this paper we replace the relaxation to the mean horizontal heat transport mechanism used in the models of Budyko and Sellers with diffusive heat transport. We couple the resulting surface temperature equation with an equation for movement of the edge of the ice sheet (called the ice line), recently introduced by E. Widiasih. We apply the spectral method to the temperature-ice line system and consider finite approximations. We prove there exists a stable equilibrium solution with a small ice cap, and an unstable equilibrium solution with a large ice cap, for a range of parameter values. If the diffusive transport is too efficient, however, the small ice cap disappears and an ice free Earth becomes a limiting state. In addition, we analyze a variant of the coupled diffusion equations appropriate as a model for extensive glacial episodes in the Neoproterozoic Era. Although the model equations are no longer smooth due to the existence of a switching boundary, we prove there exists a unique stable equilibrium solution with the ice line in tropical latitudes, a climate event known as a Jormungand or Waterbelt state. As the systems introduced here contain variables with differing time scales, the main tool used in the analysis is geometric singular perturbation theory.

Energy balance models (EBM) of climate are a type of conceptual model that focuses on major climate components and their interactions. Incoming solar radiation (or insolation), a planet's outgoing longwave radiation, heat transported by the atmosphere and oceans, and planetary albedo (or reflectivity) provide examples of such components. Early and important models of M. Budyko [3] and W. Sellers [33] were introduced to investigate the effect of ice albedo on temperature. An expanding ice sheet increases albedo, for example, decreasing the absorption of insolation and thereby decreasing temperature, so that the ice sheet expands further. Conversely, a retreating ice sheet lowers planetary albedo, increasing temperature and further reducing the size of the ice sheet. Might this "positive ice albedo feedback" have been the cause of snowball Earth events in our planet's past, during which it is believed glaciers advanced to the equator? Would a retreating ice sheet inevitably result in an ice free planet?
The models of Budyko and Sellers consider average annual surface temperature in latitudinal zones, so that there is one spatial dimension. The horizontal (i.e., meridional) transport of heat was modeled in [3] via the relaxation of zonal temperature to the global mean temperature (equation (9), section 2.1). Budyko's analysis indicated that a bifurcation, or climate tipping point, occurs: If the annual global average insolation decreases by roughly 5%, a small stable ice cap is lost and the climate system heads to a snowball Earth.
It should be noted the analyses of early EBM, such as Budyko's equation, focused on equilibrium solutions and their sensitivity to changes in parameters. There is no mechanism, for example, by which the edge of the ice sheet (called the ice line) is allowed to move in response to changes in temperature. In [40], E. Widiasih coupled zonal surface temperature with a dynamic ice line (with heat transport as in equation (9)). From a dynamical systems perspective, this coupled model was infinite-dimensional, comprised of a function space (temperature) crossed with an interval (the ice line). Using Hadamard's graph transform method, Widiasih proved the existence of an invariant, locally exponentially attracting, one-dimensional manifold on which the coupled system was well-approximated by a single ODE. In particular, she proved the existence of a stable equilibrium solution with a small ice cap, and an unstable equilibrium solution with a large ice cap. Widiasih's work further indicated the presence of a sufficiently large ice sheet would trigger a runaway snowball Earth event, while a very small ice cap would expand and limit on the stable small ice cap position. In particular, an ice sheet retreating to an ice free Earth state is not a possibility in this coupled model, assuming meridional heat transport is given by equation (9).
There is ample evidence indicating that ice sheets flowed into the ocean in tropical latitudes during two extensive glacial periods in the Neoproterozoic Era (roughly 630 and 715 million years ago; see [29] and references therein). This has led many to posit the existence of a snowball Earth during these times. Evidence also exists, however, indicating that photosynthetic organisms and more biologically complex species of sponge thrived before and after these glacial episodes [1]. This has led others to believe that a stable climate state with ice sheets in tropical latitudes, yet with a strip of open ocean water about the equator, existed. Abbot et al called this a Jormungand climate state in [1], proposing an adjustment in the latitude-dependent albedo function as a mechanism by which this stable state might exist (and finding a Jormungand state in an idealized, appropriately parametrized general circulation model). A variant of this albedo function was incorporated into the coupled Budyko-ice line model in [38], in which it was shown the infinite-dimensional dynamical system has a stable equilibrium solution with the ice line within 20 • of the equator, again assuming horizontal heat transport is modeled as in equation (9).
Returning to Budyko's original albedo function and continuing to use the relaxation to the mean heat transport term (9), a much simpler quadratic approximation to the coupled temperature-ice line model was introduced in [24]. This approach resulted in a reduction to a two-dimensional system of ODE that was analyzed using more elementary dynamical systems techniques. As with the original Budyko-ice line model, a small stable ice cap and a larger unstable ice cap were shown to exist. A similar quadratic approximation approach in the case of the Jormungand temperature-ice line model was followed in [37], where it was shown the quadratic approximation yielded an equilibrium solution with the ice edge in tropical latitudes. Notably, and due to the relative complexity of the Jormungand albedo function, the analysis in [37] required tools from the theory of discontinuous differential equations, as well as geometric singular perturbation theory. The latter theory will play an important role in the present work as well.
In this paper the relaxation to the mean heat transport term is replaced by diffusive heat transport in Budyko's surface temperature equation (12), leading to a variant of Budyko's model discussed in several previous works, including [17], [20], [25], [26], [28], and Chapter 10 in [14]. In each of these previous studies the ice line remains fixed in time. Budyko's zonal temperature PDE with diffusive heat transport is coupled here with Widiasih's dynamic ice line equation, and the resulting system is analyzed using the spectral method, most significantly in the case of the extensive glacial episodes of the Neoproterozoic Era. This approach produces a continuous vector field that is, however, nonsmooth on a "switching boundary" [8]. We show the vector field is locally Lipschitz, guaranteeing uniqueness of solutions. Geometric singular perturbation theory is used to show the existence of an attracting Jormungand state, with ice line in the tropics, in the case of diffusive heat transport.
For completeness, we apply the spectral method to the diffusive heat transport variant of the standard albedo Budyko-ice line coupled model as well. We also revisit work in [24] and [37] that makes use of the relaxation to the mean heat transport, albeit incorporating higher-order approximations of the temperature function to derive the governing equations. The results in this paper might then be viewed as completing the analysis of the Budyko-ice line coupled model begun in [40], [38], [24] and [37], an analysis informed by a dynamical systems perspective.
In the following section we present the Budyko-ice line coupled model. In section 3 we apply the spectral method to the diffusive heat transport model when using Budyko's standard albedo function. Section 4 contains an analysis of the Jormungand temperature-ice line model in the case of diffusive heat transport. In section 4 we prove uniqueness of solutions for the resulting nonsmooth system of ODE, as well as the existence of a stable equilibrium with ice cap extending to the tropics. The Budyko and Jormungand "relaxation to the mean" horizontal heat transport models are analyzed in section 5, albeit with higher order approximations than those used in [24] and [37]. We summarize this work in the concluding section.
2 The coupled temperature-ice line model

Budyko's conceptual model
The EBM introduced by Budyko focuses on the annual mean surface temperature in latitudinal zones. As mentioned above, analyses of EBM historically centered on equilibrium solutions and their sensitivity to changes in parameters (typically, annual globally averaged insolation). Energy balance was achieved when at a given latitude, where E in represents insolation (modulo the albedo), E out is the outgoing longwave radiation (OLR), and E tran is a meridional heat transport term.
One assumes a symmetry about the equator in Budyko's model, so that latitude extends from θ = 0 • to θ = 90 • . Additionally, one introduces the variable y = sin θ (the area of an infinitesimal band at "latitude" y has area proportional to dy), so that y = 0 at the equator and y = 1 at the North Pole. Any function dependent upon y is thus an even function of y, given the symmetry assumption.
The E in -term has the form where Q is the annual global mean insolation in Watts per meter squared (W/m 2 ), and s(y) provides for the distribution of insolation as a function of latitude (see Figure 1). Thus, Qs(y) is the average annual insolation at latitude y. The explicit formulation where β denotes the obliquity (or tilt) of the Earth's spin axis, was derived in [23]. We set β = 23.5 • , the current value of Earth's obliquity. Anticipating the use of the spectral method to follow, one can express (3) in terms of the even Legendre polynomials p 2n (y), writing where s 2n is given by One finds the insolation components fall off rapidly, to the extent the quadratic s(y) = s 0 p 0 (y) + s 2 p 2 (y), s 0 = 1, s 2 = −0.477, p 0 (y) = 1 and p 2 (y) = 1 2 (3y 2 − 1), approximates equation (3) uniformly to within 3% error (see Figure 1).
Including higher-order terms from (4) in the analysis did not alter the nature of the results in all of our investigations. Hence we replace equation (3) with equation (6) in all that follows, effectively setting s 2n = 0 for n > 1 in (4) (an exception to this approach will occur in section 5, where at times we include higher-order s 2n -terms).
As the function s(y)p 2n (y) occurs frequently below, for ease of notation we set q 2n (y) = s(y)p 2n (y) in all that follows. The function α(y, η) in equation (2) represents the albedo, which we assume depends upon the position of the ice line y = η. The albedo function used by Budyko has the form with α 1 < α 2 [3]. The larger value α 2 represents the albedo of ice, with ice assumed to exist at all latitudes poleward of the ice line. The value α 1 is the correspondingly lower albedo of the ice free surface equatorward of η. Each of s(y) and α(y, η) is dimensionless. The Earth radiates energy into space at longer wavelengths; following [3], we set where T (y, t) ( • C) is the annual mean surface temperature at latitude y, and A (W/m 2 ) and B (W/(m 2 • C) are determined empirically (see, for example, [15]). We note Sellers chose a formulation of the E out -term comprised of an empirically determined scaling of the Stephan-Boltzmann Law for blackbody radiation [33]. The E tran -term in equation (1) models the heat flux resulting from horizontal redistribution by the circulation of the ocean and the atmosphere. In each of [3] and [33], the horizontal heat transport term was given by a relaxation to the mean expression of the form where C > 0 is an empirical constant and is the global mean surface temperature. A latitude at which the temperature is below the global average will gain heat flux from neighboring latitudes, while heat will be transferred away horizontally if the local temperature exceeds the global mean.
In other variants of the Budyko-Sellers EBM ( [20], [25], [26], [28], for example), meridional heat transport is modeled by where D > 0 is an empirical constant; the form of the spherical diffusion operator (11) arises when assuming no radial or longitudinal dependence. This diffusive heat transport is akin to an eddy diffusion approach to dispersion by macroturbulence in the entire system [20]. In using either equation (9) or equation (11), note the net global effect of horizontal heat transport is zero as in each case the integral of E tran over the unit interval is zero.
The diffusive heat transport approach (11) has the advantage that the Legendre polynomials p n (y) are eigenfunctions of the spherical diffusion operator: (1 − y 2 ) d dy p n (y) = −n(n + 1)p n (y).
Hence the use of Legendre polynomials provides for a powerful technique in the analysis to follow. We consider the PDE version of Budyko's model where E tran is of the form (9) or (11). The heat capacity of the Earth's surface R (J/(m 2 • C)) is assumed to equal the global average value; the units on each side of equation (12) are W/m 2 .

Widiasih's dynamic ice line
Note the lack of any mechanism in Budyko's model by which the ice sheet (represented in (12) by the ice line η) is allowed to expand or retreat as the local surface temperature varies. In the case E tran = C(T − T ), a dynamic equation for the ice line was introduced in [40] and coupled to equation (12), resulting in the integro-differential system The parameter ε > 0 governs the relaxation time of the ice sheet. T c is a critical temperature, meaning that if T (η, t) < T c ice forms and the ice line moves toward the equator, whereas the ice line retreats toward the pole if T (η, t) > T c . System (13) has been analyzed in [40] and [38] when assuming the albedo varies smoothly across the ice line (in this case the model is said to be of "Sellers-type"). Using a smooth variant of albedo function (7) and assuming ε satisfies certain bounds, Widiasih proved that a stable equilibrium temperature profile-ice line pair (T * (y), η * ) exists with η * at roughly 70 • N (a small ice cap). The analysis in [40] centered on the use of Hadamard graph transforms.
With a smooth variant of albedo function (27), chosen to model the very cold world of the great glacial epochs of the Neoproterozoic Era (see section 4), Walsh and Widiasih proved the existence of a stable equilibrium for system (13) with the ice line in tropical latitudes in [38] (with constraints placed on ε). The choice of albedo function was motivated by [1], in which such a stable equilibrium was found using an idealized, appropriately parametrized general circulation model, and subsequently dubbed a "Jormungand" state. That the temperature-ice line coupled model (13) can be adjusted to find a stable climate appropriate for current times, as well as a stable climate appropriate for the deepest of ice ages, points to the power and utility of this relatively simple coupled EBM.
The present work considers discontinuous albedo functions-in which case equation (12) is said to be of "Budyko-type"-and diffusive heat transport. In the following two sections we present the analysis of the coupled temperature-ice line model assuming diffusive heat transport (11) and albedo functions (7) and (27), respectively. For completeness, in section 5 we then revisit and expand upon the quadratic approximation approaches to analyzing system (13) presented in [24] and [37].

Diffusive heat transport: standard albedo
In this section we consider the system with albedo function α(y, η) given by (7). We note equation (14a) comes with the boundary conditions that the gradient of the temperature profile T (y, t) equals zero at the equator and at the North Pole.

The spectral method
Recall functions appearing in system (14) are even functions of y, due to the assumed symmetry about the equator. We approximate system (14) by expanding each function of y in terms of the first N + 1 even Legendre polynomials, arriving at a system of N + 2 ODE as follows. Consider the expression and first note dp 2n dθ = dp 2n dy cos θ = 0 at the equator and North Pole. Hence expression (15) satisfies the prescribed boundary conditions. Substitute (15) and equation (6) into equation (14a) to arrive at where the dot represents differentiation with respect to time, s 0 and s 2 are given in (6) (recall we set s 2n = 0 for n > 1), and Recalling q 2n (y) = s(y)p 2n (y) is a polynomial of degree 2n + 2, we note for future reference α 2n is a polynomial in η of degree 2n + 3. Also note that substituting y = η into expression (15) gives the temperature at the ice line Equating coefficients of p 2n in equation (16) and simplifying yields an approximation of the infinite-dimensional system (14) given by the N + 2 ODE Note system (19) can be rewritteṅ f 2n (η) = 1 (B + 2n(2n + 1)D) (Q(s 2n − α 2n (η)), n = 1, ..., N.
We analyze system (20) via geometric singular perturbation theory.

Fenichel's Theorem
We assume the ice line moves slowly relative to the evolution of the temperature. Thus η will be a slow variable, while T 0 , ..., T 2N are fast variables. If ε = 0 then η remains fixed, in which case T 2n → f 2n (η) exponentially fast for n = 0, 1, ..., N . Hence there is an attracting manifold of rest points when ε = 0. We thus analyze system (20) via geometric singular perturbation theory, for which we now recall Fenichel's Theorem ( [11], [12], [19]). The case in which the attracting invariant manifold when ε = 0 is the graph of a function is sufficient for our purposes, so this is the setting we present here.
Consider a system of the formẋ where x ∈ R n , u ∈ R m and ε is a real parameter. Assume f and g are each C ∞ on I is an open interval containing 0. Assume the set of rest points when ε = 0, is the graph of a C ∞ function x = h 0 (u). Assume further that Λ 0 is normally hyperbolic with regard to system (22) when ε = 0.
Theorem 3.1 [11] 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 system (22) and C r (including in ε) for any r < ∞.
Remark 1. (i) In this setting, Λ ε is the graph of a function x = h ε (u) that is C r , in both u and ε, for any r < ∞.
(iv) The function x = h 0 (u) (and x = h ε (u) as well) is defined on a compact subset K of R n . In the applications of Theorem 3.1 to follow, K can be any compact subset of R.

Standard albedo model dynamics
Returning to system (20), we see that when ε = 0 there is a globally attracting curve of with f 2n as in (21), n = 0, 1, ..., N . By Fenichel's Theorem, system (20) has an attracting invariant manifold for sufficiently small ε, on which the dynamics are well approximated by the single ODEη Recalling f 2n is a polynomial of degree 2n+3, we see h(η) is a polynomial of degree 4N +3.
We use the following parameter values [24], chosen to align with the modern climate: With the above parameters set, the existence of a small stable ice cap then depends upon the value of the diffusion constant D. A larger D-value corresponds to more efficient heat transport from the warmer tropical region to northern latitudes, eliminating the possibility of any small ice cap (see Figure 2(a)). Smaller D-values induce a steeper equator-to-pole temperature gradient, and the possibility of a small stable ice cap exists (Figure 2  Proof. This result follows from Fenichel's Theorem and the analysis of the degree 7 polynomial (25) h Figure 2(b)). One can show h(η) has two real roots η 1 < η 2 in [0, 1], with h (η 1 ) > 0 and h (η 2 ) < 0. System (19) then has a stable equilibrium solution well approximated by (h 0 (η 2 ), η 2 ) and an unstable equilibrium solution near (h 0 (η 1 ), η 1 ), with h 0 as in equation (24). Figure 2(c) are higher-order approximations (N = 2, 5 in equation (15)). As noted in [28], each succeeding term in (15) contains system information pertaining to ever smaller spatial scales. Though not apparent in Figures  2(b) and 2(c), including higher-order terms serves to slightly reduce the size of the ice cap.

Remark 2. (i) Also plotted in
We also see the quadratic approximation with N = 1 suffices to capture the dynamics of system (19).
(ii) It is interesting to note the lack of any unstable equilibrium with ice line lying poleward of the stable ice line position when D = 0.35 (Figure 2(b)). Previous analyses of (uncoupled) equation (14a) did find such a stable-unstable pair of equilibria, with the existence of the unstable ice cap known as small ice cap instability (SICI) (see [27] for an excellent discussion of SICI and Budyko's model with diffusive heat transport). The value D = 0.35 is evidently sufficiently small in the coupled model (19) to preclude the existence of a small unstable ice cap. We note the SICI phenomenon does arise when D is increased to 0.394 (see Figure 3; compare with Figure 3 in [27]). Thus with more efficient meridional heat transport, the coexistence of a small stable ice cap and a "stable" ice free Earth becomes a possibility, in contrast to the results for the relaxation to the mean heat transport model ( [40], [24]; also see Figure 6).
(iii) Integrating equation (15) with respect to y over the unit interval yields as 1 0 p 2n (y)dy = 0 for n > 1. Hence T 0 represents the global mean surface temperature, a quantity that varies with the ice line η in the coupled model (19). Given the evolution of η depends on the D-value in system (19), the global mean temperature depends upon D as well. At an equilibrium point with η = η * , the global mean temperature is  (recall equations (20) and (21)). With N = 1 and D = 0.35, the stable small ice cap is positioned near η * = 0.837 (57 • N), for which T * 0 = 10.9 • C. In our present climate, η * ≈ 0.940 (70 • N) and T * 0 ≈ 15 • C, so that the model ice cap is too large and the model world is too cold when D = 0.35. Were one to require the ice line for the small stable ice cap to lie near η * = 0.940, one could compute the corresponding D-value by solving h(0.94) = f 0 (0.94) + f 2 (0.94)p 2 (0.94) − T c = 0 for D (recall f 2 depends upon D). In this case one finds D = 0.394, for which T * 0 = 14.6 • C, a good approximation to our current climate.
We now turn to diffusive heat transport and the spectral method when system (14) is adjusted to model the great glacial episodes of the Neoproterozoic Era. We incorporate an albedo function, modeled after a similar function introduced in [1] and appropriate for near snowball Earth events, into system (14). If the ice sheet were to extend down to tropical latitudes, what might serve to counteract positive ice albedo feedback to the extent a runaway snowball Earth event never occurs, with the ice line stabilizing in tropical latitudes?
A possible mechanism for this dynamic was presented in [1]. Using an idealized general circulation model (GCM) parametrized for the extensive ice ages of the Neoproterozoic Era, Abbot et al found that evaporation exceeded precipitation in a latitudinal band centered around 15-20 • . Any ice forming in these lower latitudes would have no (or significantly less) snow cover, and hence this "bare" ice would have an albedo value α i less than albedo value α 2 , but with α i greater than the albedo value α 1 of the surface south of the ice line. This reduced albedo, relative to the albedo of snow covered ice, was sufficient to enable Abbot et al to find a stable equilibrium solution with the ice line in tropical latitudes when using the GCM in [1]. This stable "Jormungand" state, with a narrow strip of ocean about the equator remaining ice free, provides a possible explanation for the survival of photosynthetic organisms through these tremendous ice ages.

The nonsmooth model equations
The use of the function (27) in system (14) leads to a nonsmooth vector field as follows.
While V − is a smooth vector field on R N +2 , in the model we restrict V − to the set Similarly, let V + be the vector field given by the model equations when ρ < η, that is, V + corresponds to the system of N + 2 ODĖ While V + is a smooth vector field on R N +2 , we restrict V + to the set Vector fields V − and V + differ only in the expressions a 2n (η) (30) and b 2n (η) (29). Note that V − is defined at points (x, ρ), and further that a 2n (ρ) = α 2 s 2n − (4n + 1)(α 2 − α 1 ) ρ 0 q 2n (y)dy) = b 2n (ρ), that is, V − and V + agree on the set The set Σ is known as a switching boundary (or discontinuity boundary) [8] for the vector field The vector field V is smooth on R N +2 \Σ and continuous on R N +2 . Standard ODE theory then implies solutions to the initial value problem exist. However, V fails to be C 1 on Σ, as we now show. Let J − and J + denote the Jacobian matrices of vector fields V − and V + , respectively.
While V is not C 1 at points in Σ, we do have that V is locally Lipschitz, thereby guaranteeing that the model ODE (37) has unique solutions.
Proposition 3. The vector field V (x, η) given by system (36) is locally Lipschitz, that is, for any z = (x, η) ∈ R N +2 , there exists a constant L(z) and δ > 0 such that for all Proof. Each of V − and V + is smooth on R N +2 , implying that V − and V + are each locally Lipschitz on R N +2 . It follows that V is locally Lipschitz on R N +2 \Σ.
Recall that for all (x, ρ) ∈ Σ, V + (x, ρ) − V − (x, ρ) = 0. By Proposition 2, J − (x, ρ) = J + (x, ρ). In this scenario the switching boundary Σ is said to be uniform with degree 2 (section 2.2 in [8]). In particular, given z 0 = (x 0 , ρ) ∈ Σ, there exists a δ > 0 and a smooth function F : B δ (z 0 ) → R N +2 satisfying, for all z = (x, η) ∈ B δ (z 0 ), This follows from consideration of the Taylor series of the function V + − V − expanded about z 0 = (x 0 , ρ). (Note if we disregard connections to the model, each of V − and V + is defined on R N +2 , and so each is defined on both sides of Σ.) The quantity η − ρ may be chosen as one of the coordinates in regions close to Σ and, since V + − V − vanishes on Σ, it follows that each term in the Taylor series expansion of each component function of V + − V − has a factor of η − ρ. The expression η − ρ is raised to the first power in equation (38) since the Jacobian matrix of V + −V − does not vanish on Σ by Proposition 2.
As each of V − and V + is locally Lipschitz at z 0 , pick positive numbers δ 1 , δ 2 , L 1 and L 2 such that for all z 1 , z 2 ∈ B δ 1 (z 0 ), and for all z 1 , z 2 ∈ B δ 2 (z 0 ), Without loss of generality, assume z 1 = (x 1 , η 1 ) ∈ S − and z 2 = (x 2 , η 2 ) ∈ S + . Then and we have that V is locally Lipschitz on Σ. Combined with the fact V is locally Lipschitz on R N +2 \Σ, we have the desired result.

Jormungand diffusion model dynamics
We see via Proposition 3 that solutions to the nonsmooth system (36) exist and are unique. However, we cannot directly apply geometric singular perturbation theory to system (36) due to Proposition 2: Fenichel's Theorem applies to smooth vector fields (at least C 1 ), and system (36) is not smooth on the switching boundary Σ. As we are interested in a stable equilibrium solution with ice edge to the south of y = ρ, we first apply the theory to system (31), subsequently restricting the domain of η to return to the model interpretation. We repeat this process for η > ρ (vector field (33)), and then use Proposition 3 to analyze the behavior of system (36). As in section 3, system (31) has a globally attracting curve of rest points when ε = 0. By Theorem 3.1, for ε sufficiently small system (31) has an attracting invariant manifold on which the dynamics are well-approximated by the ODĖ As each polynomial f − 2n has degree 2n + 3, we have that h − (η) is a polynomial of degree 4N + 3.
It is known that 600-700 million years ago the solar constant Q was roughly 94% of its current value of 343 W/m 2 . Many other parameter values are taken from [1], where the main criterion for the existence of a stable Jormungand state was found to be a significant difference between snow-covered ice albedo α 2 and bare ice albedo α i . Relative to that used in section 3, we have chosen the smaller value D = 0.25 as Abbot et al argue that meridional heat transport was less efficient during these periods of extensive glacial cover. Finally, the drawdown of atmospheric CO 2 via silicate weathering would have been greatly reduced due to the lack of significant ice free continental surfaces. As atmospheric CO 2 absorbs outgoing longwave radiation, we have reduced the A+BT -term in equation (14) by decreasing A. (40) and N = 1, system (36) has a stable Jormungand equilibrium solution for ε sufficiently small. That is, there is a stable equilibrium solution that does not correspond to a snowball Earth event and for which the ice line lies within 20 • of the equator.
One can apply Theorem 3.1 to system (33) as well, and then restrict to the set S + . For ε sufficiently small, system (33) is well-approximated by the ODEη = εh + (η), where when ρ < η. We plot h + (η) in albedo function (27) and Hadamard graph transform techniques (see Figure 8 in [38]). Note the graph of the function is continuous at η = ρ. As with the agreement of V − and V + on Σ, this follows from the fact a 2n (η) and b 2n (η) agree when η = ρ. The graphs of the invariant manifolds for vector fields (31) and (33), however, are only within O(ε) of the graphs of h − and h + plotted in Figure 4(a). Nonetheless, system (36) does admit unique solutions by Proposition 3. Hence for ε sufficiently small, a trajectory of system (36) with η(0) above η = ρ, but below the unstable equilibrium ice line η = η 2 , will see η decrease as the trajectory passes through the switching boundary Σ. As seen in Figure 4(a), h − (ρ) < 0, and the η-value along this trajectory will continue to decrease, limiting on the large, stable ice line position near η = η 1 . Interestingly, only the large stable ice cap remains when higher-order approximations are used (N = 2, 5 in equation (15); see Figure 4(b)). The incorporation of smaller spatial scale information into the model makes for a colder world; in this case every trajectory converges to the stable Jormungand state.

Bifurcations
We note that several bifurcations occur in the Jormungand model with diffusive heat transport. With several parameters to choose from, the bifurcation parameter used here is the constant A (recall the A+BT -term in equation (12) models the outgoing longwave radiation). Atmospheric greenhouse gases such as CO 2 absorb OLR, so a decrease in A might then be viewed as an increase in atmospheric CO 2 , while an increase in A can be associated with a decrease in CO 2 . In this way the parameter A can be viewed as a proxy for atmospheric greenhouse gas concentrations, although the rise and fall of the former corresponds to the fall and rise of the latter.
In Figure 5 we plot the position of the ice line at equilibrium as the parameter A is varied (with N = 5 in (15)). Not surprisingly, if A is large enough (CO 2 concentrations are sufficiently low) the climate system heads to a snowball Earth as η → 0. If A is sufficiently small, so that OLR is absorbed to a great extent, the model tends toward an ice free Earth with η → 1.
A first bifurcation occurs as A increases through roughly 150, with the creation of an unstable equilibrium solution with a small ice cap, and a stable equilibrium solution with a slightly larger ice cap. This unstable solution is a manifestation of the SICI phenomenon mentioned above. For A-values less than roughly 157 there is no stable Jormungand state as every solution with η(0) near 0 converges to the small stable ice line position.
For A between roughly 157 and 161.5 the system exhibits bistability, with a stable equilibrium solution with a moderately sized ice cap and a stable Jormungand state. As A increases through 161.5 there is a saddle node bifurcation, leaving the stable Jormungand state as the sole equilibrium solution.
As A continues to increase, so that the planet is radiating to space ever more efficiently, the stable ice line position heads to tropical latitudes. A final bifurcation occurs at A ≈ 187, where the Jormungand state disappears in another saddle node bifurcation, with the system now experiencing a runaway snowball Earth event.
It is interesting to note the bifurcation diagram plotted in Figure 5 exhibits many of the characteristics of the bifurcation diagram for the full infinite-dimensional system (13) when using a smooth variant of the Jormungand albedo function (27) (see Figure 9 in [38]). One exception is the appearance of SICI in the Jormungand diffusion model; by contrast, in [38] the ice free Earth state is unstable in the sense there is no equilibrium between the small stable ice line position and the North Pole.

Relaxation to the mean heat transport
In this section we revisit and expand upon approximations of system (13) presented in [24] and [37], investigating the effect (or lack thereof) the inclusion of higher-order terms in equation (15) has on results previously obtained. We will also allow for the use of higher-order s 2n -terms in equation (4).

Budyko's albedo function
Consider system (13) with albedo function (7); recall T denotes the global mean surface temperature. It can be shown that equilibrium temperature profiles for uncoupled equation (13a) have the form where α = 1 0 α(y)s(y)dy = α 1 η 0 s(y)dy + α 2 1 η s(y)dy (see [24], e.g.). Assuming s(y) as in (4) is truncated at n = N , T * (y, η) is a piecewise degree 2N polynomial in y with discontinuity at y = η, for each fixed η. Hence one might consider the use of Legendre polynomial approximations, given the form of the equilibrium solutions of equation (13a), even though their role as eigenfunctions is no longer relevant. This idea was originally suggested in a preprint version of [24] in the case N = 1, that is, when using quadratic approximations. Here we carry out the analysis for arbitrary N , while omitting the more tedious computational details.
Hence the long term behavior of system (48) reduces to that of equation (48a) when assuming variables other than u and η are at equilibrium. Coupled with equation (13b), the resulting two equations can be placed in the forṁ When ε = 0 system (51) has a globally attracting curve of rest points (u, η) = (f (η), η).
for ε sufficiently small. In general, h is a polynomial in η of degree 2N + 1.
In the case N = 1 the model development and analysis above reduces to that presented in section 3.2 in [37]. Using quadratic approximations for s(y) and piecewise for T (y, t) as well, the function h(η) given by equation (52) is a cubic polynomial with two zeros in [0,1] when using parameter values (26) (and C = 3.09). One is a stable equilibrium with η near η 2 ≈ 0.94, while the other is an unstable equilibrium with η near η 1 ≈ 0.24 (see Figure 6(a)).
If we restrict to the use of quadratic polynomials in the approximations, the behaviors of the diffusion and relaxation to the mean models are qualitatively the same for appropriately chosen D-values (compare Figures 6(a) and 2(b)). This is not surprising at first glance, for if η is kept fixed and equation (15) with N = 1 is substituted into equation (11), the result is −6DT 2 (t)p 2 (y). Note in the N = 1 case T (y, t) = T 0 (t) + T 2 (t)p 2 (y), We see that for a second mode approximation, and with the ice line fixed, the heat transport in the diffusion model is treated the same as that in relaxation to the mean model, varying only by the choice of constants D and C (see [26]). We further note that if we set s 2n = 0 for n > 1 as in previous sections, T 2n → 0 and V 2n → 0 as n → ∞ for n > 1, implying a return to the model behavior when using a quadratic approximation for U (y, t) and V (y, t) (and so for T (y, t)).
The only remaining question pertains to the use of higher order terms in expression (4) for s(y). The s 2n -terms appear only in the equilibrium point expressions (50). That is, T 2n and V 2n still converge to equilibrium points but, if s 2n = 0 for n > 1, T * 2n and V * 2n are no longer 0. This in turn has an effect on h(η) in equation (52). This effect is negligible, however, as can be seen in Figure 6(b). The qualitative analysis of the corresponding ODE (52) remains unchanged from the case in which s 2n = 0 for n > 1. The analysis presented in [37] using quadratic polynomials, following that presented in [24], is sufficient to understand any approximation to system (13) via the use of Legendre polynomials.

The Jormungand albedo function
In this section we consider approximations to system (13) in the case where the albedo function is given by (27). We seek to establish the existence of a stable equilibrium solution with a large ice cap, as was the case with the diffusion model in section 4. Note when η > ρ, albedo functions (7) and (27) are identical, and hence the analysis follows that presented in the previous section. However, we now let h + (η) denote the function previously defined as h(η) in equation (52), the plus sign indicating the restriction of η > ρ in the Jormungand setting. The plot of h + (η) can be seen in Figure 7.
We are lead to consider the one-dimensional differential inclusioṅ to which the theory of Filippov [13] can be applied. The analysis of this Filippov flow in the case N = 1 was presented in section 4.3 in [37], to which we refer the reader for details. In Figures 7(a) and 7(b), parameters were chosen as in Figure 13 in [37]. Referring to Figures 7(a) and 7(b), system (65) has a stable equilibrium solution with ice line near η = η 3 , and an unstable equilibrium solution with ice line near η = η 2 , for sufficiently small ε. We note a similar (η 2 , η 3 ) unstable-stable pair arose in the N = 1 case in the Jormungand diffusion model (recall Figure 4(a)).
Referring to the plot of h − (η) and h + (η) in Figure 7(a), a solution with η(0) ∈ (0, ρ) will see η increase and reach η = ρ in finite time. A solution with η(0) just to the right of η = ρ will see η decrease and again reach η = ρ in finite time. In this case the stable Jormungand state has the ice line sitting at η = ρ.
In Figure 7(c) we use the parameter values given in (40). In this case every solution has η(t) reaching η = ρ in finite time. This behavior is the discontinuous analog of the behavior of trajectories in the diffusion model when higher-order terms in (15) were included in the analysis (recall Figure 4(b)).
If we set s 2n = 0 for n > 1, then each of T 2n , V 2n and W 2n approaches 0 over time for n = 1, ..., N . Hence the long term behavior of system (54) in this case is equivalent to the N = 1 case (quadratic approximations). As with Budyko's albedo function in section 5.1, including higher-order terms from expression (4) in the analysis does not alter the model behavior in any qualitative sense (see Figure 7(b)). The use of quadratic polynomials is sufficient when approximating the temperature-ice line system (13) with   Figure 6). Parameters for (a) and (b) as in Figure 13 in [37]. (c) N = 1, parameters as in (40).
Legendre polynomials when using the Jormungand albedo function.

Conclusion
The energy balance models introduced by Budyko [3] and Sellers [33] were used to investigate the role positive ice albedo feedback plays in climate dynamics. Budyko used a relaxation to the mean horizontal heat transport model component and the straightforward albedo function (7). E. Widiasih enhanced the model by coupling Budyko's temperature equation to a dynamic ice line [40]. Using a smooth variant of albedo function (7) and working in the function space setting, Widiasih proved the existence of a stable equilibrium with a small ice cap, and an unstable equilibrium with a large ice cap, provided the time constant ε for the evolution of the ice line met certain conditions. In particular, she found that positive ice albedo feedback did not lead to an ice free Earth, although this feedback did lead to a snowball Earth if the ice line were ever positioned sufficiently far south. By adjusting the smooth albedo function, Walsh and Widiasih proved that a stable equilibrium solution with the ice line in tropical latitudes could be realized in the coupled model in the function space setting [38], again assuming certain restrictions on ε. This stable "Jormungand" climate state had previously been realized using an idealized GCM in [1], a work in which the aim was to model certain extensive glacial episodes from Earth's past.
A second approach to studying the temperature-ice line model was presented in [24], in the case where the albedo was given by the piecewise constant function (7) (and still incorporating relaxation to the mean meridional heat transport). McGehee and Widiasih used quadratic polynomial approximations for the temperature profile and deduced that the corresponding system dynamics reduced to that of a pair of ODE. Once again an unstable-stable pair of equilibria were found, similar in spirit to the result of Widiasih in the function space setting.
A similar quadratic approximation approach was used in [37], albeit with piecewise constant albedo functions (including (27)) modeled after an albedo function used in [1]. Although tools from nonsmooth dynamical systems theory had to be appealed to in the analysis, the results were akin to those found in [38], namely, that the system arrived at via quadratic approximations exhibited a stable Jormungand climate state for ε sufficiently small.
In this paper we have replaced the relaxation to the mean heat transport with a diffusive process. Many authors have previously analyzed Budyko's equation with diffusive heat transport, in each case, however, with the ice line fixed. Due to the form of the spherical diffusion operator, approximations with Legendre polynomials is a natural tool to use in the analysis. We thus coupled the temperature equation with the ice line equation and applied the spectral method, with diffusion constant D. As we assume the ice line moves slowly relative to changes in the temperature, the main instrument used to analyze the approximating systems of ODE is geometric singular perturbation theory.
When using Budyko's albedo function (7) we find the diffusive model, in essence, recreates the dynamics seen in [40] and [24] for a range of D-values. If the diffusive transport is too efficient (D is too high), however, the small stable ice cap is lost and a retreating ice line tends to the North Pole. For other ranges of D-values, a small stable equilibrium ice cap and an even smaller unstable equilibrium ice cap coexist. This small ice cap instability (SICI), previously known to occur in Budyko's (uncoupled) model, can therefore arise in the coupled model as well. Incorporating modes of order greater than 2 in the analysis did not change the qualitative nature of model behavior.
We have applied the spectral method to the coupled model with diffusive transport using the piecewise constant Jormungand albedo function (27). While the resulting vector field is nonsmooth due to the existence of a switching boundary, the vector field has been shown to be locally Lipschitz. This spectral approach produces a stable Jormungand state for a range of D-values, akin to the results in [38] and [37]. Incorporating higher modes into the approximation, however, changes the dynamics in the Jormungand model. For example, the model can exhibit bistability when using a 2-mode approximation, while the stable Jormungand state becomes the sole equilibrium solution when using a 4-mode approximation.
Each of the results for the coupled model with diffusive heat transport were arrived at via geometric singular perturbation theory and, as such, ε is assumed to be sufficiently small.
For completeness we revisited both the Budyko and Jormungand coupled models in the case of relaxation to the mean heat transport, albeit using higher mode approximations than those used in [24] and [37]. In every case we found approximations by quadratic polynomials was sufficient to discern the behavior of solutions to the corresponding systems of ODE.
We view this work as completing the analysis of the dynamical system that, originating with Budyko's seminal surface temperature equation [3] and its coupling with Widiasih's dynamic ice line equation [40], was further analyzed in [38], [24] and [37].
Acknowledgment. The author recognizes and appreciates the support of the Mathematics and Climate Research Network (http://www.mathclimate.org).
Appendix. We derive model equations (57). Substituting U, V and W from (53) into (13a) yields the three equations Equating coefficients of p 2n (y) in each of the above equations then gives system (54). The associated expression T in (55) is a function of T 2n , V 2n , W 2n and η, while T (η, t) is given by (56).