Abstract

Record lows in Arctic sea ice extent have been making frequent headlines in recent years. The change in albedo when sea ice is replaced by open water introduces a nonlinearity that has sparked an ongoing debate about the stability of the Arctic sea ice cover and the possibility of Arctic “tipping points.” Previous studies identified instabilities for a shrinking ice cover in two types of idealized climate models: (i) annual-mean latitudinally varying diffusive energy balance models (EBMs) and (ii) seasonally varying single-column models (SCMs). The instabilities in these low-order models stand in contrast with results from comprehensive global climate models (GCMs), which typically do not simulate any such instability. To help bridge the gap between low-order models and GCMs, an idealized model is developed that includes both latitudinal and seasonal variations. The model reduces to a standard EBM or SCM as limiting cases in the parameter space, thus reconciling the two previous lines of research. It is found that the stability of the ice cover vastly increases with the inclusion of spatial communication via meridional heat transport or a seasonal cycle in solar forcing, being most stable when both are included. If the associated parameters are set to values that correspond to the current climate, the ice retreat is reversible and there is no instability when the climate is warmed. The two parameters have to be reduced by at least a factor of 3 for instability to occur. This implies that the sea ice cover may be substantially more stable than has been suggested in previous idealized modeling studies.

1. Introduction

Arctic sea ice is undergoing a striking, closely monitored, and highly publicized decline. A recurring theme in the debate surrounding this decline is the question of how stable the ice cover is, and specifically whether it can become unstable. This question is often linked to the ice–albedo feedback, which is expected to play a key role in the observed sea ice retreat.

The ice–albedo feedback has been studied since at least the nineteenth century, when Croll (1875) investigated its importance for the climate system, and it has had a major role in climate science since then. The feedback is associated with a nonlinearity in the climate system due to the jump in local albedo between ice-free and ice-covered surface conditions. This nonlinearity has long been expected to affect the stability of the climate system in the sense that it can potentially trigger abrupt transitions between ice-free and ice-covered regimes (e.g., Budyko 1966). The bifurcations associated with such abrupt transitions are often referred to as “tipping points.”

The idea of an irreversible jump from one stable state to another gained momentum when studies using idealized latitudinally varying diffusive energy balance models (EBMs) of the annual-mean equilibrium state of the global climate (Budyko 1969; Sellers 1969) encountered such bistability in realistic parameter regimes (Budyko 1972; Held and Suarez 1974). These early studies were primarily concerned with the possibility of a catastrophic transition to a completely ice-covered planet, typically called the “snowball earth instability.” In addition, however, many studies have shown that small polar ice covers that terminated poleward of about 75° latitude are also typically unstable in EBMs, a feature referred to as the “small ice cap instability” (SICI) [see review in North (1984)]. More recent studies of the SICI in EBMs have, for example, compared this behavior with global climate model (GCM) results (Winton 2008) and identified additional bifurcations when a more complex representation of energy transport is used (Rose and Marshall 2009).

A bistability analogous to the SICI was found by Thorndike (1992) in an idealized single-column model (SCM) of Arctic sea ice and climate. In contrast with EBMs, this model had no representation of spatial variations, but it included seasonal variations. Modeling the sea ice seasonal cycle required representations of thermodynamic processes associated with sea ice thickness changes. More recently, a number of studies using a range of similar SCMs identified bistability and bifurcations associated with the loss of the sea ice cover (Flato and Brown 1996; Björk and Söderkvist 2002; Eisenman 2007; Eisenman and Wettlaufer 2009; Müller-Stoffels and Wackerbauer 2011; Abbot et al. 2011; Moon and Wettlaufer 2011; Eisenman 2012; Moon and Wettlaufer 2012; Müller-Stoffels and Wackerbauer 2012; Björk et al. 2013).

The results of these EBMs and SCMs stand in contrast with results from most comprehensive GCMs, which typically have found no evidence of bistability (e.g., Armour et al. 2011) and simulate a smooth transition from modern conditions to a seasonally ice-free Arctic (Winton 2006; Ridley et al. 2008; Winton 2008; Tietsche et al. 2011; Ridley et al. 2012). In these simulations, the Arctic and Antarctic sea ice areas typically scale approximately linearly with global-mean temperature (Gregory et al. 2002; Winton 2011; Armour et al. 2011), displaying no indication of instability. The loss of the wintertime-only sea ice cover in a very warm climate (i.e., the transition from seasonally ice-free to perennially ice-free conditions) also typically occurs in a smooth and linear fashion in GCMs (Ridley et al. 2008; Winton 2008; Armour et al. 2011). One GCM showed some evidence of nonlinearity during this transition (Winton 2008), but further analysis suggests that the nonlinearity is not sufficient to cause bistability in this case (Li et al. 2013). In other words, even though the ice–albedo feedback is operating in comprehensive GCMs, they simulate strikingly linear sea ice changes. In contrast to comprehensive GCMs, an atmosphere-only GCM with no seasonal cycle was found to simulate bistability resembling the SICI (Langen and Alexeev 2004), and a GCM with simplified atmospheric physics and idealized ocean geometries was found to simulate bistability with the ice edge of the cold state located at considerably lower latitudes than the SICI (Ferreira et al. 2011).

The discrepancy between the instabilities found in idealized models and the smooth ice retreat found in most comprehensive GCMs raises a conundrum: Is the disagreement between the two approaches the result of a fundamental misrepresentation of the underlying physics in GCMs, or is it rather the result of some aspect of the simplifications used in the idealized models? In more general terms, what physical processes dictate whether there are multiple sea ice states under a given forcing? These questions are the focus of this study.

Our approach is essentially to add a standard SCM representation of seasonally varying sea ice thickness to a standard EBM, thereby combining the two classes of modeling approaches. We therefore study a system that is akin to the idealized models discussed above, but it incorporates both spatial and seasonal variations, taking a small step from conventional idealized models toward comprehensive GCMs and ultimately nature. We investigate how the inclusion of both variations influences the stability of the ice cover, focusing specifically on whether or not unstable solutions occur such that warming and cooling the climate leads to a hysteresis loop.

Note that seasonally varying EBMs have been employed previously to study how idealized land–ocean configurations dictate the glaciation of polar regions (North and Coakley 1979; North et al. 1983; Lin and North 1990; Huang and Bowman 1992). Furthermore, several EBM studies have included some representation of a seasonal cycle and sea ice thickness changes (Suarez and Held 1979; Morales Maqueda et al. 1998; Bendtsen 2002; Bitz and Roe 2004). However, these previous studies have not considered how the inclusion of a seasonal cycle affects the stability of the system.

The article is structured as follows: in section 2 we formulate the diffusive energy balance sea ice model. The simulated climate in the default parameter regime and its response to forcing are presented in section 3. Section 4 discusses how the model reduces to a standard EBM and a standard SCM as limiting cases in the parameter space. Section 5 investigates the presence of bistability as parameters are varied. Conclusions are presented in section 6.

2. Model of sea ice and climate

Here we derive the model used in this study, which is an idealized representation of sea ice and climate with seasonal and latitudinal variations in a global domain. The surface is an aquaplanet with an ocean mixed layer that includes sea ice when conditions are sufficiently cold. We consider only zonally uniform climates (Fig. 1).

Fig. 1.

Schematic of the global model of climate and sea ice described in section 2, showing the fluxes included in the model: insolation (yellow), OLR (red), horizontal heat transport (green), ocean heating (dark blue), and vertical heat flux through the ice (blue). The temperature of the ice is given by T at the surface and Tm at the base.

Fig. 1.

Schematic of the global model of climate and sea ice described in section 2, showing the fluxes included in the model: insolation (yellow), OLR (red), horizontal heat transport (green), ocean heating (dark blue), and vertical heat flux through the ice (blue). The temperature of the ice is given by T at the surface and Tm at the base.

a. Representation of model state

SCMs typically compute the seasonally varying sea ice thickness h(t) as a prognostic variable. The surface temperature of the ice T is calculated diagnostically based on h and the surface energy flux. When the ice thickness reaches zero, SCMs often evolve the temperature of an ocean mixed layer of fixed depth (Thorndike 1992, and later studies). This temperature is taken to be vertically uniform in the ocean mixed layer and hence equal to the surface temperature T. To succinctly account for both ice thickness when ice is present and the temperature of the ocean mixed layer when ice is not present, an SCM can evolve the surface enthalpy E(t), which is defined as (Eisenman and Wettlaufer 2009)

 
formula

where Lf is the latent heat of fusion of sea ice and cw is the heat capacity of the ocean mixed layer, which is equal to the product of the ocean mixed layer specific heat capacity, density, and depth. The melting point Tm is approximated to be constant at its freshwater value. The default values for all model parameters are discussed in section 2d and summarized in Table 1.

