Nonsmooth Invariant Manifolds in a Conceptual Climate Model

TL;DR

This paper investigates a conceptual climate model to understand the stability of tropical ice sheets and the conditions under which open water can persist at the equator during ice ages, using advanced mathematical techniques.

Contribution

It introduces a novel analysis of a coupled temperature-ice line model, proving the existence of a stable tropical ice line state through nonsmooth invariant manifolds.

Findings

Existence of a stable tropical ice line state with open water at the equator.

Identification of bifurcations related to CO2 levels and heat transport efficiency.

Application of spectral and Hadamard graph transform methods to analyze nonsmooth dynamics.

Abstract

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 CO$_2$ concentrations and the efficiency of meridional heat transport, respectively, are varied.

Tables

Table 1. Table 1: Values of parameters and constants.

Parameter	Value	Units	Parameter	Value	Units
$R$	20	Wyr${(}^{\circ }$C m${}_{}{}^2){}^{-1}$	${\alpha }_i$	0.4	dimensionless
$Q$	321	Wm-2	${\alpha }_2$	0.8	dimensionless
$A$	164	Wm-2	$s_0$	1	dimensionless
$B$	1.9	W${(}^{\circ }$C m${}_{}{}^2){}^{-1}$	$s_2$	-0.477131	dimensionless
$D$	0.25	W${(}^{\circ }$C m${}_{}{}^2){}^{-1}$	$s_4$	-0.045029	dimensionless
$\beta$	23.4	degrees	$s_6$	0.007937	dimensionless
$T_c$	0	∘C	$s_8$	0.013859	dimensionless
${\alpha }_1$	0.30	dimensionless	$s_{10}$	0.008663	dimensionless

Equations

$$R\cdot \frac{T_{n+1}(y)-T_n(y)}{(n+1)-n}=Qs(y)(1-\alpha (y,\eta ))-(A+B{T}_n(y))+D\frac{\partial }{\partial y}(1-y^2)\frac{\partial T_n}{\partial y}.$$

$$s(y)=\frac{2}{{\pi }^2}{\int }_0^{2\pi }\sqrt{1-(\sqrt{1-y^2}\sin \beta \cos \gamma -y\cos \beta {)}^2}d\gamma ,$$

$$D\frac{\partial }{\partial y}(1-y^2)\frac{\partial T}{\partial y}$$

$$-C\left(T(y)-{\int }_0^{1}T(y)dy\right),C>0.$$

$$T_{n+1}(y)=T_n(y)+\frac{1}{R}\left(Qs(y)(1-\alpha (y,\eta ))-(A+B{T}_n(y))+D\frac{\partial }{\partial y}(1-y^2)\frac{\partial T_n}{\partial y}\right).$$

$$T_{n+1}(y)$$

$${\eta }_{n+1}$$

$$\alpha (y,\eta )=\{\begin{array}{ll}{\alpha }^{-}(y,\eta ),& \eta <\rho \\ {\alpha }^{+}(y,\eta ),& \rho \le \eta ,\end{array}$$

$$\begin{array}[]{cc}\alpha^{-}(y,\eta)=\begin{cases}\alpha_{1},&y<\eta\\ \frac{1}{2}(\alpha_{1}+\alpha_{i}),&y=\eta\\ \alpha_{i},&\eta<y<\rho\\ \frac{1}{2}(\alpha_{i}+\alpha_{2}),&y=\rho\\ \alpha_{2},&\rho<y,\end{cases}&\hskip 7.22743pt\mbox{and}\hskip 10.84006pt\alpha^{+}(y,\eta)=\begin{cases}\alpha_{1},&y<\eta\\ \frac{1}{2}(\alpha_{1}+\alpha_{2}),&y=\eta\\ \alpha_{2},&y>\eta.\end{cases}\end{array}$$

$$\frac{d}{dy}(1-y^2)\frac{d}{dy}{p}_n(y)=-n(n+1)p_n(y),n=0,1,....$$

$$T_n(y)=\sum _{i=0}^N{x}_{n,2i}{p}_{2i}(y),$$

$$s(y)=\sum _{i=0}^N{s}_{2i}{p}_{2i}(y),$$

$$s_{2i}=(4i+1){\int }_0^{1}s(y)p_{2i}(y)dy.$$

$$s(y){\alpha }^{\pm }(y,\eta )=\sum _{i=0}^N{a}_{2i}^{\pm }(\eta )p_{2i}(y)dy,$$

$$a_{2i}^{+}(\eta )$$

$$a_{2i}^{-}(\eta )$$

$$={\alpha }_2{s}_{2i}-(4i+1)\left(({\alpha }_2-{\alpha }_i){\int }_{\eta }^{\rho }s(y)p_{2i}(y)dy+({\alpha }_2-{\alpha }_1){\int }_0^{\eta }s(y)p_{2i}(y)dy\right).$$

$${\sum }_{i=0}^N{x}_{n+1,2i}{p}_{2i}(y)$$

$$-\left(A+B{\sum }_{i=0}^N{x}_{n,2i}{p}_{2i}(y)\right)-D{\sum }_{i=0}^{N}2i(2i+1)x_{n,2i}{p}_{2i}(y)].$$

