On the nature of the sea ice albedo feedback in simple models.

We examine the nature of the ice-albedo feedback in a long-standing approach used in the dynamic-thermodynamic modeling of sea ice. The central issue examined is how the evolution of the ice area is treated when modeling a partial ice cover using a two-category-thickness scheme; thin sea ice and open water in one category and "thick" sea ice in the second. The problem with the scheme is that the area evolution is handled in a manner that violates the basic rules of calculus, which leads to a neglected area evolution term that is equivalent to neglecting a leading-order latent heat flux. We demonstrate the consequences by constructing energy balance models with a fractional ice cover and studying them under the influence of increased radiative forcing. It is shown that the neglected flux is particularly important in a decaying ice cover approaching the transitions to seasonal or ice-free conditions. Clearly, a mishandling of the evolution of the ice area has leading-order effects on the ice-albedo feedback. Accordingly, it may be of considerable importance to reexamine the relevant climate model schemes and to begin the process of converting them to fully resolve the sea ice thickness distribution in a manner such as remapping, which does not in principle suffer from the pathology we describe.

1. Introduction

The original thermodynamic theory coupling sea ice and climate dealt with the system as a column of atmosphere, ice, and ocean [Maykut and Untersteiner, 1971]. This approach is the cornerstone of contemporary theoretical studies [e.g., Eisenman, 2012, and references therein] and it underlies the thermodynamics of the sea ice components of all contemporary climate models [e.g., AOMIP, 2013]. We understand that the ice cover presents a distribution of ice thicknesses g(h) to the atmosphere and ocean that force its growth, decay, and deformation [Thorndike et al., 1975]. Although the treatment of this distribution as a continuous differentiable function is based on clear reasoning, its practical implementation in either simple or complex models is a major challenge. Hibler [1979] (H79 throughout this paper) developed an implementation scheme for the theory of Thorndike et al. [1975] wherein both ice and open water are considered as part of a grid cell. In such a so-called two-category-thickness scheme, one category consists of thin sea ice and open water and the other category is “thick” sea ice. The areal fraction of both categories is computed at each time step. The scheme emerged at a time when the perennial ice state was not questioned. However, it does not conserve the total area of ice in a grid cell and hence, due to the nature of the ice-albedo feedback, is of particular importance as the state of the ice cover changes from perennial to seasonal. Here, we demonstrate this in a simple model. The question of how, and how rapidly, the ice cover may decay toward the seasonal state is the main implication of the analysis that follows.

It is important to note that a substantial literature on the simulation of fully resolved sea ice thickness distributions g(h) began about 20 years ago [e.g., Flato and Hibler, 1995], an important approach being the application of the Dukowicz and Baumgardner [2000] remapping scheme to g(h) by Lipscomb [2001]. Such approaches are not in principle influenced by the particular problem we discuss that is associated with a two-category-thickness scheme. However, Earth System Models that continue to use a two-category-thickness scheme [e.g., Watanabe et al., 2011; Notz et al., 2013], or any area-thickness scheme (e.g., 15 of 29 models in the Arctic Ocean Model Intercomparison Project) [AOMIP, 2013] would suffer from the pathologies we examine here.

1.1. Multiple Sea Ice Cover States

A main focus of the attempt to discern the origin of the decline of the Arctic sea ice cover is the evolution of the summer sea ice minimum [e.g., Kwok and Untersteiner, 2011] and the associated question of whether future summers will be ice free, so that there is ice only in winter. The approaches to the problem range from theoretical treatments [e.g., Thorndike, 1992; Eisenman and Wettlaufer, 2009; Mueller-Stoffels and Wackerbauer, 2011; Abbot et al., 2011; Moon and Wettlaufer, 2012; Eisenman, 2012; Stranne and Björk, 2012; Moon and Wettlaufer, 2013] and global climate model simulations [e.g., Holland et al., 2006; Winton, 2008; Tietsche et al., 2011] to interpretation of observations [e.g., Serreze, 2011; Stroeve et al., 2012; Agarwal et al., 2012].