Table 1.

Default parameter values.

Default parameter values.
Default parameter values.

In contrast to SCMs, standard EBMs compute the equilibrium value of the spatially varying surface temperature T(x). Here, x ≡ sinθ with latitude θ, such that x = 0 at the equator and x = 1 at the North Pole.

The model developed here evolves the surface enthalpy as a function of both time and latitude E(t, x). A standard SCM is represented in each latitudinal grid box, and the grid boxes are coupled via latitudinal diffusion of T, as is done in standard EBMs. The temperature T represents either (i) the temperature of the ocean mixed layer in the absence of ice or (ii) the ice surface temperature when ice is present. Below, we begin by discussing the seasonally varying EBM formulation that will be used to give spatial dependence, and then we add an SCM representation at each latitudinal grid box to evolve sea ice thickness. The reader who is primarily interested in the model results may skip to section 3.

b. EBM formulation with seasonal variations

The time evolution of E(t, x) is determined at each latitude by the net energy flux into the atmospheric column and surface below:

 
formula

which includes fluxes from top-of-atmosphere net solar radiation aS, outgoing longwave radiation (OLR) L, meridional heat transport in the atmosphere D2T, and heat flux into the model domain from the ocean below Fb (Fig. 1). Each of these terms is discussed below. Last, F in Eq. (2) represents a spatially uniform, seasonally constant climate forcing that is varied during some model runs. Increasing F can be interpreted as an increase in atmospheric CO2, for example, analogous to global warming scenarios in GCM simulations.

Note that most EBMs compute only the equilibrium climate state, ∂E/∂t = 0 in Eq. (2), whereas here we consider the time evolution of the system. Since we are considering an aquaplanet, and the ocean mixed layer has an effective heat capacity that is more than an order of magnitude higher than the atmosphere (North and Coakley 1979), we neglect the heat capacity of the atmospheric column. This implies that the vertical temperature profile in the atmosphere at a given time and latitude is fully determined by T.

1) Solar radiation

Following North and Coakley (1979), we use an approximate representation of the insolation,

 
formula

where S0 is the annual-mean insolation at the equator, S1 determines the amplitude of seasonal insolation variations, ω = 2π yr−1 is the annual frequency, and S2 determines the equator-to-pole insolation gradient (Table 1). Note that the spatial dependence of S(t, x) is written here without using Legendre polynomials, which have often been used in the EBM literature as they enable analytical solutions. We choose this equivalent simpler representation because the model presented here is solved numerically.

Some of the insolation is reflected back to space and the rest is absorbed. The fraction of incident solar radiation that is absorbed is called the planetary, or top-of-atmosphere (TOA), coalbedo (i.e., one minus albedo). It depends on factors including the solar zenith angle, clouds, and the presence of reflective ice at the surface (Lian and Cess 1977). Similar to North (1975b), we write the planetary coalbedo as

 
formula

where a0, a2, and ai are empirical parameters (see Table 1). It is the discrete jump in a(x, E) at the ice edge that introduces nonlinearity and therefore the possibility of multiple states. The addition of ice thickness evolution in the next section, however, will add further nonlinearity. Some studies have considered various smoothed albedo transitions between ice and water (e.g., Cahalan and North 1979; Eisenman and Wettlaufer 2009), which we eschew here for simplicity.

2) OLR

Following Budyko (1969), we represent the OLR as a linear function of the surface temperature,

 
formula

where A and B are empirical parameters (Table 1).

3) Heat transport

Following North (1975b), the meridional heat flux in the climate system is approximated to be proportional to the meridional gradient of the surface temperature, thereby influencing the surface temperature evolution like a diffusive process with constant diffusivity D. Because of converging meridians on a spherical earth, meridional diffusion takes the form (North 1975b)

 
formula

The seasonal cycle introduces a small seasonally varying heat transport across the equator. For both conceptual transparency and computational convenience, we approximate this cross-equatorial heat flux to be zero, allowing us to restrict the model domain to a single hemisphere. Neglecting this in the EBM described in this section causes the temperature to change by no more than 0.3 K at any latitude at any time of year. Inserting Eq. (1) with E > 0 (open water) into Eq. (6) at x = 0 gives the equatorial boundary condition

 
formula

At the pole (x = 1), the heat transport in Eq. (6) goes to zero for any distribution of T because of converging meridians, and no boundary condition is necessary.

4) Ocean heat flux

The upward heat flux (Fb) from the deep ocean into the base of the sea ice or ocean mixed layer, which is associated with heat flux convergence within the ocean column below the mixed layer, is crudely approximated to be constant in space and time. We choose a value of Fb appropriate for the Arctic Ocean, where there is net horizontal ocean heat flux convergence (see section 2d). Note that while this choice helps to accurately simulate conditions in the Arctic Ocean, and treating Fb as constant retains simplicity in the model formulation, this makes the global ocean a heat source, albeit a weak one. This heating is ultimately compensated by increasing the default value of A to get an observationally consistent default climate.

Note that in the absence of sea ice, ∂E/∂t = cwT/∂t on the left-hand side of Eq. (2). Hence the model so far is a fairly standard EBM with the addition of time dependence, seasonally varying forcing, and specified ocean heat flux convergence.

c. Addition of SCM physics

Next, we add a representation of sea ice thickness evolution to the model and compute the resulting effect on the surface temperature. This is done by coupling Eq. (2) to an SCM representation of thermodynamic sea ice growth and melt (Fig. 1). Note that energy fluxes are represented in the present model at the top of the atmosphere, as in typical EBMs, in contrast to typical SCMs, which consider surface fluxes.

We approximately follow the SCM developed in Eisenman (2012). This SCM can be directly derived (cf. Eisenman and Wettlaufer 2009) as an approximate representation of the SCM of Maykut and Untersteiner (1971), which represents processes including vertical energy flux between the ice and the atmosphere, vertical conduction of heat within the ice, vertical heat flux between the ice and the ocean below (Fb), basal congelation during the winter growth season, and surface and basal ablation during the summer season.

Vertical heat conduction within the ice is governed by a diffusion equation for the internal ice temperature. Here we assume that the energy associated with phase changes is much larger than the energy associated with temperature changes within the ice. Equivalently, the relevant Stefan number NS ≡ (LfΔh)(cihΔT)−1 (Eisenman 2012), satisfies NS ≫ 1, where ci is the specific heat of the ice. In this approximation, the vertical diffusion gives a linear temperature profile, with temperature varying from TTm at the top surface to Tm at the bottom surface. The heat flux upward through the ice is then equal to k(TmT)/h, where k is the ice thermal conductivity, which we treat as constant.

If the ocean mixed layer cools to Tm, ice begins to grow. As long as ice is present, the temperature of the ocean mixed layer is set to remain constant at Tm. In this regime, the ocean heat flux Fb goes directly into the bottom of the ice without losing heat to the ocean mixed layer.

Hence, when ice is present, the change in energy per unit area is equal to the change in E, and the ice thickness (h = −E/Lf) evolves according to Eq. (2). In other words, Eq. (2) applies equivalently for both ice-covered and ice-free conditions. In the model results, we identify the ice edge latitude θi as the latitude where E changes sign, and we similarly define xi ≡ sinθi.

There are two possible regimes when ice is present, separated by the surface melt threshold: either (i) T < Tm, in which case the net surface flux must be zero, or (ii) T = Tm, in which case surface melt occurs proportionally to the net surface flux. We determine which case applies by computing the surface temperature T0 that would balance the surface energy flux and determining whether it is above or below the melting point. Since we treat the heat capacity of the atmosphere as negligible compared with the ocean mixed layer, the net surface energy flux is equal to the top-of-atmosphere energy flux plus the horizontal energy convergence in the atmosphere. Then T0 is computed from

 
formula

If T0 < Tm, then the surface flux is balanced by a subfreezing surface temperature, and the frozen regime (i) applies, with T = T0. On the other hand, T0 > Tm indicates that the surface flux cannot be balanced by a subfreezing surface temperature, and surface melt (ii) occurs. In this case the surface temperature remains at the melting point, T = Tm, and the surface melt is calculated using a Stefan condition, which states that the ice thinning from surface melt is equal to the net surface flux divided by Lf. Note that here the value of T0 at a given location depends on T at other locations because of the ∇2T term, making this model considerably less straightforward to solve than a typical SCM.

Ice-free surface conditions can be concisely included by adding a third regime for the surface temperature:

 
formula

Equations (2), (8), and (9) determine the model, which represents the meridionally and seasonally varying global surface temperature and allows sea ice of varying thickness to form when the ocean mixed layer temperature reaches the melting point.

d. Default parameter values