$$x_{n+1,0}$$

$$x_{n+1,2i}$$

$${\eta }_{n+1}$$

$$f_0^{\pm }(\eta )=\frac{1}{B}(Q(s_0-a_0^{\pm }(\eta ))-A),f_{2i}^{\pm }(\eta )=\frac{1}{B+2i(2i+1)D}Q(s_{2i}-a_{2i}^{\pm }(\eta )),i\ge 1.$$

$$\Sigma =\{(x,\eta ):x\in {\mathbb{R}}^{N+1},\eta =\rho \};$$

$$\begin{array}[]{cc}f_{2i}(\eta)=\begin{cases}f^{-}_{2i}(0),&\eta<0\\ f^{-}_{2i}(\eta),&0\leq\eta<\rho\\ f^{+}_{2i}(\eta),&\rho\leq\eta\leq 1\\ f^{+}_{2i}(1),&1<\eta,\end{cases}&\hskip 14.45377pt\mbox{and}\hskip 14.45377ptq_{2i}(\eta)=\begin{cases}p_{2i}(0),&\eta<0\\ p_{2i}(\eta),&0\leq\eta\leq 1\\ p_{2i}(1),&1<\eta.\end{cases}\end{array}$$

$$x_{n+1,2i}$$

$${\eta }_{n+1}$$

$${\gamma }_i=\frac{B+2i(2i+1)D}{R},i=0,...,N.$$

$$F:{\mathbb{R}}^{N+1}\times \mathbb{R}\to {\mathbb{R}}^{N+1},(F(x,\eta ){)}_i=-{\gamma }_i((x{)}_i-f_{2i}(\eta )),i=0,...,N,$$

$$G:{\mathbb{R}}^{N+1}\times \mathbb{R}\to \mathbb{R},G(x,\eta )=ϵ\left({\sum }_{i=0}^N(x{)}_i{q}_{2i}(\eta )-T_c\right).$$

$$(X_{n+1},{\eta }_{n+1})=(X_n+F(X_n,{\eta }_n),{\eta }_n+G(X_n,{\eta }_n)).$$

$$H:{\mathbb{R}}^{N+1}\times \mathbb{R}\to {\mathbb{R}}^{N+1}\times \mathbb{R},H(x,\eta )=(x+F(x,\eta ),\eta +G(x,\eta )).$$

Full text

Nonsmooth Invariant Manifolds in a Conceptual Climate Model

James Walsh

Department of Mathematics, Oberlin College

(May 8, 2017)

Abstract.  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.

1 Introduction