The rudiments of the ice-albedo feedback provide the framework for examining the nature of transitions from the perennial ice state to either a seasonal or ice-free sea ice state. In the framework of simplified versions of the column model of Maykut and Untersteiner [1971], the ice-albedo feedback treats the sea ice albedo as a function of ice thickness h, transitioning continuously from that of sea ice to that of the ocean [Eisenman and Wettlaufer, 2009; Eisenman, 2012; Moon and Wettlaufer, 2012]. As the greenhouse gas forcing (modeled as an additional surface heat flux ) increases, these theoretical approaches capture the nature and general conditions of the transitions between perennial, seasonal, and ice-free states. Eisenman [2012] provides a recent summary of the models and methods used to predict four general scenarios under which ice retreat may occur as increases.

2. Partial Ice Cover and the Ice-Albedo Feedback

2.1. Column Models

Due to the strength of the ice-albedo feedback even in the simplest of models, it is important to attempt to model partial ice cover, which requires an ocean mixed layer that is in communication with the atmosphere unless the ocean is completely ice covered. The most common two-category-thickness methodology for ice area A evolution appears to have originated from Hibler [1979]. It is prudent to take the minimalist approach to demonstrate the key matters at hand, which we do by implementing the H79 approach in the simple column model of Eisenman and Wettlaufer [2009], which is derived from that of Maykut and Untersteiner [1971]. We first summarize the relevant aspects of Eisenman and Wettlaufer [2009] and Eisenman [2007], and then in section 2.2, we describe the implementation of H79.

When the temperature of the ice , it evolves along with the ice thickness h according to and where L is the latent heat of fusion per unit volume, c is the specific heat capacity of ice at constant pressure, and k is the thermal conductivity. The sum of sensible, latent, downward and upward longwave, and shortwave heat fluxes at the surface is given as F, and all but the upward longwave flux are specified from observed radiation climatology as in Eisenman and Wettlaufer [2009]. The flux from the ocean mixed layer into the base of the ice is F.

When , the temperature evolves along with the ice thickness according to

The ocean mixed layer is treated as a thermodynamic reservoir with a typical observationally based characteristic depth of H = 50 m and heat flux entrained through the bottom of the mixed layer of F = 0.5 W m−2. The turbulent heat flux between the ocean and the underside of the ice is a complex quantity modeled crudely here as being proportional to the elevation of the mixed layer temperature T above freezing by , in which ρ and c are the density and specific heat capacity at constant pressure of seawater, is the heat transfer coefficient, and is the square root of kinematic stress at ice-ocean interface, also known as the friction velocity [e.g., Maykut and McPhee, 1995]. A typical observational value of leads to with . Measurements show that while in the upper summer mixed layer T can be as much as above freezing, it is an order of magnitude smaller in winter, giving a seasonally averaged F of about 5 W m−2 [Maykut and McPhee, 1995]. Therefore, according to whether , a continuously evolving seasonal cycle is captured, using (1) and (2) or (3) and (4).

2.2. Modeling Partial Sea Ice Cover

Now we proceed to the issue of modeling partial ice cover using column models. We begin by summarizing the commonly used approach to this problem developed by Hibler [1979]. Such a method can be rationalized physically for the perennial ice cover for which it was developed, but we discuss its behavior when the ice fraction decreases, such as is relevant during the transition to seasonal ice as the observed state of the ice cover changes [Perovich and Polashenski, 2012].

The H79 methodology to determine ice concentration A, or the fraction of a grid cell covered by ice, is the core focus here. This requires a form of homogenization over the subgrid scale to account for open water. As noted in section 1, this method provides the framework for parameterizations in two-category-thickness distributions for area-thickness modeling schemes. Of the 29 models participating in the Arctic Ocean Model Intercomparison Project [AOMIP, 2013], 15 use area-thickness schemes, as do Earth System Models such as MPI-ESM-LR [Notz et al., 2013] and MIROC [Watanabe et al., 2011]. For clarity, but without loss of generality, we discuss the H79 approach in terms of a model that includes a single grid cell. The approach applies to either the area of ice in the grid cell or, as is done in H79, the fraction of the grid cell covered by ice A. Although in many of the equations that follow these can be used interchangeably, we use areal fraction for consistency with H79. This means that in the ice-covered fraction of the grid cell ice thickness h becomes the volume per area V, which has units of length. Variables and constants are defined in Table 1.

Table 1

Descriptions and Values of State Variables and Model Parameters