The default values of all model parameters are given in Table 1. Values for D, B, and cw are adopted from North and Coakley (1979), where cw corresponds to a 75-m ocean mixed layer depth. The values of D typically used in the literature range from D = 0.4 (Lin and North 1990) to 0.66 W m−2 K−1 (Rose and Marshall 2009). The insolation parameters S0 and S2 are adopted from the Legendre polynomial coefficients in North and Coakley (1979). Here we set the parameter dictating the seasonal amplitude of the solar forcing S1 to be 25% larger than the astronomically based value used in North and Coakley (1979) in order to better match the observed seasonal range of sea ice thickness in the central Arctic. Note that this increase in S1 may help account for factors such as seasonally varying cloud cover and water vapor that are not included in the model. Since North and Coakley (1979) do not consider an explicit albedo jump, we adopt coalbedo parameter values (a0, a2, and ai) from North (1975b), converting the Legendre polynomial coefficients into our formalism.

For the heat flux upward from the ocean Fb, Maykut and Untersteiner (1971) use a value of Fb = 2 W m−2 based on central Arctic observations. Thorndike (1992) adopts a possible range for Fb of 0–10 W m−2. Other observations suggest an annual-mean value of Fb = 5 W m−2 (Maykut and McPhee 1995). We adopt a value within the range of these previous estimates, Fb = 4 W m−2.

We use values for the thermal conductivity k and latent heat of fusion Lf corresponding to pure ice.

Last, we let A = 193 W m−2, similar to the value in Mengel et al. (1988), in order to simulate an annual-mean hemispheric-mean temperature consistent with the observed present-day climate when F = 0.

In the model simulations presented in this study, we vary the climate forcing F, as well as D and S1. For notational convenience, we denote the default values of the latter two parameters as D* and (Table 1).

e. Numerical solution

Numerically solving the system of Eqs. (2), (8), and (9) is not straightforward because of the implicit representation of T0 in Eqs. (8) and (9). This represents a nonlinear ordinary differential equation in x that involves a free boundary between melting and freezing ice surfaces that would need to be solved at each time step in a standard forward Euler time-stepping routine. Furthermore, using forward Euler time-stepping for this diffusive system would require very high temporal resolution for numerical stability.

To efficiently integrate the model, we instead employ a numerical approach that solves a system of equations that can be visualized in terms of two separate layers. Horizontal diffusion occurs in a “ghost layer” with heat capacity cg, all other processes occur in the main layer, and the temperature of the ghost layer is relaxed toward the temperature of the main layer with time scale τg. This method bypasses the need to solve a differential equation in x at each time step and allows us to use an implicit time-stepping routine with coarse time resolution in the ghost layer. We emphasize that the ghost layer does not represent a separate physical layer such as the atmosphere, which would add physical complexity to the model. This “two layer” system reduces to the system of Eqs. (2), (8), and (9) above in the limit that cg → 0 and τg → 0. In the model simulations, cg and τg are chosen to be sufficiently small that further reducing their values does not substantially influence the numerical solution (values given in Table 1). Further details are given in  appendix A.

The model is solved numerically in 0 ≤ x ≤ 1 using 400 meridional grid boxes that are equally spaced in x and 1000 time steps per year. The spatial resolution is required to sufficiently resolve the evolution of the sea ice cover, especially in parameter regimes where the SICI region is small. The temporal resolution is chosen to ensure numerical stability of the solution.

Code to numerically solve the present model, as well as the standard EBM described in section 2b but with S1 = 0, is available online at http://eisenman.ucsd.edu/code.html.

3. Model results in the default parameter regime

In this section we discuss the simulated climate in the parameter regime (D = D* and S1 = ), first with F = 0, then in the case where the climate is warmed by increasing F, and finally when F is ramped back down.

a. Simulated modern climate (F = 0)

The simulated climate with all parameters at their default values (Table 1) and F = 0 is illustrated in Fig. 2, where the seasonal cycle of the equilibrated climate is plotted. Figure 2a shows the seasonal cycle of E(t, x), which fully represents the model state since E is the only prognostic variable and the forcing varies seasonally. The associated surface temperature (Fig. 2b) and ice thickness (Fig. 2c) are roughly consistent with present-day climate observations in the Northern Hemisphere. The simulated surface temperature varies approximately parabolic with x from an equatorial temperature of about 30°C throughout the year to a polar temperature of −10°C in winter and −3°C in summer (Fig. 2d). There is perennial sea ice around the pole (Fig. 2e) with ice thickness at the pole varying seasonally between 3.1 and 3.4 m (Fig. 2f). The edge of the sea ice cover ranges from θi = 58° in winter to 76° in summer (Fig. 2g), which can be compared with the observed zonal-mean ice edge latitude in the Northern Hemisphere (Eisenman 2010).

Fig. 2.

Simulated climate in the default parameter regime. Contour plot of the seasonal cycle of (a) surface enthalpy E(x, t), (b) surface temperature T(x, t), and (c) sea ice thickness h(x, t). The black curve in (a)–(c) indicates the ice edge. (d) Surface temperature T in summer and winter, corresponding to dashed and solid vertical lines in (c). (e) Ice thickness h in summer and winter where x > 0.7. (f) Seasonal cycle of ice thickness at the pole hp. (g) Seasonal cycle of the latitude of the sea ice edge θi.

Fig. 2.

Simulated climate in the default parameter regime. Contour plot of the seasonal cycle of (a) surface enthalpy E(x, t), (b) surface temperature T(x, t), and (c) sea ice thickness h(x, t). The black curve in (a)–(c) indicates the ice edge. (d) Surface temperature T in summer and winter, corresponding to dashed and solid vertical lines in (c). (e) Ice thickness h in summer and winter where x > 0.7. (f) Seasonal cycle of ice thickness at the pole hp. (g) Seasonal cycle of the latitude of the sea ice edge θi.

The presence of a substantial seasonal sea ice cover in the simulated climate (Fig. 2c), which is qualitatively consistent with observations, stands in contrast with typical SCM results. It has been a conundrum in SCM studies that they have often struggled to simulate a stable seasonal sea ice cover. The SCM in Thorndike (1992) was forced by a square-wave seasonal cycle, with insolation and longwave radiation switching between a constant summer value for half the year and a constant winter value for the other half. It did not feature any stable seasonal ice. The SCM in Eisenman and Wettlaufer (2009), which had smoothly varying insolation and longwave radiation, did feature stable seasonal ice, but only in a relatively small range of climate forcings. In Eisenman and Wettlaufer (2009, see their Fig. S5), this stable seasonal ice was hypothesized to arise because of a larger seasonal amplitude in the longwave forcing than was used by Thorndike (1992). Moon and Wettlaufer (2012) proposed the contrasting hypothesis that the square-wave seasonal cycle in solar forcing was the essential factor preventing the Thorndike (1992) SCM from simulating a stable seasonal ice cover. They suggested that resolving the seasonal cycle in solar forcing with at least two different values during summer was necessary for stable seasonal ice.

The results presented here suggest an alternative resolution, namely that the inclusion of a horizontal dimension is an essential ingredient for accurately simulating a stable seasonal ice cover. In contrast with Moon and Wettlaufer (2012), who found that square-wave solar forcing did not allow seasonal ice, the model developed here still simulates a large seasonal migration of the ice edge when the cos ωt factor in Eq. (3) is replaced with a square wave (not shown). In contrast with Eisenman and Wettlaufer (2009), who found a narrow SCM forcing range that allowed for stable seasonal ice, seasonal ice occurs in this spatially varied model over a wide range of latitudes.

b. Forced warming (increasing F)

Next, we numerically simulate the equilibrium response of the model to changes in the climate forcing by slowly ramping up F. Beginning with a 200-yr spinup simulation at the initial forcing value (F = −10 W m−2), we increase F in increments of 0.2 W m−2 and integrate the system for 40 years at each value of F. This approach is used in all of the following simulations of the model response to increasing or decreasing climate forcing.

As shown in Fig. 3, the climate steadily warms and the seasonally varying sea ice cover steadily recedes until the pole is ice free throughout the year. The summer ice disappears at F = 2.5 W m−2 and the winter ice at F = 11 W m−2 (see also Fig. 4a).

Fig. 3.

Simulated warming and sea ice loss in the default parameter regime with increasing climate forcing, F. Shown are the annual-mean (a) surface temperature T, (b) surface enthalpy E, and (c) ice thickness h, as functions of location x and forcing F. Black dashed lines show the summer (above) and winter (below) locations of the sea ice edge.

Fig. 3.

Simulated warming and sea ice loss in the default parameter regime with increasing climate forcing, F. Shown are the annual-mean (a) surface temperature T, (b) surface enthalpy E, and (c) ice thickness h, as functions of location x and forcing F. Black dashed lines show the summer (above) and winter (below) locations of the sea ice edge.