The theory of nonsmooth dynamical systems is playing an ever expanding role in the analysis of conceptual climate models. Nonsmooth models arising in the study of the overturning ocean circulation include [16], [29], [30] and [34]. Several works have appeared concerning nonsmooth models of sea ice loss in the Arctic ([7], [8], [11], 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 [2] or the release of carbon dioxide from the deep ocean [13]; with changes in Antarctic ice volume [24]; or with changes in glacial ablation rates [33].

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 [19] 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 (mya) ago, respectively, ice sheets flowed into the ocean at sea level in tropical latitudes [3]. Many believe these extensive glaciations were snowball Earth episodes, with ice covering the entire planet, triggered in part by positive ice albedo feedback (see [12] 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 [1] and references therein). The actual extent of the ice cover during these ice ages remains controversial [3].

Conceptual climate models, focusing on notions of energy balance, were introduced by M. Budyko [5] and W. Sellers [28] 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 [12], 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 [1]. For example, photosynthetic organisms and other life forms are know 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 [1]. 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 [1], after the sea serpent from Norse mythology.

We incorporate the assumptions from [1] 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 ([14] 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 Hadamard’s graph transform method.

The components of the model, including an albedo function designed for the Neoproterozoic ice ages under consideration and following that presented in [1], are introduced in Section 2. In Section 3 the piecewise-defined 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 $\epsilon$ representing the time constant for movement of the ice sheet, in Section 4. In addition, the unperturbed invariant manifold is shown to be Lipschitz continuous and the function spaces to be used for the graph transforms are introduced in Section 4. In Section 5 the existence of an attracting invariant manifold is proved to exist, for sufficiently small $\epsilon$, and in Section 6 we present the model dynamics and we discuss local smooth and nonsmooth bifurcations. We provide concluding remarks in Section 7.

2 The temperature–ice line model

Let $T_{n}(y)\ (^{\circ}$C) denote the average annual surface temperature at $y=\sin\theta$, where $\theta$ 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 [6], [20], [21], [22] and [23]. Meridional heat transport refers to the heat flux across latitude circles associated with atmospheric and oceanic currents.

Consider the difference equation

[TABLE]

The model assumes a symmetry across the equator, so we take $y\in[0,1]$. The units on each side of (1) are Wm*-2*.  $R$ is the heat capacity of the Earth’s surface, with units Wyr$(^{\circ}$C m${}^{2})^{-1}$, while $Q$ (Wm*-2*) is the average annual insolation for the entire Earth. The function

[TABLE]

in which $\beta$ denotes the tilt (or obliquity) of the Earth’s spin axis, serves to incorporate insolation as a function of latitude [17]. 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 $y=\eta$ and called the ice line, with no ice at latitudes equatorward of $y=\eta$. The function $\alpha(y,\eta)$ represents the surface albedo, which depends upon the position of the ice line. Hence $Qs(y)(1-\alpha(y,\eta))$ 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 [10]. The diffusive meridional heat transport term is assumed to be proportional to the surface temperature gradient, and

[TABLE]

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 [5] and Sellers [28] each utilized a discrete time approach in their seminal 1969 papers. We also note Budyko used an alternative formulation of the heat transport term, namely,

[TABLE]

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

[TABLE]

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 [35] 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 used here are similar in spirit to those in [35].

We consider the coupled system

[TABLE]

where $\epsilon>0$ is a small parameter and $T_{c}$ is a critical temperature. Note that $T_{n}(\eta_{n})>T_{c}$ implies $\eta_{n+1}>\eta_{n}$, so that the ice line moves poleward. If $T_{n}(\eta_{n})<T_{c}$, then $\eta_{n+1}<\eta$ and the ice line descends equatorward. Equilibria of system (5) are pairs $(T^{*}(y),\eta^{*})$ for which equation (1) equals zero, with the temperature at the ice line additionally satisfying $T^{*}(\eta^{*})=T_{c}$. The ice line quation (5b) was introduced in [35].

2.1 The albedo function

The albedo function used here is motivated by [1] , 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 $y=\rho$ would not have a snow cover, resulting in an albedo for this “bare” ice lying between the high albedo of snow covered 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 [1]. Working in the infinite dimensional setting and using heat transport term (3), a dynamically stable Jormungand equilibrium solution was shown to exist in [32] when using a smooth version of albedo function (6).

Incorporating the lower albedo of bare ice into (5), we fix a latitude $y=\rho$ as in [1] and we define the albedo function as follows:

[TABLE]

where

[TABLE]

Here, $\alpha_{1}<\alpha_{i}<\alpha_{2}$ represent the albedos of ice free surface, bare ice, and snow covered ice, respectively. The use of the superscript + will indicate $\rho\leq\eta$, while the superscript - will indicate $\eta<\rho$, in all that follows.

3 The model difference equations

The use of Legendre polynomial expansions to analyze the temperature equation (4) is quite common ([15], [20], [21], [22], [23]), due to the fact the Legendre polynomials $p_{n}(y)$ are eigenfunctions of the spherical diffusion operator:

[TABLE]

We incorporate finite Legendre expansion approximations in the coupled temperature-ice line model (5). Recall $T_{n}(y),s(y)$ and $\alpha(y,\eta)$ are each even functions of $y$, given the assumed symmetry of the model about the equator. We thus write

[TABLE]

with $x_{n,2i}$ undetermined coefficients, $i=0,...,N$. As the even Legendre polynomials each have zero gradient at the equator and at the North Pole [21], $T_{n}(y)$ satisfies the prescribed boundary conditions for any $N$.

We also write

[TABLE]

where

[TABLE]

Finally, we express

[TABLE]

where

[TABLE]

and

[TABLE]

Again, $a^{+}_{2i}(\eta)$ corresponds to $\rho\leq\eta$, while $a^{-}_{2i}(\eta)$ corresponds to $\eta<\rho$, via definition (6).

Substitution of (8), (9) and (10) into equation (5a) yields

[TABLE]

Equating the coefficients of $p_{2i}(y)$ in (11), we have

[TABLE]

Adding the ice line equation (5b) and doing a bit of rearranging, one finds the $N+2$ piecewise-defined equations

[TABLE]

where

[TABLE]

(Recall $f^{+}_{2i}(\eta)$ corresponds to $\rho\leq\eta$, and $f^{-}_{2i}(\eta)$ corresponds to $\eta<\rho$.) We note each $f^{\pm}_{2i}(\eta)$ is a polynomial of degree $2N+2i+1$.

The state space for (12) is $\mathbb{R}^{N+1}\times[0,1]$. We observe that $a^{+}_{2i}(\rho)=a^{-}_{2i}(\rho),$ so that $f^{+}_{2i}(\rho)=f^{-}_{2i}(\rho),\,i=0,...,N$, and hence (12) is everywhere continuous. Notably, however, system (12) is not smooth at any point in the set

[TABLE]

for example, $\frac{\partial f^{-}_{0}}{\partial\eta}=\frac{Q}{B}(\alpha_{i}-\alpha_{1})s(\eta)$ and $\frac{\partial f^{+}_{0}}{\partial\eta}=\frac{Q}{B}(\alpha_{2}-\alpha_{1})s(\eta)$ do not agree on $\Sigma.$ 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 $\mathbb{R}^{N+1}\times\mathbb{R}$ as follows. (For the persistence of smooth, noncompact normally hyperbolic invariant manifolds of bounded geometry, see [9].) Extend expression (8) to $y\in\mathbb{R}$ by setting $T_{n}(y)=\sum^{N}_{i=0}x_{n,2i}p_{2i}(0)$ for $y<0$, and by setting $T_{n}(y)=\sum^{N}_{i=0}x_{n,2i}p_{2i}(1)$ 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) is well defined for $y,\eta\in\mathbb{R}$. For $i=0,...,N$, set

[TABLE]

These choices of $f_{2i}(\eta)$ and $q_{2i}(\eta)$ ensure that

[TABLE]

is a continuous system of difference equations for $(x,\eta)\in\mathbb{R}^{N+1}\times\mathbb{R}$. We note this system is nonsmooth at any point $(x,\eta)$ with $\eta=0,\rho$ or 1.

For $v\in\mathbb{R}^{N+1}$, let $(v)_{i}$ denote the $i$th coordinate of $v$, where we take $i=0,1,...,N$. For ease of notation, set

[TABLE]

Consider the functions

[TABLE]

and

[TABLE]

Letting $X_{n}=(x_{n,0},...,x_{n,2N}),$ system (15) is equivalent to

[TABLE]

We therefore are interested in understanding the behavior of orbits for the mapping

[TABLE]

4 The case $\boldsymbol{\epsilon=0}$

4.1 The invariant manifold is Lipschitz continuous

The position of the ice line remains fixed for all time when $\epsilon=0$. In this case (15) admits a manifold $M_{0}=\mbox{graph}(h^{0})$ of fixed points, where

[TABLE]

(see Figure 1). We assume throughout that $N$ is (maximally) chosen so that

[TABLE]

Note (21) implies $|1-\gamma_{i}|\leq 1-\gamma_{0}<1$ for $i=0,...,N$. With $R,B$ and $D$ values as in Table 1, $N=5$ satisfies requirement (21), ensuring that expansion (8) provides a good approximation to the temperature distribution function $T_{n}(y)$. See Figure 2, for example, in which the equilibrium temperature profiles

[TABLE]

are plotted when $\epsilon=0$ for $\eta=0.3$ and $N=4,5$.

Assumption (21) guarantees that for $\epsilon=0$,   $M_{0}$ is a globally attracting curve of fixed points. As noted above, the existing theory for the persistence of normally hyperbolic invariant manifolds under perturbation cannot be applied due to the lack of smoothness of $M_{0}$. We will use Hadamard’s graph transform technique to prove the existence of an attracting invariant manifold $M_{\epsilon}$ that is within $O(\epsilon)$ of $M_{0}$, for sufficiently small $\epsilon$. This will allow for the analysis of the behavior of orbits of the function $H(x,\eta)$ (19), and provide for an understanding of bifurcations for system (15).

We first note that, while nonsmooth, $M_{0}$ is the graph of a Lipschitz continuous function.

Proposition 1. The function $h^{0}(\eta)$ is Lipschitz continuous.

Proof.  Let $h^{0}_{\pm}(\eta)=(f^{\pm}_{0}(\eta),...,f^{\pm}_{2N}(\eta))$. For $\rho\leq\eta\leq 1,$ we have

[TABLE]

[TABLE]

Note we used $s(0)$ and 1 as upper bounds for $s(\eta)$ and $q_{2i}(\eta)$, respectively, on $[0,1]$.

For $0\leq\eta\leq\rho,$ a similar computation yields

[TABLE]

Suppose $0<\eta_{1}<\rho<\eta_{2}<1$. As $h^{0}$ is not smooth at points for which $\eta=\rho$, the derivative cannot be used to bound Lip$(h^{0}$) in this case. Alternatively, and using $h^{0}_{-}(\rho)=h^{0}_{+}(\rho)$, we have

[TABLE]

Let

[TABLE]

One easily checks that $L_{0}$ serves as an upper bound for $\|h^{0}(\eta_{2})-h^{0}(\eta_{1})\|/|\eta_{2}-\eta_{1}|$ if $\eta_{1}<0$ or $1<\eta_{2}$ as well, given the simplicity of the extension (15) of (12). Hence, Lip($h^{0})\leq L_{0}.$  □

4.2 Defining the function space ${\cal B}_{\boldsymbol{L}}$

We 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}$ (22) is a Lipschitz constant for $h^{0}(\eta)$. Let

