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.
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).
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)
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.
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:
which includes fluxes from top-of-atmosphere net solar radiation aS, outgoing longwave radiation (OLR) L, meridional heat transport in the atmosphere D∇2T, 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,
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
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.
Following Budyko (1969), we represent the OLR as a linear function of the surface temperature,
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)
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
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 = cw∂T/∂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 T ≤ Tm at the top surface to Tm at the bottom surface. The heat flux upward through the ice is then equal to k(Tm − T)/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
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:
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.
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).
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).
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
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).
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
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.
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.
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,
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.
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).
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
where κ ≡ [1 + (ωcw/B)2]−1/2 (blue line in Fig. 9a).
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
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).
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
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.
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.
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.
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.
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
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
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 τg ≪ cgD−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.
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 cg ≪ cw in ice-free locations and cg ≪ Lfδ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/cg ≪ D−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 D∇2Tg is computed using the finite difference operator for the Laplacian, where is a matrix with nonzero elements
and diagonal end terms
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.
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
Note that there is no spatial derivative in this case, since heat transport is set to zero. The solution to this linear system is
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(S0 − S2) − A + Fb + F and A2 = aiS1. Equation (14) is then given from Fw − Fc.
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 Tp ≡ T(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
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
When this temperature field is subjected to the forcing F, the polar temperature initially evolves as
The potential V(Tp) for a given forcing F is then computed by integrating
over a range of initial values of with associated equilibrium forcing values Fi. This gives
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 Tp ≡ T(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.