Symbol	Description	Units/Value
h	Ice thickness	m
V	Ice volume per unit grid cell area	m
A	Ice areal fraction in grid cell
Ti	Ice surface temperature
Tml	Ocean mixed layer temperature
h0	Equivalent thickness for newly formed ice	0.5 m
L	Latent heat of fusion per unit volume	3 × 108 J/m3
k	Thermal conductivity of ice	2 W/m/K
ρw	Density of water at constant pressure	1 × 103 kg/m3
cpi	Specific heat capacity of ice at constant pressure	2 × 106 J/m3/K
cpw	Specific heat capacity of water at constant pressure	4 × 106 J/m3/K
αi	Albedo of ice	0.65
αml	Albedo of ocean mixed layer	0.20
Fsw	Shortwave radiation at ice or ocean surface	Seasonal; W m−2
Flw	Longwave radiation at ice or ocean surface	Seasonal; W m−2
Ftop	Net surface sensible, latent, and radiative heat flux	W m−2
	Greenhouse gas forcing	0–30 W m−2
ch	Ocean-ice heat transfer coefficient	0.006
	Ocean-ice friction velocity	0.5 cm s−1
γ	Ocean-ice heat exchange coefficient =	120 W m−2/K
Hml	Mixed layer depth	50 m
Fentr	Heat flux entrained into mixed layer base	0.5 W m−2
FB	Seasonal average ice-ocean heat flux	5 W m−2
Fml	Total heat flux into the mixed layer (equation (10))	Seasonal; W m−2

The ice concentration increases when T reaches zero and continues to cool so that the mixed layer flux imbalance F drives the creation of new ice as

An “equivalent thickness” h0 is assigned to the new area ascribing volume to it. Thus, area increases only when the mixed layer freezes, but once it does so, the new volume of that ice increases only by increasing the ice thickness at fixed area. Because, within the framework of column models, sea ice growth rate is calculated (or specified) as a function of ice thickness and season, the value of h0 controls the rate at which the ice cover grows. The value of h0 used in H79 is 50 cm. Although the growth rate in winter decreases by a factor of four as open water solidifies to a thickness of 50 cm, the ice concentration in H79 increases based on the growth rate for open water ∼12 cm d−1 (see Figure 3 in H79). Importantly, in this and similar two-category models, the open water fraction is not meant to represent an entirely ice-free region. Rather, the model domain is split into a fraction containing thick ice, with the rest covered by a mixture of open water and thin ice, such as in leads. The volume of this thin ice is assumed to be negligible compared to the thick ice volume, which as we shall see is one of the problems in dealing quantitatively with processes such as ice-albedo feedback.

Energy balance dictates that area decays in this model when volume ablates and hence

The proportionality between volume and area rates of change is based on an argument about the ice thickness distribution in the model domain under the following assumptions [Hibler, 1979]: (a) the ice is linearly distributed in thickness between 0 and , thereby giving a mean thickness of V/A, and (b) all of this ice melts at the same rate. As illustrated schematically in Figure 1, this gives a rate of area decay as the rate of thickness decay times the inverse slope of the thickness distribution;

Figure 1

Schematic illustrating the proportionality between the rate of change of ice area A and the thermodynamic decrease of volume following Hibler [1979].

We note here that ice growth is a nonlinear function of thickness and here it is computed under the assumption that all ice within A is of the mean thickness V/A as opposed to the linear distribution between 0 and used for ablation.

Finally, the persistent convergence and divergence of the wind field results in an observed net average annual export of v0 = 10% of the ice area. Thus, the ice dynamics are represented in such a model by requiring that , and a term is added to the area evolution equation (which accounts for volume export). Export is included in the results shown in Figure 2, but to avoid the clutter in the theoretical development, we omit the term in the equations that follow because it has no effect on the main points.

Figure 2

The three plots show two sets of solutions for the evolution of the dimensionless ice volume, ice area , and the mixed layer energy as a function of greenhouse gas forcing (W m−2). The dash-dot lines (open diamonds) show the solutions of the equations above in which area evolution is complete (area evolution is incomplete). In both cases, red shows the end of the summer (late August) and blue the end of the winter (late March). The differences are discussed in the text in more detail, but we note here that a principal feature is that when area evolution is incomplete, substantially higher values of are required before seasonal or perennial ice is lost.