[TABLE]

and set

[TABLE]

We note $d\geq 1$ and, given the parameters in Table 1, $L_{0}\geq 80$, so that $L\geq 80.$ Let

[TABLE]

We will make use of various inequalities involving the expressions above, including the following.

Proposition 2. (a) $-\gamma_{0}+\frac{\gamma_{N}}{d}\leq 0$.

(b) $\gamma_{0}[L((1+\gamma_{N})(N+1)+\gamma_{0}K)]^{-1}\leq[L(N+1+K)]^{-1}\leq\gamma_{0}(N+1)^{-1}$.

(c) If  $0<\epsilon<\gamma_{0}[L((1+\gamma_{N})(N+1)+\gamma_{0}K)]^{-1}$, then

[TABLE]

Proof.  (a) Note $-\gamma_{0}+\frac{\gamma_{N}}{d}\leq 0$ is equivalent to $d\geq\frac{\gamma_{N}}{\gamma_{0}}=1+2N(2N+1)\frac{D}{B}$, which is assumed in (23).

(b) Note $\frac{1+\gamma_{N}}{\gamma_{0}}\geq\frac{1}{\gamma_{0}}=\frac{R}{B}\geq 1$. Hence

[TABLE]

We also have

[TABLE]

with the latter inequality following from $\frac{1}{L}\leq\frac{1}{80}\leq\frac{B}{R}=\gamma_{0}$.