Fig. 4.

Simulated sea ice loss and recovery in the default parameter regime with increasing and decreasing climate forcing, F. Shown are the ice area Ai vs (a) the forcing F and (b) the annual-mean hemispheric-mean surface temperature anomaly. The solid thick red and thin blue lines represent the annual-mean ice area in the warming and cooling scenarios, respectively; note that they coincide, implying no hysteresis. The faint solid and dashed lines show the winter and summer ice areas in the warming scenario.

Fig. 4.

Simulated sea ice loss and recovery in the default parameter regime with increasing and decreasing climate forcing, F. Shown are the ice area Ai vs (a) the forcing F and (b) the annual-mean hemispheric-mean surface temperature anomaly. The solid thick red and thin blue lines represent the annual-mean ice area in the warming and cooling scenarios, respectively; note that they coincide, implying no hysteresis. The faint solid and dashed lines show the winter and summer ice areas in the warming scenario.

c. Test for hysteresis

It is noteworthy that the sea ice declines smoothly, with no jumps occurring during the transition from perennial sea ice to seasonally ice-free conditions and then to perennially ice-free conditions. Rather, the summer and winter sea ice edges both respond fairly linearly to F (dashed lines in Fig. 3). This can similarly be considered in terms of the ice area, given by Ai = AH(1 − xi), where AH is the surface area of the hemisphere (Fig. 4a).

Figure 4b shows that the loss of summer and winter ice also relates linearly to the annual-mean hemispheric-mean temperature. The model becomes seasonally ice free at 2°C of warming, and it becomes perennially ice free at 6°C (Fig. 4b). These values compare with comprehensive GCM results from Armour et al. (2011), who find the complete loss of summer and winter ice to occur at 4° and 7°C, respectively.

After the climate has become perennially ice free, we slowly ramp F back down again. We find that the ice recovers during cooling along the same trajectory as the ice retreat during warming, with no hysteresis (Fig. 4). The linearity and reversibility of the response in the present model is consistent with results from most comprehensive GCMs (e.g., Winton 2006; Armour et al. 2011; Winton 2011), and it is in contrast with previous results from EBMs and SCMs.

4. Reduction to EBM and SCM

The two standard types of idealized models discussed above, EBMs and SCMs, both exhibit bistability of the sea ice cover in ostensibly realistic parameter regimes. These results can both be reproduced in the present model. As we show in this section below, the present model reduces to a standard EBM in the limit of no seasonal cycle (S1 = 0), and the present model reduces to standard SCMs at each spatial grid point in the limit of no horizontal heat transport (D = 0).

a. EBM regime (S1 = 0)

When the seasonal amplitude is set to zero in the present model (S1 = 0), the steady-state solution is equivalent to the solution of a standard annual-mean EBM. In this limit, the equilibrium state of the system becomes time independent, with Eq. (2) reducing to

 
formula

Here, the time evolution of the enthalpy E is no longer included, so the surface temperature can be determined directly from this ordinary differential equation in x, without consideration of the vertical heat conduction within the ice.

Equation (10) is equivalent to the formulation of a standard annual-mean EBM [e.g., Eq. (22) in North et al. (1981)]. Analytical solutions can be found by making use of the property that Legendre polynomials are eigenfunctions of the diffusion operator Eq. (6). They are, however, complicated by the dependence of the coalbedo on T. Following Held and Suarez (1974), we approximately solve Eq. (10) for F as a function of xi [e.g., Eqs. (28), (29), and (37) in North et al. (1981)], rather than vice versa. The approximate analytical solution involves a sum of even Legendre polynomials, and we retain terms up to the 40th degree, similar to previous studies (e.g., Mengel et al. 1988).

Figure 5a shows both the analytical solution and numerical model results. It illustrates the evolution of the sea ice edge xi under varied climate forcing F with S1 = 0 and D = D*. The figure illustrates sea ice retreat in the numerical simulation as F is increased (red) until an ice-free state is reached, followed by the onset and advance of the sea ice cover when F is subsequently decreased (blue).

Fig. 5.

(a) Steady-state locations of the simulated ice edge (xi) under varied climate forcing (F) in the parameter regime where the model reduces to an EBM (D = D*, S1 = 0). The numerical solution is indicated for F being slowly ramped up (red) and then back down (blue). Stable and unstable states from the analytic EBM solution are also shown as solid and dashed black lines, respectively. (b) As in (a), but in terms of the enthalpy at the pole (Ep) rather than the ice edge location. The vertical axis on the right indicates corresponding values of the surface temperature at the pole Tp (for E > 0; green) and ice thickness at the pole hp (for E < 0; orange). Both panels also indicate ΔF that is defined as the difference between the bifurcation points Fc (blue square) and Fw (red square).

Fig. 5.

(a) Steady-state locations of the simulated ice edge (xi) under varied climate forcing (F) in the parameter regime where the model reduces to an EBM (D = D*, S1 = 0). The numerical solution is indicated for F being slowly ramped up (red) and then back down (blue). Stable and unstable states from the analytic EBM solution are also shown as solid and dashed black lines, respectively. (b) As in (a), but in terms of the enthalpy at the pole (Ep) rather than the ice edge location. The vertical axis on the right indicates corresponding values of the surface temperature at the pole Tp (for E > 0; green) and ice thickness at the pole hp (for E < 0; orange). Both panels also indicate ΔF that is defined as the difference between the bifurcation points Fc (blue square) and Fw (red square).

The hysteresis loop corresponds to the SICI. Note that a second hysteresis loop corresponding to the snowball earth instability occurs in substantially colder climates, as in standard EBMs. Figure 5a highlights the close agreement between numerical results (blue/red) and the analytical solution (black). It can be shown that the climates indicated in Fig. 5 are stable where the slope is positive (solid black) and unstable where the slope is negative (dashed black) (North 1975a). The system therefore does not support an ice cover with an equilibrium ice edge poleward of xi = 0.98, or 79° latitude.

A measure of the extent of bistability in the system can be obtained by considering the width of the hysteresis loop. We define two critical values of the forcing F: (i) Fw is the value at which the system first transitions to a perennially ice-free pole in a warming scenario, and (ii) Fc is the value at which the wintertime ice cover first reappears in a cooling scenario. These two values of F are indicated in Fig. 5 by a red square and a blue square, respectively. A saddle-node bifurcation occurs at each of these values in the parameter regime shown in Fig. 5. The width of the hysteresis loop is then defined as

 
formula

Note that ΔF may be seen as a societally relevant measure of instability and associated irreversibility, since it indicates how much the radiative forcing would need to be reduced for the sea ice to return after crossing a tipping point during global warming, although it should be noted that this requires long time scales for the climate system to equilibrate.

It is useful to also consider the hysteresis loop in terms of the enthalpy at the pole Ep. Computing Fc and Fw from Ep (as defined in Fig. 5b) gives equivalent results to using xi (as defined in Fig. 5a), since the bifurcations are associated with a transition between an ice-covered pole and an ice-free pole. Here we use both methods interchangeably to compute the bifurcation points, for example using Ep in the numerical results when D = 0.

b. SCM regime (D = 0)

The present model reduces to an array of independent SCMs when meridional heat transport is turned off by setting D = 0. In this case no communication occurs between individual locations and the state of each grid point is determined by the balance of vertical heat fluxes alone. For a given value of x, Eq. (2) reduces to an ordinary differential equation in time representing a standard SCM [e.g., Eq. (2) in Eisenman (2012)]. As in typical SCMs, we find bistability in the parameter regime (D = 0, S1 = ), with ΔF = 7.0 W m−2.

5. Dependence of bistability on parameter values

In the previous section, we showed that the present model reduces to an EBM when S1 = 0 and to an SCM when D = 0, exhibiting bistability during global warming in both regimes. However, the results in section 3b show that there is no such bistability when D = D* and S1 = . This raises the question of how spatial and seasonal variations impact the stability of the system when varying amounts of both are included in the model.

a. Full (D, S1) parameter space

To address this, we performed 441 model simulations that explore how the hysteresis width ΔF varies in the (D, S1) parameter space. We use 21 evenly spaced values for D in the range [0, 0.76] W m−2 K−1 and 21 evenly spaced values for S1 in the range [0, 351] W m−2. In each simulation, F is ramped up and then down, and ΔF is computed. The results are shown in Fig. 6.

Fig. 6.

Dependence of the extent of bistability ΔF on horizontal communication D and seasonal amplitude S1 summarizing the main result of this study. Parameter values associated with typical SCMs and EBMs are indicated, as well as the default parameter values in the present model.