Using such a scheme, one can derive a partial ice cover model from the column treatment of section 2.1 as follows. We determine ice volume rather than ice thickness. In the ice-covered fraction of the model domain A, the vertical thermodynamic growth of the ice is represented by rewriting (2) and (4) as and

The total heat flux into the mixed layer is written as where is the net surface longwave radiation balance, which depends on the mixed layer temperature T, the shortwave radiative balance is F, and the albedo of the mixed layer is α [Eisenman and Wettlaufer, 2009; Eisenman, 2007]. Therefore, if T > 0, this leads to heating or cooling according to and no new ice area is formed, F = 0. However, when the mixed layer reaches the freezing temperature (T = 0), supercooling is prohibited such that , and any additional heat loss is available to form new ice ( ).

2.3. Area Evolution and Ice-Albedo Feedback

The principal issue here is that equations (8) and (9) do not correctly capture the evolution of the areal fraction of ice. The nub of the matter is the appropriate grid homogenization of the mixture theory. We write the volume of ice in a grid cell of area A as and use the logic of H79, that there are two ice categories, thick (h) and thin (h0), the latter intended to represent both open water and ice up to a value of h0. Then, one would naturally write the average ice thickness as and hence the growth rate in a grid cell would be , where the dots denote differentiation with respect to time. However, the scheme described in section 2.2 implements the assumption that the mass of thin ice is small relative to that of thick ice in the following manner. On the one hand, this assumption imposes in order to compute the average “thick ice” thickness as . On the other hand, the reality that thin ice grows rapidly is implemented by computing the growth rate as with denoting the growth rate of ice with thickness h/A and the term set to zero (see equation (15) of H79). This is the main point. Although such an approach can be physically rationalized in the central perennial pack ice cover when the mass of thin ice is small relative to that of thick ice, it becomes questionable as the area of thin ice increases in magnitude and the rate of change of its area becomes significant. Indeed, as the ice cover transitions to a seasonal state undergoing rapid and large changes in areal coverage, and during a part of such a seasonal cycle when A becomes vanishingly small, the concept of an average thickness defined as , exhibits obvious pathologies.

The appropriate conservation law requires the addition of the term to the left-hand side of equations (8) and (9). To assess the importance of such a term and facilitate simple analysis and interpretation of general features of the freezing and melting process, it is prudent to render these equations dimensionless through the introduction of the following scalings; where is the threshold thickness mentioned above, c is again the specific heat capacity of the ice, is the temperature difference over ice of thickness h0 and hence F is the associated conductive heat flux. The dimensionless ratio is a Stefan number, which represents the relative importance of latent heat to the specific heat in the ice, and is large ( ) here. These scalings lead to dimensionless versions of equations (8) and (9) appropriately modified to include the area evolution as and where the fluxes are just the dimensional fluxes scaled by F. Because , the balance in both equations (13) and (14) requires that showing that the neglected term is of the same order of magnitude as that kept in the scheme. Importantly, in transitioning to a seasonal ice state driven by the ice-albedo feedback, correctly capturing the rapidly changing evolution of ice area is crucial, otherwise it is found that ice loss is only expected in untoward parameter regimes. Indeed, when Eisenman [2007] neglected this term and used the correct value of the latent heat of fusion of ice he could not simulate a realistic ice cover. For example, he found that in order to obtain multiple sea ice states under greenhouse gas forcing of 1.5 times the present value, he had to artificially decrease L by factors ranging from 4 to 10. This reflects the fact that, for a given radiation balance, the artificially large latent heat flux associated with the neglect of the extra term in equations (13) and (14) can only be balanced by positing an artificially low value of L. Moreover, when using the correct thermophysical constants, he had to use 5 times (16 times) the present greenhouse gas concentration to transition from perennial to seasonal ice (multiple sea ice states).

2.4. Energy Flux Conservation

Having demonstrated the size of the missing term, we return to dimensional variables in this section. During the melt season, the contribution to the volume evolution of the lateral melting term is insufficient to conserve energy flux. In the H79 scheme, this lateral melting is calculated indirectly in order to maintain the functional form of the model sea ice thickness distribution, save for the constant difference owing to the change of the mean. Because greater than 90% of the incident sunlight is absorbed by open water, we understand that the partitioning of the ablation of the ice cover between top, bottom, and lateral boundaries is, among other factors, a complex function of the open water fraction and ice floe perimeter [e.g., Maykut and Perovich, 1987; Steele, 1992]. The impasse faced by the H79 approach is discussed by Steele [1992] in terms of the lack of an explicit equation for ice floe perimeter. It is thus natural to ask for the origin of the heat source for lateral melting within this model framework.