(c) Note the hypothesis of part (c), together with part (b), imply $c>0$. Showing $c<1$ is equivalent to showing

[TABLE]

or

[TABLE]

an expression equivalent to $\epsilon<\gamma_{0}[L((1+\gamma_{N})(N+1)+\gamma_{0}K)]^{-1}$.  □

We will show the desired invariant manifold $M_{\epsilon}$ is the graph of a function $g^{*}(\eta)$, for $\epsilon>0$ and sufficiently small. To that end, we utilize the function spaces

[TABLE]

and

[TABLE]

We turn to the main result in the following section.

5 The invariant manifold $\boldsymbol{M_{\epsilon}}$

In this section we prove the following result.

Theorem 1.  For $\epsilon>0$ and sufficiently small, there exists $g^{*}\in{\cal B}_{L}$ such that $M_{\epsilon}=\mbox{graph}(g^{*})$ is a locally attracting invariant manifold for mapping (19). In addition, $M_{\epsilon}$ lies within $O(\epsilon)$ of $M_{0}$.

That Theorem 1 holds true follows from the following propositions.

Proposition 3. Assume $0<\epsilon<(L(N+1+K))^{-1}$. Let $g\in{\cal B}_{L}$. Then for all $\eta\in\mathbb{R}$, there exists $\beta=\beta(\eta,g)\in\mathbb{R}$ such that

[TABLE]

Proof.  Let $g\in{\cal B}_{L}$.  Consider the mapping

[TABLE]

Let $b_{1},b_{2}\in{\cal B}_{L}$, let $w\in\mathbb{R}$ and, for ease of notation, let $u=w+b_{2}(w),\ v=w+b_{1}(w)$. We have

[TABLE]

Using the fact $g\in{\cal B}_{L}$ and $\|q_{2i}\|_{\infty}=1$, we have, for $i=0,...,N$,

[TABLE]

Since $w$ was arbitrary, we have

[TABLE]

Give our assumed bound on $\epsilon$, we have shown $\Phi$ is a contraction map on the complete space (${\cal B},\|\cdot\|_{\infty})$, implying $\Phi$ has a unique fixed point $b^{*}\in{\cal B}$ satisfying, for all $\eta\in\mathbb{R}$,

[TABLE]

Letting $\beta=\eta+b^{*}(\eta)$ completes the proof.  □

Remark. We note $\beta$ in Proposition 3 depends upon $g\in{\cal B}_{L}$, given that the function $b^{*}$ in the proof of Proposition 3 depends on $g$. As in [35], we call $\beta$ the preimage of the ice line $\eta$.

The next step is to produce a graph that is invariant under the function $H(x,\eta)$. To that end, consider the mapping

[TABLE]

where $\beta$ is the preimage of the ice line $\eta$ whose existence is guaranteed by Proposition 3. We will show that $\Gamma$ fixes a function $g^{*}\in{\cal B}_{L}$, and that the graph of $g^{*}$ provides an invariant manifold $M_{\epsilon}$ for $H(x,\eta)$, provided $\epsilon$ is sufficiently small.