Fig. 6.

Dependence of the extent of bistability ΔF on horizontal communication D and seasonal amplitude S1 summarizing the main result of this study. Parameter values associated with typical SCMs and EBMs are indicated, as well as the default parameter values in the present model.

Figure 6 indicates how the inclusion of both spatial and seasonal variations influences the presence of bistability. Beginning in the EBM regime (D = 0, S1 = ) and adding a diffusivity of just D = 0.1D* suffices to eliminate any hysteresis in the model. Beginning in the SCM regime (D = D*, S1 = 0) and adding a seasonality of just S1 = 0.2 similarly eliminates any hysteresis. Last, if both parameters are reduced by equal factors from their default values (D = D*, S1 = ), bistability does not occur until (D, S1) < 0.3(D*, ). In other words, including even a relatively small amount of both meridional and seasonal variations destroys the bistability that occurs in SCMs and EBMs.

The results in Fig. 6 summarize the main finding of this study. In the following subsections we discuss the underlying mechanisms. We begin by considering the value of ΔF at the origin in Fig. 6, followed by the limiting cases along the horizontal axis (S1 = 0) and vertical axis (D = 0).

b. No seasonal cycle and no heat transport (D = 0, S1 = 0)

In the simplest version of the model there is no seasonal cycle (S1 = 0) and no heat transport (D = 0). In this case, Fc is the forcing associated with the pole being ice free at T = Tm, and Fw is the forcing associated with the pole being ice covered at T = Tm. Hence the width of the hysteresis loop (ΔF) is equal to the jump in solar forcing between ice-covered and ice-free conditions at the pole,

 
formula

as can be seen from Eqs. (2)(4). It is interesting to note in Fig. 6 how the region of bistability connects this point (D = 0, S1 = 0), the standard EBM regime (D = D*, S1 = 0), and the standard SCM regime (D = 0, S1 = ). This may be taken to imply that these three limiting cases of bistability all effectively represent the same physical process, namely multiple stable states due to the ice–albedo feedback.

c. Stabilization from horizontal transport (S1 = 0)

Here we consider how ΔF varies along the horizontal axis of Fig. 6 (i.e., in the EBM regime). Polar temperatures are typically higher for larger meridional heat transport. As D increases, the sea ice cover therefore disappears and also reappears at lower values of F (i.e., Fw and Fc both decrease). This is illustrated in Fig. 7a, where the critical forcings, Fw and Fc, are plotted versus D for both the numerical solution and the analytical solution discussed in section 4a. For the range of forcings F < Fc, the model does not feature SICI bistability, with the only stable state being a finite ice cover. Similarly, the range F > Fw only allows stable states with an ice-free pole. In the range Fc < F < Fw both ice-covered and ice-free conditions are stable, resulting in SICI bistability. These three regimes are illustrated in the inset of Fig. 7a.

Fig. 7.

Dependence of the critical forcing values on D in the EBM regime (S1 = 0). (a) Fc (blue) and Fw (red) vs D for numerical results (squares) and analytical results (solid line). The inset indicates the bistable regime. (b) Hysteresis width ΔF vs D for numerical (squares) and analytical (solid line) approximate solutions. Also indicated are the default value D = D* (dashed vertical line) and the value at which the SICI merges with the snowball earth instability D = Dmax (solid vertical line).

Fig. 7.

Dependence of the critical forcing values on D in the EBM regime (S1 = 0). (a) Fc (blue) and Fw (red) vs D for numerical results (squares) and analytical results (solid line). The inset indicates the bistable regime. (b) Hysteresis width ΔF vs D for numerical (squares) and analytical (solid line) approximate solutions. Also indicated are the default value D = D* (dashed vertical line) and the value at which the SICI merges with the snowball earth instability D = Dmax (solid vertical line).

The dependence of ΔF on D is shown in Fig. 7b. Note that the numerical results approach the analytical results when the spatial resolution increases. We find that ΔF typically decreases with increasing D, with the greatest sensitivity of ΔF occurring when D is small. This indicates that the introduction of diffusive communication between different latitudes has a stabilizing effect on the polar sea ice cover.

At a critical diffusivity Dmax (marked by a solid vertical line in Fig. 7), the small ice cap instability merges with the snowball earth instability. For D > Dmax, there is no stable climate with 0 < xi < 1. This property of EBMs for large values of D has been noted previously (Lindzen and Farrell 1977). For these values of D the global surface temperature becomes relatively isothermal, allowing stable solutions only with T < Tm everywhere or T > Tm everywhere.

The stabilizing effect of increased diffusivity when D < Dmax can be visualized by considering a potential V associated with the model evolution. We identify the model state by the surface temperature at the pole Tp and we define the potential such that dTp/dt = −dV/dTp. Hence valleys in V are stable equilibria and peaks in V are unstable equilibria. The computation of V is described in  appendix C.

We compare two potentials in Fig. 8. One has a smaller diffusivity (D = 0.12 W m−2 K−1) and the other has a larger diffusivity (D = 0.14 W m−2 K−1). We choose values of F such that the unstable equilibria occur at the same value of Tp, using F = 55 and 51 W m−2, respectively. As one might expect intuitively from the effect of diffusion, the higher diffusivity is associated with a smoother potential. In other words, there is a smaller barrier between the two stable equilibria, which suggests that a smaller change in forcing is necessary to transition between the states (i.e., ΔF is reduced).

Fig. 8.

Potential (V) in the EBM regime (S1 = 0) with D = 0.12 W m−2 K−1, F ≡ 55 W m−2 (solid) and with D = 0.14 W m−2 K−1, F ≡ 51 W m−2 (dashed), illustrating how increasing D reduces the barrier between stable equilibria by smoothing the potential wells.

Fig. 8.

Potential (V) in the EBM regime (S1 = 0) with D = 0.12 W m−2 K−1, F ≡ 55 W m−2 (solid) and with D = 0.14 W m−2 K−1, F ≡ 51 W m−2 (dashed), illustrating how increasing D reduces the barrier between stable equilibria by smoothing the potential wells.

Note that explanations involving the diffusive length scale associated with Eq. (10) have also been suggested for the influence of D on the extent of the SICI (Lindzen and Farrell 1977; North 1984). Furthermore, the qualitative relationship we find between D and ΔF is consistent with the conclusion in Eisenman (2012) that that parameter changes in climate models that give rise to thinner ice and warmer ice-free ocean surface temperatures make the system less prone to bistability.

d. Stabilization from seasonal cycle (D = 0)

The influence of the amplitude of the seasonal cycle (S1) in the SCM regime (D = 0) is illustrated in Fig. 9. In this regime, Fc can be readily obtained analytically because it is the point of transition from a perennially ice-free state, in which case the model equations are linear. We find

 
formula

where κ ≡ [1 + (ωcw/B)2]−1/2 (blue line in Fig. 9a).

Fig. 9.

Dependence of the critical forcing values on S1 in the SCM regime (D = 0). (a) Numerical results for Fw (red squares) and Fc (blue squares) vs S1 and the exact analytical solutions for Fc (blue line). The exact analytical solutions for Fw when ice thickness evolution is neglected is also included (red line). (b) Hysteresis width ΔF vs S1 for numerical results (squares), as well as the analytical solution that neglects ice thickness evolution (solid line). The default parameter value S1 = is indicated in both panels by a dashed vertical line.

Fig. 9.

Dependence of the critical forcing values on S1 in the SCM regime (D = 0). (a) Numerical results for Fw (red squares) and Fc (blue squares) vs S1 and the exact analytical solutions for Fc (blue line). The exact analytical solutions for Fw when ice thickness evolution is neglected is also included (red line). (b) Hysteresis width ΔF vs S1 for numerical results (squares), as well as the analytical solution that neglects ice thickness evolution (solid line). The default parameter value S1 = is indicated in both panels by a dashed vertical line.

The analogous relationship between Fw and S1 is nonlinear, since Fw represents the transition from the ice-covered regime, in which the equations are nonlinear. However, we can consider the model as developed in section 2a, neglecting the nonlinear ice thickness evolution added in section 2b. In this case, only the surface albedo changes when T drops below Tm, and Fw (red line in Fig. 9a) takes a form similar to Fc. This leads to a hysteresis width (black line in Fig. 9b) of

 
formula

We therefore find that even in the simple case with no ice thickness evolution, increasing S1 reduces the presence of instability. The derivation of Eqs. (13) and (14) is given in  appendix B.

