Neoproterozoic glacial ages in a nonsmooth conceptual model

There is widespread agreement that ice sheets flowed into the ocean in tropical latitudes at sea level during the Earth’s past. Whether these extreme ice ages were snowball Earth events, with the entire surface covered in ice, or whether ocean water remained ice free in regions about the equator, continues to be controversial. For the latter situation to occur, the effect of positive ice albedo feedback would have to be damped to stabilize an advancing ice sheet shy of the equator. In this paper we analyze a conceptual model comprised of a zonally averaged surface temperature equation coupled to a dynamic ice line equation. This difference equation model is aligned with the cold world of these great glacial episodes through an appropriately chosen albedo function. Using the spectral method, the analysis leads to a nonsmooth singular perturbation problem. The Hadamard graph transform method is applied to prove the persistence of an invariant manifold, thereby providing insight into model behavior. A stable climate state with the ice line resting in tropical latitudes, but with open water about the equator, is shown to exist. Also presented are local smooth and nonsmooth bifurcations as parameters related to atmospheric CO2 concentrations and the efficiency of meridional heat transport, respectively, are varied.


Introduction
The theory of nonsmooth dynamical systems is playing an ever-expanding role in the development and analysis of conceptual climate models. Nonsmooth models arising in the study of the overturning ocean circulation include [1][2][3] and [4]. Several works have appeared concerning nonsmooth models of sea ice loss in the Arctic ( [5][6][7], for example). Moreover, nonsmooth models often appear when modeling the glacial cycles, where the lack of smoothness may be associated with the dynamic between ice sheets and atmospheric carbon dioxide [8] or the release of carbon dioxide from the deep ocean [9]; with changes in Antarctic ice volume [10]; or with changes in glacial ablation rates [11].
The present work focuses on surface albedo, which plays a significant role in the evolution of a planet's climate. The term albedo refers to the extent to which the surface reflects incoming solar radiation (or insolation). Snow and ice, for example, have higher albedos than does ice-free surface. Notably, ice albedo interacts with climate in a positive feedback loop. Were an existing ice sheet to expand, the surface albedo would increase, thereby lowering temperature and triggering further ice sheet growth. On the other hand, were an ice sheet to retreat surface albedo would decrease, leading to greater absorption of insolation and an increase in temperature. This would lead to further reduction in the size of the ice cover. (See [12] for the effect the recent and ongoing loss of ice in the Arctic has had on Arctic temperatures, for example.) There is widespread agreement that during two ice ages in the Neoproterozoic Era, roughly 715 and 630 million years ago (mya), respectively, ice sheets flowed into the ocean at sea level in tropical latitudes [13]. Many believe these extensive glaciations were snowball Earth episodes, with ice covering the entire planet, triggered in part by positive ice albedo feedback (see [14] and references therein). There are others, however, who argue the advance of the ice sheet likely stopped shy of the equator, leaving a strip of open ocean water (see [15] and references therein). The actual extent of the ice cover during these ice ages remains controversial [13].
Conceptual climate models, focusing on notions of energy balance, were introduced by M Budyko [16] and W Sellers [17] in 1969 to investigate positive ice albedo feedback. Having latitude as the sole spatial variable, each model focused on the effect ice albedo has on zonally averaged surface temperature. Sellers' model indicated a reduction of the solar constant by 2%-5% was sufficient to trigger an ice age. Budyko's work suggested that once the ice sheet extended down below a critical latitude, a snowball Earth event inexorably followed. The solar constant 600-700 mya was 6% less than today's value, so that insolation absorbed at the surface was significantly lower during these times relative to our present climate. The modeling of Sellers and Budyko, together with physical evidence [14], support the proposition that the Neoproterozoic ice ages mentioned above were snowball Earth events.
There is also evidence, however, supporting the existence of open water near the equator during these ice ages [15]. For example, photosynthetic organisms and other life forms are known to have survived these glaciations. A mechanism by which positive ice albedo feedback might be sufficiently damped so as to ensure the advance of a large ice sheet stopped prior to reaching the equator was introduced in [15]. Abbot et al argued that a stable climate state, with glaciers extending to the tropics but with open water encircling the equator, could exist given certain assumptions concerning the albedo, to be described below. With seasonal changes causing the open ocean to snake around the equator, this stable climate state was called a Jormungand state in [15], after the sea serpent from Norse mythology.
We incorporate the assumptions from [15] into an energy balance surface temperature model coupled with a dynamic ice line, assuming diffusive meridional heat transport. The discrete-time dynamical systems approach to the analysis of the model assumed in this work leads to consideration of a singular perturbation problem for which the unperturbed invariant manifold is nonsmooth. Hence the existing theory for normally hyperbolic invariant manifolds ( [18] and references therein) cannot be directly applied. The main result presented herein is the proof of the persistence of an invariant manifold under singular perturbation and its consequences for understanding model behavior and bifurcation scenarios. The technique used in the proof is the Hadamard graph transform method [18,19].
The components of the model, including an albedo function designed for the Neoproterozoic ice ages under consideration and following that presented in [15], are introduced in section 2. In section 3 the piecewisedefined model difference equations are described. While the system of equations is everywhere continuous, the system is shown to be nonsmooth on a certain hyperplane in phase space.
The model is placed within the context of singular perturbation problems, with singular parameter  representing the time constant for movement of the ice sheet, in section 4. Additionally, the main result and subsequent analysis of model behavior, including the long-term dynamics and bifurcation scenarios, are presented in section 4. We provide concluding remarks in section 5, while the proof of our main result, and an explication of its role in determining the temperature-ice line dynamics, can be found in the appendix.