Proposition 4. Assume $0<\epsilon<\gamma_{0}[L((1+\gamma_{N})(N+1)+\gamma_{0}K)]^{-1}$. There exists $c\in[0,1)$ such that for all $g,h\in{\cal B}_{L}$,  $\|\Gamma g-\Gamma h\|_{\infty}\leq c\|g-h\|_{\infty}$.

Proof.  Let $g,h\in{\cal B}_{L}$ and let $\eta\in\mathbb{R}$. Via Proposition 3, there exist $\beta_{1},\beta_{2}\in\mathbb{R}$ such that

[TABLE]

By definition, $(\Gamma g)(\eta)-(\Gamma h)(\eta)=g(\beta_{1})+F(g(\beta_{1}),\beta_{1})-h(\beta_{2})-F(h(\beta_{2}),\beta_{2})$, the $i$th component of which has the form

[TABLE]

$i=0,...,N.$ Recalling assumption (21), we have

[TABLE]

Note that

[TABLE]

Moreover for $i=0,...,N$,

[TABLE]

Returning to (31), we have

[TABLE]

Invoking part (b) of Proposition 2, inequality (32) simplifies to

[TABLE]

Substituting (33) into (30) yields

[TABLE]

Note

[TABLE]

by Proposition 2, part (c). Hence $\|\Gamma g-\Gamma h\|_{\infty}\leq c\|g-h\|_{\infty}$, as desired.  □

We highlight the conditions imposed on $\epsilon$ in Proposition 4:

[TABLE]

Proposition 5. For $\epsilon$ satisfying (34), the mapping $\Gamma:{\cal B}_{L}\rightarrow{\cal B}$ has a unique fixed point $g^{*}\in{\cal B}_{L}$.

Proof.  We show that for all $g\in{\cal B}_{L},\ \Gamma g\in{\cal B}_{L}$, and then use the fact ${\cal B}_{L}$ is a closed subset of the complete space $({\cal B},\|\cdot\|_{\infty})$. The conclusion then follows by invoking Proposition 4.

Let $g\in{\cal B}_{L},\ \eta\in\mathbb{R}$. Letting $\beta\in\mathbb{R}$ denote the preimage of $\eta$ for $g$, as in Proposition 3, we have $(\Gamma g)(\eta)=g(\beta)+F(g(\beta),\beta)$. Note

[TABLE]

It then follows from (35) and assumption (21) that

[TABLE]

using the fact $L\geq dM$. That $1-\gamma_{0}+\textstyle{\frac{\gamma_{N}}{d}}\leq 1$ follows from Proposition 2, part (a). Since $\eta$ was arbitrary, we have $\|\Gamma g\|_{\infty}\leq L$.

Using the fact the Lipschitz constant of the sum of two functions is the sum of the Lipschitz constants, together with (35) and (21), we have

[TABLE]

using the fact $L\geq dL_{0}$. We have shown $\Gamma g\in{\cal B}_{L}$.  □

Proposition 6. For $\epsilon$ satisfying (34), $M_{\epsilon}=\mbox{graph}(g^{*})$ is invariant under the mapping $H$, where $g^{*}\in{\cal B}_{L}$ is the unique fixed point of $\,\Gamma$.

Proof.  Let $\beta\in\mathbb{R}$. Let $\eta\in\mathbb{R}$ satisfy $\beta=\eta+b_{g^{*}}(\eta)$, where $b_{g^{*}}$ is the fixed point of the mapping $\Phi$, corresponding to $g^{*}$, introduced in the proof of Proposition 3. Then $\eta=\beta+G(g^{*}(\beta),\beta)$, and by Proposition 5

[TABLE]

Hence  $H(g^{*}(\beta),\beta)=(g^{*}(\beta)+F(g^{*}(\beta),\beta),\ \beta+G(g^{*}(\beta),\beta))=(g^{*}(\eta),\eta).$  □

Proposition 7. Assume $\epsilon$ satisfies (34) and let $M_{\epsilon}$ be as in Proposition 6. Then

• (a)

$M_{\epsilon}$ is within $O(\epsilon)$ of $M_{0}$, and

• (b)

For $x\in\mathbb{R}^{N+1}$ with $\|x\|\leq L$, and for $\beta\in\mathbb{R}$,  $\|H^{j}(x,\beta)-M_{\epsilon}\|\rightarrow 0$ as $j\rightarrow\infty$, where $H^{j}$ denotes the $j$th iterate of $H$.

Proof.  To prove (a), let $\eta\in\mathbb{R}$ and let $\beta=\eta+b_{g^{*}}(\eta)$, where

[TABLE]

as in Proposition 3. Recall

[TABLE]

By definition of $F$,

[TABLE]

Once again invoking (21), we have

[TABLE]

from which it follows

[TABLE]

where $\omega$ is constant. Hence $M_{\epsilon}$ is within $O(\epsilon)$ of $M_{0}$.