Including the ice thickness evolution enhances this stabilizing effect. In this case Fw, and hence ΔF, is more sensitive overall to S1 (red squares in Fig. 9a and black squares in Figs. 9b). Note that the diminished slope when S1 > 220 W m−2 in Figs. 9a and 9b is associated with the onset of seasonally ice-free conditions, analogous to Fig. 9d of Eisenman (2012). Since there is no smoothing of the albedo transition in the present model, further bifurcations associated with the transitions to and from seasonally ice-free conditions also occur (not shown), as in Fig. 9h of Eisenman (2012). There the presence of such bifurcations was suggested to be an artifact of a single-column representation. Consistent with this, they cease to occur when D > 0.05 W m−2 K−1 (not shown). However, since Fc and Fw are defined in terms of the presence of ice in winter, seasonally ice-free conditions do not influence the analysis here.

Equation (13) indicates that Fc increases linearly with S1 (Fig. 9a). This can be understood by considering that F = Fc is defined to be the point during a cooling scenario when ice first appears at the pole. This occurs when the winter minimum surface temperature crosses the freezing point (i.e., when ice appears on the coldest day of the year), at which point the ice–albedo feedback causes an abrupt transition to a perennially ice-covered state. In the perennially ice-free regime, the annual-mean value of E is linearly related to F and the seasonal amplitude of E is linearly related to S1 because the governing equations are linear. The larger the seasonal amplitude of E, the warmer the annual-mean value of E can be and still have the winter minimum satisfy E < 0. Hence, larger values of S1 are associated with larger values of Fc. This point is illustrated schematically in Fig. 10 (blue shading).

Fig. 10.

Schematic illustrating the dependence of ΔF on S1 in the SCM regime (D = 0). Turning on the seasonal cycle causes Fc (critical value of F in a cooling scenario; blue) to increase and Fw (critical value of F in a warming scenario; red) to decrease. The solid blue and red lines represent the enthalpy at the pole when S1 = 0 and the colored shading represents the seasonal range at the pole when S1 > 0. Ice is present when E < 0. The decrease in hysteresis width associated with the inclusion of seasonal variations is also indicated (from dark gray ΔF to light gray ΔF).

Fig. 10.

Schematic illustrating the dependence of ΔF on S1 in the SCM regime (D = 0). Turning on the seasonal cycle causes Fc (critical value of F in a cooling scenario; blue) to increase and Fw (critical value of F in a warming scenario; red) to decrease. The solid blue and red lines represent the enthalpy at the pole when S1 = 0 and the colored shading represents the seasonal range at the pole when S1 > 0. Ice is present when E < 0. The decrease in hysteresis width associated with the inclusion of seasonal variations is also indicated (from dark gray ΔF to light gray ΔF).

An analogous argument applies to Fw (red shading in Fig. 10), which indicates the transition point when an initially ice-covered pole is subjected to increasing F. When the seasonal amplitude (S1) increases, the warmest day of summer becomes ice free at an earlier point, that is, at a smaller value of F (red shading in Fig. 10). This argument for Fw applies exactly when ice thickness evolution is neglected. When it is included, the relationship becomes more complicated but the monotonic dependence remains (red squares in Fig. 9a).

Thinking in terms of the potential as in Fig. 8, one can visualize that a larger seasonal amplitude in forcing causes the system to climb higher onto either side of the well around a stable equilibrium. If the seasonal amplitude is sufficiently large, then the system will cross back and forth over the hump that separates the two stable equilibria during the course of each year, thereby having only a single seasonally varying equilibrium solution. A caveat with this explanation, however, is that the potential well associated with the time evolution of the system actually varies seasonally when S1 > 0 and cannot strictly be visualized as remaining constant throughout the year.

Another feature of the climate system that is represented in comprehensive GCMs but not in idealized models is high-frequency weather variability, which could be viewed as noise in the context of the model presented here. Although not included in this analysis (cf. Moon and Wettlaufer 2013), it is likely that the addition of a simple representation of noise would stabilize the model in a manner similar to increasing S1, further diminishing the presence of bistability.

e. Constructing full ΔF(D, S1) space from limiting cases

We can consider the influence of both D > 0 and S1 > 0 together by visualizing the combined influence of D and S1 on the potential, adopting the crude approximation that the potential remains constant throughout the year when S1 > 0, as discussed above. In this picture, larger values of S1 will remove bistability in the presence of a larger hump in the potential well. This indicates that the influence of increasing D (decreasing the height of the hump) and increasing S1 (increasing how large a hump the system can cross seasonally) can be added together to estimate ΔF. Hence we simply add the dependence on D shown in Fig. 7b to the dependence on S1 shown in Fig. 9b to create a conceptually more tractable estimate of the main result of Fig. 6.

In other words, using the functional notation ΔF(D, S1) as a shorthand, we approximate that

 
formula

Negative values of ΔF0 are treated as zero, because they imply that S1 is more than large enough for the system to cross the potential hump seasonally and hence that no bistability occurs. The first term on the right-hand side of Eq. (15) is given by the analytical solution from Eq. (10), the second term is the numerical result because no analytical solution of the nonlinear ice thickness evolution was obtained for Fw with D = 0, and the third term is the analytical expression Eq. (12). The result is shown in Fig. 11.

Fig. 11.

Approximation to ΔF(D, S1) generated by adding the results from Figs. 7b and 9b.

Fig. 11.

Approximation to ΔF(D, S1) generated by adding the results from Figs. 7b and 9b.

It is not surprising that Fig. 11 does not exactly match Fig. 6, giving the approximation that was made by adding the two limiting cases, but the resemblance between the two suggests that the dependence of ΔF on D and S1 can be qualitatively explained by considering the EBM and SCM limiting cases in isolation. Note that this approximation to ΔF (Fig. 11) is smaller than the full numerical solution for ΔF (Fig. 6) at all (D, S1) points, which suggests that the two effects reinforce each other in reducing the bistability.

6. Conclusions

Previous studies using seasonally varying SCMs and spatially varying EBMs have found instabilities in the sea ice cover associated with the ice–albedo feedback. Studies using comprehensive GCMs, however, have typically not found such instabilities. Here we developed a model of climate and sea ice that includes both seasonal and spatial variations. The model accounts for the evolution of a spatially varying surface temperature and sea ice thickness as well as a diffusive representation of meridional heat transport. The governing equation is a partial differential equation in time and latitude.

When the diffusivity was removed (D = 0), we showed that the model reduced to a standard SCM represented as an ordinary differential equation in time with no spatial dimension. When seasonality was removed (S1 = 0), we showed that the model reduced to a standard EBM represented as an ordinary differential equation in latitude with no temporal dimension. When we varied the parameters, we found that including representations of both seasonal and spatial variations causes the stability of the system to substantially increase; that is, any instability and associated bistability was removed (Fig. 6).

Understanding of the underlying physical mechanisms was developed by considering the effects of spatial and temporal variations separately. To this end, we studied the parameter regimes for which the model reduces to a traditional EBM and SCM, respectively. We constructed a potential for the system with two wells separated by a hump, and we showed that horizontal diffusion smooths the potential, reducing the height of the hump. Meanwhile, we argued that the seasonal cycle in the forcing can be approximately visualized as periodically pushing the system up the sides of the well. The larger the seasonal cycle, the higher the system is pushed, and at a critical value of the seasonal amplitude the ball will pass over the hump and seasonally transition between the two wells, thus creating a single seasonally varying stable equilibrium. The combination of smoothed potential wells and seasonal variations eliminates the possibility of any bistability when both parameters are increased to even relatively small values.

This result may help to reconcile the discrepancy between low-order models and comprehensive GCMs in previous studies. Specifically, it suggests that the low-order models overestimate the likelihood of a sea ice “tipping point.”

It is worth emphasizing that the present model contains substantial nonlinearity, including the ice–albedo effect and factors associated with sea ice thickness changes. Such nonlinearity has been shown in SCMs to lead to both accelerating sea ice loss and bifurcations. Nonetheless, the present model simulates sea ice loss that is not only reversible but also has a strikingly linear relationship with the climate forcing as well as with the global-mean temperature. This is in contrast with SCMs and EBMs, and it is consistent with GCMs. The results presented here indicate that the nonlinearities in the model are essentially smoothed out when latitudinal and seasonal variations are included.

These results may have bearing on other processes in the climate system. Previous studies have found other phenomena that feature bifurcations in low-order models but not in most comprehensive GCMs, such as the Atlantic Ocean meridional overturning circulation (e.g., Stommel 1961; Stouffer et al. 2006). The findings presented here raise the possibility that these disparities may be reconciled by similar mechanisms.

Acknowledgments

We are grateful to Kyle Armour and two anonymous reviewers for helpful comments on an earlier version of this article. The authors acknowledge ONR Grant N00014-13-1-0469 and NSF Grant ARC-1107795.

APPENDIX A

Numerical Integration of Model