We examine the partitioning of vertical and lateral oceanic-heat fluxes to account for the contribution . When the average thickness is used to determine the heat flux required to balance the volume change originating in vertical ablation , the result differs from the analogous procedure in which the ice thickness is distributed evenly from 0 to . Thus, we conclude that part of this heat flux difference over an ice area A is used for lateral melting. Hence, if is the net heat flux over sea ice of thickness h, we write as

For simplicity of exposition, we view the principal contribution to as arising from shortwave radiative fluxes as

Therefore, the lateral oceanic-heat flux ablating the ice is . Finally, in order to conserve energy flux balance at each time step during the melt season, we subtract the lateral oceanic-heat flux from the evolution of ocean sensible heat, thereby avoiding an anomalous increase in ocean heat content.

3. Discussion

Now that the essential point has been made using this simple analysis, in Figure 2, we show the dramatic effect of employing a scheme that deals with the area evolution as discussed above. The missing-area (and hence ice-mass) term, as described in the argument leading to equation (15), has the basic effect of neglecting a leading-order latent heat flux in the energy balance; because and have the same sign, not including the missing area leads to a larger latent heat flux. Under the same radiative forcing, the consequences of this missing latent heat flux can be clearly laid bare. First, because the effective latent heat flux is larger than it should be, the volume of ice in steady state with the same radiative balance can be larger as seen in the top plot. The most distinct case is for , where the maximum sea ice thickness is 3.5 m (2 m) when this latent heat flux is ignored (included). Second, when this latent heat flux (and hence the associated ice areal change) is ignored, a larger value of ΔF0 (∼5 W m−2) is required before both the perennial and seasonal ice states vanish. Third, the range of over which seasonal sea ice exists in a stable state is infinitesimal (practically nonexistent) when the area evolution is incomplete whereas it is ∼5 W m−2 when the area evolution and hence latent heat flux is complete. Finally, since the ice cover vanishes at smaller values of greenhouse forcing when area evolution is treated completely, and hence the ice-albedo feedback is appropriately captured, the heat content of the exposed mixed layer is larger. We note that in the original H79 treatment there was no mixed layer; the ocean temperature was constrained to lie on the freezing point and any excess heat absorbed was immediately added to a basal heat flux and applied to the underside of the thick ice. This treatment is clearly unrealistic in the limit of a vanishing ice cover when the missing term in such a scenario becomes all the more important, and a direct application of this scheme simply amplifies the differences shown in Figure 2, so we do not include these figures here.

It is important to note that although we have focused on the simplest (two-category-thickness) schemes, our arguments can be generalized. In a multithickness model, and . Hence, regardless of the scheme, one must provide some form of rule for , which may for example include the lateral heat flux to sea ice. Because one must construct some rule to evolve the open water or thin ice fraction, our main point is general, whereas using a remapping scheme [Dukowicz and Baumgardner, 2000] to solve a thickness distribution equation [Lipscomb, 2001] is in principle free from this problem.

4. Conclusion

We describe how a long standing approach used in the thermodynamic modeling of sea ice does not treat the complete evolution of the ice area and thus cannot capture the influence of the ice-albedo feedback. The missing-area term, as described in the argument leading to equation (15), has the effect of neglecting a leading-order latent heat flux in the energy balance. By deriving energy balance models for partial ice cover with and without the appropriate area evolution, we have demonstrated the sensitivity of the results to the missing area. It is found to be particularly important in a decaying ice cover approaching seasonally ice-free conditions. Although we have not independently analyzed how this erroneous treatment of area evolution has propagated through the range of GCMs used, our analysis indicates the possibility that it could in fact be one of the underlying features responsible for the observed recent Arctic sea ice decline being more rapid than is forecast using a multimodel ensemble mean [e.g., Stroeve et al., 2012]. Thus, it is suggested that it may be of considerable importance to reexamine the relevant climate model schemes and to begin the process of converting them to fully resolve the sea ice thickness distribution in a manner such as remapping [Dukowicz and Baumgardner, 2000; Lipscomb, 2001], which does not in principle suffer from the pathology we have described here.