For part (b), let $x\in\mathbb{R}^{N+1}$ with $\|x\|\leq L$. Let $\beta\in\mathbb{R}$, and set $p=(x,\beta)$. Let $g:\mathbb{R}\rightarrow\mathbb{R}^{N+1},\ g(\beta)=x$, so that $p=(g(\beta),\beta).$ Note $g\in{\cal B}_{L}$. Let $q=(g^{*}(\beta),\beta)$. Again as in the proof of Proposition 3, let $r,s\in\mathbb{R}$ satisfy $\beta=r+b_{g^{*}}(r)=s+b_{g}(s)$, that is,

[TABLE]

We have

[TABLE]

Moreover, recalling definitions (16) and (27), one finds

[TABLE]

from which it follows

[TABLE]

We then have

[TABLE]

The assumed bound on $\epsilon$ implies $\epsilon<\frac{\gamma_{0}}{N+1}$ by Proposition 2(b), from which it follows $1-\gamma_{0}+\epsilon(N+1)<1$. This implies $\|H^{j}(p)-M_{\epsilon}\|\rightarrow 0$ exponentially fast as $j\rightarrow\infty$.  □

This completes the proof of Theorem 1. We now return to the model and use Theorem 1 to gain insight into the behavior of orbits of $H$ and, in particular, to show the existence of a stable fixed point for which the ice line lies in tropical latitudes (the Jormungand climate state). We also investigate system bifurcations.

6 Model behavior

6.1 Existence of the Jormungand climate state

Consider the attracting invariant manifold $M_{\epsilon}=\mbox{graph}(g^{*})$ for mapping (19), whose existence follows from Theorem 1. Using the fact $M_{\epsilon}$ is parametrized by $\eta$ and lies within $O(\epsilon)$ of $M_{0}$, the ice line equation on $M_{\epsilon}$ becomes

[TABLE]

(noting $\sum^{N}_{i=0}q_{2i}(\eta)$ is bounded for $\eta\in\mathbb{R}$). Hence the dynamics of $H(x,\eta)$ on $M_{\epsilon}$ are well-approximated by orbits of the function

[TABLE]

which we write

[TABLE]

Iteration of (36) thus allows for the analysis of the behavior of system (15) for $\epsilon$ positive and sufficiently small.

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 [10]. During the extremely cold climate of the Neoproterozoic Era ice ages under consideration here, the drawdown of atmospheric CO2 via silicate weathering would be greatly reduced. The associated buildup of atmospheric CO2 would serve to decrease OLR, leading to a smaller value for $A$ relative to that provided in [10] for the present climate. The value $A=164$ Wm*-2* in Table 1 is in line with values used in [1].

The value $B=1.9$ W$(^{\circ}$C m${}^{2})^{-1}$ is taken from [10]. The critical temperature $T_{c}=0^{\circ}$C is from [1], while $R=20$ Wyr(∘C m2)-1 is taken from [27]. The solar constant is known to have been 94% of its current value 600-700 mya; hence we use $Q=321$ Wm*-2*.

From a conceptual modeling perspective, the diffusion constant appropriate for the present climate is $D\approx 0.5$ W/(∘C m2) [26]. Given that the cold climates discussed here may reduce the efficiency of the heat transport, we use $D\approx 0.25$, as indicated in [26]. Assuming the obliquity was not significantly different 600-700 mya, we use the current value $\beta=23.4^{\circ}$, for which the coefficients $s_{2i},\ i=0,1,...,5$ in Table 1 and arising from expansion (9) are provided in [18].

The albedo values are difficult to estimate ([25], chapter 3). The choice of the albedo of bare ice $\alpha_{i}=0.4$ being slightly greater than that of ice free surface ($\alpha_{1}=0.3$) is important; for sufficiently large $\alpha_{i}$-values the ice albedo feedback dominates, resulting in the loss of the Jormungand solution described below.

With parameters as in Table 1, we plot $z(\eta)$ for $N=1,...,5$ in Figure 3. The function $\phi_{\epsilon}(\eta)$ will have a fixed point $\eta=\eta^{*}$ precisely when $z(\eta^{*})=0$. Moreover, any zeros $\eta=\eta^{*}$ of $z(\eta)$ satisfying $z^{\prime}(\eta^{*})<0$ will correspond to attracting fixed points for $\phi_{\epsilon}(\eta)$, assuming $\epsilon>0$ is sufficiently small, since $\phi^{\prime}_{\epsilon}(\eta^{*})=1+\epsilon z^{\prime}(\eta^{*})$. Hence model equations (15) admit an attracting fixed point having $\eta$-value near any $\eta=\eta^{*}$ for which $z(\eta^{*})=0$ and $z^{\prime}(\eta^{*})<0$. Near any zero $\eta=\eta^{*}$ of $z(\eta)$ with $z^{\prime}(\eta^{*})>0$, system (15) will have an unstable fixed point with 1-dimensional unstable manifold and $(N+1)$-dimensional stable manifold.