Here we discuss the numerical approach that we employ to solve the model described by Eqs. (2), (8), and (9). As mentioned above in section 2e, in order to address the issue of diffusive heat transport between different surface temperature regimes in Eq. (9), we create a system of equations that can be visualized as representing two separate layers.

Diffusion takes place in a “ghost” layer with temperature Tg that evolves according to

 
formula

where cg is the heat capacity of the ghost layer, and the first term on the right-hand side causes Tg to relax toward T with time scale τg. All other processes occur in the main layer, whose surface enthalpy evolves as

 
formula

The temperature of the main layer T is as defined in Eq. (9), with Eq. (8) replaced by

 
formula

Note that unlike Eq. (8), which involves ∇2T, Eq. (A3) can be readily solved for T0.

To demonstrate the approximate equivalence between this two-layer system and the model described in section 2a, we begin by defining a diffusive length scale associated with the diffusion operator in Eq. (A1), δxD (cf. Lindzen and Farrell 1977; North 1984), and a time scale associated with changes in T, δtT. If τgcgD−1, then the second term on the right-hand side of Eq. (A1) can be neglected, and furthermore if τgδtT, then the relaxation time scale is sufficiently short that Eq. (A1) can be approximated by Tg = T.

Next, we add Eqs. (A1) and (A2) to get

 
formula

Inserting Tg = T and the definition of E from Eq. (1), we see that the second term on the left-hand side of Eq. (A4) can be neglected as long as cgcw in ice-free locations and cgLfδh/δTg in locations with ice, where δh and δTg are scales of the variations of h and Tg. However, the first constraint on τg imposes a lower limit on cg, as can be seen more clearly when the constraint is rewritten as τg/cgD−1. Hence in the limit of small cg, small τg, and small τg/cg, the two-layer system described in this section reduces to the model comprised of Eqs. (2), (8), and (9).

The two-layer system is considerably more amenable to numerical integration, so we use this for numerical solutions of the model after verifying that our values of τg and cg (see Table 1) are sufficiently small that reducing them does not substantially influence the solution.

We solve the two-layer system by integrating Eq. (A1) using the implicit Euler method, since it is effective for solving the diffusion equation, and integrating Eq. (A2) using the forward Euler method, since it is more straightforward for nonlinear systems. We use the same time step Δt for both equations, which simplifies the coupling.

Numerically, Tg is given as a vector of length n, computed at each point in the domain x = [Δx/2, 3Δx/2, …, 1 − Δx/2], where Δx = 1/n. The diffusive term D2Tg is computed using the finite difference operator for the Laplacian, where is a matrix with nonzero elements

 
formula

and diagonal end terms

 
formula

with and ≡ [0, Δx, …, 1 − Δx] (Bitz and Roe 2001). Note that is a tridiagonal matrix and hence rapidly invertible, which allows for efficient implicit Euler time stepping of Eq. (A1). Note that the inverted matrix used to solve Eq. (A1) for Tg contains not only but also T and therefore needs to be computed at each time step.

We further verified the numerical validity of this two-layer approach by removing the melting ice condition from Eq. (9). This creates an unphysical situation where the ice surface temperature is allowed to be greater than Tm, but it makes the system of Eqs. (2), (8), and (9) more amenable to numerical integration. In this case, the one-layer system of Eqs. (2), (8), and (9) can be directly solved by a single equation that represents an implicit Euler time stepping of T where E > 0 and a forward Euler time stepping of h where E < 0 (cf. Bitz and Roe 2001). We found that numerical integrations of the one-layer and two-layer representations of the system with this adjustment to Eq. (9) were consistent.

APPENDIX B

Analytic Solution for Fc when D = 0

Here we derive the analytic expression Eq. (13), giving the critical forcing at which ice first appears at the pole in a cooling scenario in the limit of no heat transport (D = 0). At the pole (x = 1), the evolution of the system of Eqs. (2), (8), and (9) under ice-free conditions is given by

 
formula

where

 
formula

Note that there is no spatial derivative in this case, since heat transport is set to zero. The solution to this linear system is

 
formula

with κ ≡ [1 + (ωcw/B)2]−1/2 and ϕ ≡ tan−1(ωcw/B) = 2.94 months. A transient term associated with the initial condition is not included in the solution Eq. (B2) since we are interested in the spun-up model behavior. Ice appears when the winter minimum value of T drops below Tm. This minimum occurs at t = ϕ/ω. Substituting this into Eq. (B2), setting T = Tm, and solving for the critical forcing F = Fc, leads to Eq. (13) in section 4. In the model without ice thickness evolution, Fw is obtained in the same way, but with A1 = ai(S0S2) − A + Fb + F and A2 = aiS1. Equation (14) is then given from FwFc.

APPENDIX C

Calculation of Potential Wells

The potential (V) in the EBM regime is calculated from the definition dV/dTp = −dTp/dt. Since we are considering the EBM regime, we consider solutions of the time-varying system discussed in section 2b, in which Tp is proportional to Ep. A complication is that the state of the EBM requires the full temperature field T(x), rather than just TpT(1), to be uniquely determined. Here we take the crude approach of considering only temperature fields T(x) that are the equilibrium solution for some value of the forcing F.

For notational convenience we define the net energy flux at the pole (x = 1) when F = 0 as

 
formula

From Eqs. (1), (2), and (C1), since we are not including ice thickness, we then have

 
formula

The potential depends on the rate at which a transient or initial value of Tp = relaxes toward equilibrium for a given value of F. Defining Fi as the forcing for which Ti is the equilibrium temperature field with at the pole, we can write 0 = Fi + fp(Ti), or

 
formula

When this temperature field is subjected to the forcing F, the polar temperature initially evolves as

 
formula

The potential V(Tp) for a given forcing F is then computed by integrating

 
formula

over a range of initial values of with associated equilibrium forcing values Fi. This gives

 
formula

where Tr is an arbitrary reference temperature and c is a related constant of integration. The analytical EBM solution [Eq. (10)] discussed in the main text is expressed in terms of F as a function of xi using Legendre polynomials to the 40th degree. The field T(x) associated with xi can be computed [Eqs. (25)–(31) in North et al. (1981)]. Here we consider only the value at the pole TpT(1) and the analytical solution is used to express a complicated function for F(Tp). We insert this for Fi = F() in Eq. (C6) and calculate the integral numerically. The resulting potential V(Tp) is shown in Fig. 8 for two different values of D and F.

REFERENCES