The temperature-ice line model
Let T n (y) (°C) denote the average annual surface temperature at y=sin θ, where θ is the latitude, at year n. The variable y, referred to as 'latitude' in all that follows, is convenient in that a latitudinal band at y has area proportional to dy. The energy balance equation for temperature has three terms, corresponding to absorbed insolation, outgoing longwave radiation, and meridional heat transport, respectively. A diffusion approach for the heat transport term is used here, as is the case, for example, in [20,[21][22][23] and [24]. Meridional heat transport refers to the heat flux across latitude circles associated with atmospheric and oceanic currents.
Consider the difference equation in which β denotes the tilt (or obliquity) of the Earth's spin axis, serves to incorporate insolation as a function of latitude [25]. For the planet Earth, s(y) monotonically decreases from a maximum s(0) at the equator to a minimum s(1) at the North Pole. The model assumes ice exists everywhere above the edge of the ice sheet, denoted η and called the ice line, with no ice at latitudes equatorward of η. The function a h ( ) y, represents the surface albedo, which depends upon the position of the ice line. Hence , models the incoming solar radiation absorbed at the surface at latitude y.
The (A+BT n )-term represents the outgoing longwave radiation, with A (Wm −2 ) and B (W(°C m 2 ) −1 ) estimated by satellite data [26]. The diffusive meridional heat transport term is assumed to be proportional to the surface temperature gradient, and is the form the spherical diffusion operator takes when assuming no radial or longitudinal dependence. The parameter D>0 (W(°C -) m 2 1 ) is a diffusion coefficient. The diffusive approach comes with the boundary conditions that the gradiant of the temperature profile T n (y) equals zero at the equator and at the North Pole.
We note that Budyko [16] and Sellers [17] each utilized a discrete time approach in their seminal 1969 papers. (Indeed, our discrete approach also conforms with [26], in which it is argued that the annual cycle is the most appropriate time scale at which to model outgoing longwave radiation linearly as a function of T.) We also note Budyko used an alternative formulation of the heat transport term, namely, In this approach it is the temperature at latitude y relative to the global average surface temperature that governs the meridional heat transport. Equation (1) leads to consideration of the temperature evolution equation Budyko, Sellers and others working with energy balance climate models focused on the equilibrium temperature profiles *( ) T y of model equations such as (4), and how these solutions varied with a parameter (typically Q). In particular, there was no consideration of an ice line allowed to move with changes in temperature. Not until E. Widiasih's paper [27] was the evolution of surface temperature coupled to a dynamic ice line, thereby introducing a dynamical systems approach to the study of the model. Although Widiasih used heat transport term (3) and worked in the infinite dimensional setting (with parameters aligned with the present climate), the methods employed here are similar in spirit to those used in [27].
We consider the coupled system where  > 0 is a small parameter and T c is a critical temperature. Note that T n (η n )>T c implies η n+1 >η n , so that the ice line moves poleward. If T n (η n )<T c , then η n+1 <η and the ice line descends equatorward. Equilibria of system (5) are pairs * * h ( ( ) ) T y , for which equation (1) equals zero, with the temperature at the ice line additionally satisfying * * h = ( ) T T c . The ice line quation (5b) was introduced in [27]. Due to the large difference between the time scales for changes in the temperature and the ice line, system (5) is an example of a fast-slow system [18].

The albedo function
The albedo function used here is motivated by [15] , in which the extensive glacial episodes of the Neoproterozoic Era were modeled. Utilizing an idealized general circulation model parametrized for these great ice ages, Abbot et al found that evaporation exceeded precipitation in a latitude band situated within the tropics. This leads to the conceptual model assumption that ice forming below a fixed latitude, denoted throughout by ρ, would not have a snow cover, resulting in an albedo for this 'bare' ice lying between the high albedo of snowcovered ice and the low albedo of ice-free surface. There is then an increase in absorbed solar radiation through the bare ice, relative to that passing through the snow covered ice. The resulting increase in local surface temperature provides a possible mechanism by which the advance of the ice sheet might be arrested.
A stable Jormungand equilibrium solution, with the ice line resting in tropical latitudes, was found by Abbot et al when running the large computer model used in [15]. Working in the infinite dimensional setting and using heat transport term (3), a dynamically stable Jormungand equilibrium solution was shown to exist in [28] when using a smooth version of albedo function (6).
Incorporating the lower albedo of bare ice into (5), we fix a latitude ρ as in [15] and define the albedo function as follows: Here, a a a < < i 1 2 represent the albedos of ice-free surface, bare ice, and snow-covered ice, respectively. The use of the superscript + will indicate ρη, while the superscript − will indicate η<ρ, in all that follows.

The model difference equations
The use of Legendre polynomial expansions to analyze the temperature equation (4) is quite common ( [29,[21][22][23][24]), due to the fact the Legendre polynomials P n (y) are eigenfunctions of the spherical diffusion operator: We incorporate finite Legendre expansion approximations in the coupled temperature-ice line model (5). Recall T n (y), s(y) and a h ( ) y, are each even functions of y, given the assumed symmetry of the model about the equator. We thus write for θ=0, π/2 [30], T n (y) satisfies the prescribed boundary conditions for any N.
We also write corresponds to η<ρ, via definition (6). Also note the appearance of ρ as a limit of integration in h Substitution of (8), (9) and (10) into equation (5a) yields Equating the coefficients of ( ) P y Adding the ice line equation (5b) and doing a bit of rearranging, one finds the N+2 piecewise-defined equations , and hence (12) is everywhere continuous. Notably, however, system (12) is not smooth at any point in the set do not agree on Σ. This fact provides for the lack of smoothness of the invariant manifolds to be described below.
For convenience with the ensuing analysis, we extend system (12) to   + N 1 as follows. (For the persistence of smooth, noncompact normally hyperbolic invariant manifolds of bounded geometry, see [31].) for y<0, and by setting for y>1. That is, for all y<0 we set T n (y)=T n (0), while for all y>1 we set T n (y)=T n (1). Note albedo function (7) . We note this system is nonsmooth at any point (x, η) with η=0, ρ or 1. For where we take i=0, 1, ... , N. For ease of notation, set Consider the functions , , , 0, ... , , 16 We therefore are interested in understanding the behavior of orbits for the mapping (so M 0 is parametrized by η; see figure 1). We assume throughout that N is (maximally) chosen so that Assumption (21) guarantees that for ò=0, M 0 is a globally attracting curve of fixed points, with orbits converging to M 0 exponentially fast. In the appendix we use the Hadamard graph transform technique [19] to prove the existence of an attracting invariant manifold  M that is within  ( ) O of M 0 in the sup norm, for sufficiently small  . This allows for the analysis of the behavior of orbits of the function H and provides for an understanding of various bifurcation scenarios, as we demonstrate in the following section.

Existence of a stable Jormungand climate state
The proof of the following theorem can be found in the appendix.  Before discussing the model dynamics, we pause to comment on the choices for parameters presented in table 1 and used for the plots below. Generally, parameters utilized in a conceptual model of a system as complex as planetary climate will not be well-constrained. As mentioned previously, the parameters A and B used in the outgoing longwave radiation (OLR) term have been estimated as functions of surface temperature using satellite data [26]. During the extremely cold climate of the Neoproterozoic Era ice ages under consideration here, the drawdown of atmospheric CO 2 via silicate weathering would be greatly reduced. The associated buildup of atmospheric CO 2 would serve to decrease OLR, leading to a smaller value for A relative to that provided in [26] for the present climate. The value A=164 Wm −2 in table 1 is in line with values used in [15].
From a conceptual modeling perspective, the diffusion constant appropriate for the present climate is D≈0.5 W/(°C m 2 ) [33]. Given that the cold climates discussed here may reduce the efficiency of the heat transport, we use D≈0.25, as indicated in [33]. Assuming the obliquity was not significantly different 600-700 mya, we use the current value β=23.4°, for which the coefficients = s i , 0, 1, ... ,5 i 2 in table 1 and arising from expansion (9) are provided in [34].
The albedo values are difficult to estimate ( With parameters as in table 1, we plot z(η) in (22) for N=1, ... , 5 in figure 3. correspond to attracting fixed points for  f h ( ), assuming  > 0 is sufficiently small, since . Hence mapping H (19) will have an attracting fixed point on  M with η-value near any 0, H will have an unstable fixed point with 1-dimensional unstable manifold and (N+1)-dimensional stable manifold.
As for the dynamics on  M more globally, the behavior of η restricted to  M is given by iteration of the function is a homeomorphism for  sufficiently small, from which it follows all  p -orbits are monotonic, also ensues from the analysis in the appendix. Hence for ,η is either fixed, strictly increasing, or strictly decreasing under iteration of H.
We conclude that  p will have three fixed points η 1 <η 2 <η 3 near the three respective zeroes of z(η) for N>1 (see figures 3(b)-(e)). Starting with an η>η 2 , the ice line will tend to the stable fixed point with a smaller ice cap, for which η=η 3 . For initial conditions with η<η 2 , the ice line will tend to the stable fixed point with a large ice cap (η=η 1 ). We note, in particular, that the model admits a stable Jormungand equilibrium state with the ice line resting in tropical latitudes below ρ for each of N=1, ... , 5. We also remark that both an ice-free Earth and a snowball Earth are unstable in the sense that  < ( ) p 1 0(so an ice line would descend to η=η 3 ), and  > ( ) p 0 0(so the ice line would retreat to η=η 1 ), for N>1.

Local bifurcations
We fix N=5 in this section, for which z(η) is plotted in figure 3(e). We note the parameter A can be thought of as a proxy for atmospheric greenhouse gas concentrations, in the following sense. A large A-value corresponds to a large OLR term, which in turn corresponds to a low value of, say, atmospheric CO 2 concentration. A smaller Avalue leads to a small OLR term, and hence a larger atmospheric concentration of CO 2 .
, the corresponding fixed point of H on  M is stable (resp., unstable).

Note the parameter A appears in
h ( ) f i 2 only when i=0. Hence A appears precisely once in the expression for z(η) in (22). While we omit the details, setting z(η)=0 and solving for A yields a correspondence between fixed points of  f and A, that is, between the position of the ice line at equilibrium and the parameter A. For sufficiently small A values-or large CO 2 concentrations-there are no equilibria, with the ice line increasing for all time. The physical interpretation is that the climate is heading to an ice-free state. A saddlenode bifurcation occurs at A≈153 Wm −2 , with a small unstable ice cap and a moderately sized stable ice cap equilibria pair appearing. In an analysis using Legendre polynomial expansions to approximate the temperature equation (4), but in which the ice line was held stationary, North argued for the existence of the unstable small ice cap equilibrium state in [23].
The stable Jormungand state comes into existence as A passes through a nonsmooth fold bifurcation [36] at roughly 159 Wm −2 . Note the system now admits four fixed points, provided by two pair of stable/unstable fixed points. The stable moderately sized ice cap equilibrium disappears in a saddle node bifurcation at A≈166 Wm −2 . As the atmospheric CO 2 concentrations continue to decrease-and the OLR term increases, cooling the planet-the Jormungand state is the sole remaining stable equilibrium.
The model exhibits an unstable equilibrium with ice line near the equator for A slightly larger than 175 Wm −2 . There is one final saddle-node bifurcation at A≈181 Wm −2 in which the Jormungand state is lost; the OLR is sufficiently large so as to dominate the warming effect the bare ice had on the absorption of insolation. If A is large enough there are no equilibria, with the ice line heading toward the equator in a runaway snowball Earth episode.
We note the range of A-values for which a stable Jormungand state exists is as large as that obtained in [15] when using heat transport formulation (3). That the Jormungand state is stable for such a wide range of CO 2 values is notable in that the Neoproterozoic Era ice ages under consideration here lasted several million years, evidently exhibiting an insensitivity to changes in CO 2 .
Fixing A=164 Wm −2 , we plot the equilibrium ice line position as a function of the diffusion coefficient D in figure 5. While there is a significant range of D-values for which the stable Jormungand state exists, the model exhibits bistability over much of this range. The Jormungand state is the sole stable equilibrium only for Dvalues roughly between 0.35 and 0.44. For very efficient heat transport (to the right in figure 5), the climate system tends toward a 'uniform' climate of either a warm ice-free world or a snowball Earth, depending on whether the ice line starts above or below an unstable equilibrium position. For highly inefficient heat transport (toward the left in figure 5), the warm tropical climate and the cold northern regions have little interaction, leading to a stable equilibrium solution with the ice line lying between 0.4 and 0.6.

Conclusion
The zonally averaged surface temperature models of M Budyko and W Sellers are energy balance models consisting of terms representing incoming solar radiation, outgoing longwave radiation, and meridional heat transport. These conceptual models were introduced to investigate the role played by positive ice albedo feedback in influencing climate. E Widiasih incorporated a dynamical systems approach by coupling the temperature profile equation to an ice line allowed to respond to changes in temperature. Working with Budyko's heat transport term (3) and in the infinite dimensional setting, with parameters aligned with the present climate, Widiasih showed there exists a stable equilibrium temperature function-ice line pair , with * h positioned roughly as it is today. She also showed there exists a second, unstable equilibrium solution with a large ice cap. While not equilibria, Widiasih argued that a snowball Earth state is stable, while the ice-free state is unstable. The heat transport term (3) used by Budyko was replaced with the diffusive transport term (2), and the resulting equation coupled to Widiasih's ice line equation, in [37]. Using the spectral method, and with parameters aligned with the present climate, results mirroring those of Widiasih were obtained.
Controversy continues to shroud the proposition that snowball Earth events have occurred in our planet's past. A means by which the effect of positive ice albedo feedback might be lessened for sufficiently large ice sheets, so that a strip of ocean water about the equator remained ice-free, was presented in [15]. Incorporating this mechanism into the Budyko-Widiasih model with heat transport term (3), and working in the infinite dimensional setting, a stable Jormungand equilibrium solution with the ice line in tropical latitudes was shown to exist in [28]. It is of interest to note a stable Jormungand state was found when running an idealized general circulation model in [15], parametrized for the extremely cold world associated with an extensive ice age.
We have incorporated the diffusive heat transport term (2), along with changes in surface albedo motivated by [15], into the surface temperature-ice line model in this work. Assuming a discrete time approach and using the spectral method, the analysis results in a nonsmooth singular perturbation problem. We proved an invariant manifold persists under perturbation, assuming the ice sheet moves sufficiently slowly relative to the evolution of temperature. This result was used to show the model exhibits a stable Jormungand equilibrium state, with appropriately chosen parameters. The model displays a rich bifurcation structure, both as atmospheric CO 2 concentrations vary and with changes in the efficiency of the meridional heat transport.
The dynamics of the model analyzed here are determined by a single equation governing the movement of the ice line, for  > 0 and sufficiently small. Coupling this one equation with either a dynamic equation for atmospheric CO 2 concentrations (as in [8], where A varies with time), or with a latitudinal dependence assumed for the diffusion constant (as suggested in [22]), would provide for further lines of inquiry.   Note we used s(0) and 1 as upper bounds for s(η) and h ( ) q i 2 , respectively, on [0, 1]. For 0ηρ, a similar computation yields Suppose 0<η 1 <ρ<η 2 <1. As h 0 is not smooth at points for which η=ρ, the derivative cannot be used to bound Lip(h 0 ) in this case. Alternatively, and using One easily checks that L 0 serves as an upper bound for or 1<η 2 as well, given the simplicity of the extension (15) of (12). Hence, Lip(h 0 )L 0 . , We now introduce various quantities that will appear in ensuing estimates, as well as the function space to be used for the graph transform. Recall that L 0 (23) is a Lipschitz constant for h 0 (η). Let We note d1 and, given the parameters in table 1, L 0 80, so that L80. Let , 0, ... , , and .
We will make use of various inequalities involving the expressions above, including the following.
We also have with the latter inequality following from   g = We will show the desired invariant manifold  M is the graph of a function * h ( ) g , for  > 0 and sufficiently small. To that end, we utilize the function spaces is continuous and , We are now in a position to prove the existence of the invariant manifold  M .
Theorem 2. For  > 0 and sufficiently small, there exists *  Î g L such that is a locally attracting invariant manifold for mapping (19). In addition,  M lies within  ( ) O of M 0 in the sup norm.
That theorem 2 holds true follows from the following propositions 3-7.
and, for ease of notation, let u=w+b 2 (w), v=w+b 1 (w). We have Since w was arbitrary, we have Given our assumed bound on  , we have shown Φ is a contraction map on the complete space ( completes the proof of proposition 3. , Remark. We note b in proposition 3 depends upon  Î g L , given that the function * b in the proof of proposition 3 depends on g . As in [27], we call b the preimage of the ice line h. The next step in the proof of theorem 2 is to produce a graph that is invariant under the function H(x, η) (Propositions 4-6). To that end, consider the mapping , , , where β is the preimage of the ice line η whose existence is guaranteed by proposition 3. We will show that Γ fixes a function *  Î g L , and that the graph of g * provides an invariant manifold  M for H(x, η), provided ò is sufficiently small.