As for the dynamics on $M_{\epsilon}$ more globally, note $p_{\epsilon}$ is injective due to Proposition 3. Given that $g^{*}$ and each $q_{2i}$ are bounded for $\eta\in\mathbb{R},\ p_{\epsilon}$ is surjective as well. Thus $p_{\epsilon}:\mathbb{R}\rightarrow\mathbb{R}$ is a homeomorphism, from which it follows all $p_{\epsilon}$-orbits are monotonic. For $\eta$ satisfying $G(g^{*}(\eta),\eta)<0$,   $p_{\epsilon}(\eta)<\eta$, and the corresponding orbit will be decreasing. For $\eta$ satisfying $G(g^{*}(\eta),\eta)>0$,   $p_{\epsilon}(\eta)>\eta$, and the corresponding orbit will be increasing.

For the parameter values in Table 1, $p_{\epsilon}$ will have three fixed points $\eta_{1}<\eta_{2}<\eta_{3}$ near the three respective zeroes of $z(\eta)$ for $N>1$ (see Figure 3(b-e)). Starting with an $\eta>\eta_{2}$, the ice line will tend to the stable fixed point with a smaller ice cap, for which $\eta=\eta_{3}$. For initial conditions with $\eta<\eta_{2}$, the ice line will tend to the stable fixed point with a large ice cap ($\eta=\eta_{1}$). We note, in particular, that the model admits a stable Jormungand equilibrium state with the ice line resting in tropical latitudes below $y=\rho$ 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_{\epsilon}(1)<0$ (so an ice line would descend to $\eta=\eta_{3}$), and $p_{\epsilon}(0)>0$ (so the ice line would retreat to $\eta=\eta_{1}$), for $N>1$.

6.2 Local bifurcations

We fix $N=5$ in this section, for which $z(\eta)$ 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 a large OLR term, which in turn corresponds to a low value of, say, atmospheric CO2 concentration. A smaller $A$-value leads to a small OLR term, and hence a larger atmospheric concentration of CO2.

Note the parameter $A$ appears in $f_{2i}(\eta)$ only when $i=0$. Hence $A$ appears precisely once in expansion (37) for $z(\eta)$. While we omit the details, setting $z(\eta)=0$ and solving for $A$ yields a correspondence between fixed points of $\phi_{\epsilon}$ and $A$, that is, between the position of the ice line at equilibrium and the parameter $A$.

For $\eta\in[0,1]$ the position of the ice line at equilibrium as a function of $A$ is plotted in Figure 4. The solid branches in Figure 4 represent stable fixed points, while dotted branches are comprised of unstable fixed points. The vertical red line in Figure 4 corresponds to the $A$-value used in Figure 3, for which one unstable fixed point was sandwiched between two stable fixed points.

For sufficiently small $A$ values—or large CO2 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 saddle-node bifurcation occurs at $A\approx 153$ Wm*-2*, with a small unstable ice cap and a moderately sized stable ice cap equilibria pair appearing. The existence of the unstable small ice cap equilibrium state was previously shown to occur when using Legendre polynomial expansions to approximate the temperature equation (4) in [22].

The stable Jormungand state comes into existence as $A$ passes through a nonsmooth fold bifurcation [4] 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\approx 166$ Wm*-2*. As the atmospheric CO2 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\approx 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 [1] when using heat transport formulation (3). That the Jormungand state is stable for such a wide range of CO2 values is notable in that the Neoproterozoic Era ice ages under consideration here lasted several million years, evidently exhibiting an insensitivity to changes in CO2.

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 $D$-values 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.

7 Conclusion

M. Budyko’s zonally averaged surface temperature model is an energy balance model comprised of terms representing incoming solar radiation, outgoing longwave radiation, and meridional heat transport. Budyko introduced his conceptual model to investigate the role played by positive ice albedo feedback in influencing climate. E. Widiasih brought a dynamical systems approach to Budyko’s model 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 $(T^{*}(y),\eta^{*})$, with $y=\eta^{*}$ 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 [31]. 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 [1]. 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 [32]. It is of interest to note a stable Jormungand state was found when running an idealized general circulation model in [1], 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 [1], into the Budyko-Widiasih 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 CO2 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 $\epsilon>0$ and sufficiently small. Coupling this one equation with either a dynamic equation for atmospheric CO2 concentrations (as in [2], where $A$ varies with time), or with a latitudinal dependence assumed for the diffusion constant (as suggested in [21]), would provide for further lines of inquiry.

Acknowledgment. The author recognizes and appreciates the support of the Mathematics and Climate Research Network (www.mathclimate.org).