REFERENCES
Abbot
,
D. S.
,
M.
Silber
, and
R. T.
Pierrehumbert
,
2011
: Bifurcations leading to summer Arctic sea ice loss. J. Geophys. Res.,116, D19120, doi:.
Armour
,
K. C.
,
I.
Eisenman
,
E.
Blanchard-Wrigglesworth
,
K. E.
McCusker
, and
C. M.
Bitz
,
2011
: The reversibility of sea ice loss in a state-of-the-art climate model. Geophys. Res. Lett.,38, L16705, doi:.
Bendtsen
,
J.
,
2002
:
Climate sensitivity to changes in solar insolation in a simple coupled climate model
.
Climate Dyn.
,
18
,
595
609
, doi:.
Bitz
,
C. M.
, and
G. H.
Roe
,
2001
: Ice and climate modeling: Course notes. [Available online at http://www.atmos.washington.edu/~bitz/.]
Bitz
,
C. M.
, and
G. H.
Roe
,
2004
:
A mechanism for the high rate of sea ice thinning in the Arctic Ocean
.
J. Climate
,
17
,
3623
3632
, doi:.
Björk
,
G.
, and
J.
Söderkvist
,
2002
:
Dependence of the Arctic Ocean ice thickness distribution on the poleward energy flux in the atmosphere
.
J. Geophys. Res.
,
107
,
3173
, doi:.
Björk
,
G.
,
C.
Stranne
, and
K.
Borenäs
,
2013
:
The sensitivity of the Arctic Ocean sea ice thickness and its dependence on the surface albedo parameterization
.
J. Climate
,
26
,
1355
1370
, doi:.
Budyko
,
M. I.
,
1966
: Polar ice and climate (in Russian). Gidrometeoizdat, 32 pp.
Budyko
,
M. I.
,
1969
:
The effect of solar radiation variations on the climate of the Earth
.
Tellus
,
21A
,
611
619
, doi:.
Budyko
,
M. I.
,
1972
:
The future climate
.
Eos, Trans. Amer. Geophys. Union
,
53
,
868
–874, doi:.
Cahalan
,
R. F.
, and
G. R.
North
,
1979
:
A stability theorem for energy-balance climate models
.
J. Atmos. Sci.
,
36
,
1178
1188
, doi:.
Croll
,
J.
,
1875
: Climate and Time in their Geological Relations. Cambridge University Press, 577 pp.
Eisenman
,
I.
,
2007
: Arctic catastrophes in an idealized sea ice model. WHOI Tech. Rep.
2007
-
2
, Woods Hole Oceanographic Institution, 161 pp.
Eisenman
,
I.
,
2010
: Geographic muting of changes in the Arctic sea ice cover. Geophys. Res. Lett.,37, L16501, doi:.
Eisenman
,
I.
,
2012
: Factors controlling the bifurcation structure of sea ice retreat. J. Geophys. Res.,117, D01111, doi:.
Eisenman
,
I.
, and
J. S.
Wettlaufer
,
2009
:
Nonlinear threshold behavior during the loss of Arctic sea ice
.
Proc. Natl. Acad. Sci. USA
,
106
,
28
32
, doi:.
Ferreira
,
D.
,
J.
Marshall
, and
B.
Rose
,
2011
:
Climate determinism revisited: Multiple equilibria in a complex climate model
.
J. Climate
,
24
,
992
1012
, doi:.
Flato
,
G. M.
, and
R. D.
Brown
,
1996
:
Variability and climate sensitivity of landfast Arctic sea ice
.
J. Geophys. Res.
,
101
, 25 767–25 777, doi:.
Gregory
,
J. M.
,
P. A.
Stott
,
D. J.
Cresswell
,
N. A.
Rayner
,
C.
Gordon
, and
D. M. H.
Sexton
,
2002
:
Recent and future changes in Arctic sea ice simulated by the HadCM3 AOGCM
.
Geophys. Res. Lett.
,
29
,
2175
, doi:.
Held
,
I. M.
, and
M. J.
Suarez
,
1974
:
Simple albedo feedback models of the icecaps
.
Tellus
,
26A
,
613
629
, doi:.
Huang
,
J.
, and
K. P.
Bowman
,
1992
:
The small ice cap instability in seasonal energy balance models
.
Climate Dyn.
,
7
,
205
215
, doi:.
Langen
,
P. L.
, and
V. A.
Alexeev
,
2004
: Multiple equilibria and asymmetric climates in the CCM3 coupled to an oceanic mixed layer with thermodynamic sea ice. Geophys. Res. Lett.,31, L04201, doi:.
Li
,
C.
,
D.
Notz
,
S.
Tietsche
, and
J.
Marotzke
,
2013
:
The transient versus the equilibrium response of sea ice to global warming
.
J. Climate
,
26
,
5624
5636
, doi:.
Lian
,
M. S.
, and
R. D.
Cess
,
1977
:
Energy balance climate models: A reappraisal of ice–albedo feedback
.
J. Atmos. Sci.
,
34
,
1058
1062
, doi:.
Lin
,
R. Q.
, and
G. R.
North
,
1990
:
A study of abrupt climate change in a simple nonlinear climate model
.
Climate Dyn.
,
4
, 253–261, doi:.
Lindzen
,
R. S.
, and
B.
Farrell
,
1977
:
Some realistic modifications of simple climate models
.
J. Atmos. Sci.
,
34
,
1487
1501
, doi:.
Maykut
,
G. A.
, and
N.
Untersteiner
,
1971
:
Some results from a time-dependent thermodynamic model of sea ice
.
J. Geophys. Res.
,
76
,
1550
1575
, doi:.
Maykut
,
G. A.
, and
M. G.
McPhee
,
1995
:
Solar heating of the Arctic mixed layer
.
J. Geophys. Res.
,
100
(
C12
), 24 691–24 703, doi:.
Mengel
,
J. G.
,
D. A.
Short
, and
G. R.
North
,
1988
:
Seasonal snowline instability in an energy balance model
.
Climate Dyn.
,
2
,
127
131
, doi:.
Moon
,
W.
, and
J. S.
Wettlaufer
,
2011
:
A low-order theory of Arctic sea ice stability
.
Europhys. Lett.
,
96
, 39001, doi:.
Moon
,
W.
, and
J. S.
Wettlaufer
,
2012
: On the existence of stable seasonally varying Arctic sea ice in simple models. J. Geophys. Res.,117, C07007, doi:.
Moon
,
W.
, and
J. S.
Wettlaufer
,
2013
:
A stochastic perturbation theory for non-autonomous systems
.
J. Math. Phys.
,
54
, 123303, doi:.
Morales Maqueda
,
M.
,
A. J.
Willmott
,
J. L.
Bamber
, and
M. S.
Darby
,
1998
:
An investigation of the small ice cap instability in the Southern Hemisphere with a coupled atmosphere–sea ice–ocean–terrestrial ice model
.
Climate Dyn.
,
14
,
329
352
, doi:.
Müller-Stoffels
,
M.
, and
R.
Wackerbauer
,
2011
:
Regular network model for the sea ice-albedo feedback in the Arctic
. Chaos,
21
, 013111, doi:.
Müller-Stoffels
,
M.
, and
R.
Wackerbauer
,
2012
:
Albedo parametrization and reversibility of sea ice decay
.
Nonlinear Processes Geophys.
,
19
,
81
94
, doi:.
North
,
G. R.
,
1975a
:
Analytical solution to a simple climate model with diffusive heat transport
.
J. Atmos. Sci.
,
32
,
1301
1307
, doi:.
North
,
G. R.
,
1975b
:
Theory of energy-balance climate models
.
J. Atmos. Sci.
,
32
,
2033
2043
, doi:.
North
,
G. R.
,
1984
:
The small ice cap instability in diffusive climate models
.
J. Atmos. Sci.
,
41
,
3390
3395
, doi:.
North
,
G. R.
, and
J. A.
Coakley
,
1979
:
Differences between seasonal and mean annual energy balance model calculations of climate and climate sensitivity
.
J. Atmos. Sci.
,
36
,
1189
1204
, doi:.
North
,
G. R.
,
R. F.
Cahalan
, and
J. A.
Coakley
Jr.
,
1981
:
Energy balance climate models
.
Rev. Geophys.
,
19
,
91
–121, doi:.
North
,
G. R.
,
J. G.
Mengel
, and
D. A.
Short
,
1983
:
Simple energy balance model resolving the seasons and the continents: Application to the astronomical theory of the ice ages
.
J. Geophys. Res.
,
88
(
C11
),
6576
6586
, doi:.
Ridley
,
J.
,
J.
Lowe
, and
D.
Simonin
,
2008
:
The demise of Arctic sea ice during stabilisation at high greenhouse gas concentrations
.
Climate Dyn.
,
30
,
333
341
, doi:.
Ridley
,
J.
,
J.
Lowe
, and
H. T.
Hewitt
,
2012
:
How reversible is sea ice loss?
Cryosphere
,
6
,
193
198
, doi:.
Rose
,
B. E. J.
, and
J.
Marshall
,
2009
:
Ocean heat transport, sea ice, and multiple climate states: Insights from energy balance models
.
J. Atmos. Sci.
,
66
,
2828
2843
, doi:.
Sellers
,
W. D.
,
1969
:
A global climatic model based on the energy balance of the Earth–atmosphere system
.
J. Appl. Meteor.
,
8
,
392
400
, doi:.
Stommel
,
H.
,
1961
:
Thermohaline convection with two stable regimes of flow
.
Tellus
,
13B
,
224
230
, doi:.
Stouffer
,
R. J.
, and Coauthors
,
2006
:
Investigating the causes of the response of the thermohaline circulation to past and future climate changes
.
J. Climate
,
19
,
1365
1387
, doi:.
Suarez
,
M. J.
, and
I. M.
Held
,
1979
:
Sensitivity of an energy balance climate model to variations in the orbital parameters
.
J. Geophys. Res.
,
84
(C8),
4825
4836
, doi:.
Thorndike
,
A. S.
,
1992
:
A toy model linking atmospheric thermal radiation and sea ice growth
.
J. Geophys. Res.
,
97
(
C6
),
9401
–9410, doi:.
Tietsche
,
S.
,
D.
Notz
,
J. H.
Jungclaus
, and
J.
Marotzke
,
2011
: Recovery mechanisms of Arctic summer sea ice. Geophys. Res. Lett.,38, L02707, doi:.
Winton
,
M.
,
2006
: Does the Arctic sea ice have a tipping point? Geophys. Res. Lett.,33, L23504, doi:.
Winton
,
M.
,
2008
: Sea ice–albedo feedback and nonlinear Arctic climate change. Arctic Sea Ice Decline: Observations, Projections, Mechanisms, and Implications, Geophys. Monogr., Vol. 180, Amer. Geophys. Union, 111–131.
Winton
,
M.
,
2011
:
Do climate models underestimate the sensitivity of Northern Hemisphere sea ice cover?
J. Climate
,
24
,
3924
3934
, doi:.