Sea-ice models for climate study: retrospective and new directions

Abstract

Scientific observations of sea ice began more than a century ago, but detailed sea-ice models appeared only in the latter half of the last century. The high albedo of sea ice is critical for the Earth’s heat balance, and ice motion across the ocean’s surface transports fresh water and salt. The basic components in a complete sea-ice model must include vertical thermodynamics and horizontal dynamics, including a constitutive relation for the ice, advection and deformational processes. This overview surveys topics in sea-ice modeling from the global climate modeling perspective, emphasizing work that significantly advanced the state of the art and highlighting promising new developments.

1. Introduction

Formed from frozen sea water rich in biological and chemical species, sea ice exists as a thin layer at the interface of the ocean and atmosphere, quite sensitive to small changes in temperature and radiative forcing. Salts and other inclusions within sea ice alter its thermal and morphological properties. Sea ice can be studied at many scales, from the microscale of individual crystals to the basin scale (thousands of kilometers). Individual ice floes range in extent from ˜1 m to tens of kilometers. Global climate models (GCMs) are mainly concerned with the basin scale, but modelers can and do include horizontal scales down to 10 km or less. Smaller-scale processes are parameterized; that is, they are modeled in terms of larger-scale variables.

Sea ice grows in up to 10% of the world ocean during some part of the year. In the Arctic, the ice reaches a maximum extent of ˜15 × 10⁶km² in late winter, covering the entire Arctic Ocean and many peripheral seas. The minimum extent occurs in September, when ice is confined to the central Arctic. The minimum extent averaged 7 × 10⁶km² during 1979−2008 but has been steadily decreasing. The past four years, 2007-10, have had the four lowest September minima during the satellite era (since 1979), with the all-time low of 4.3 × 10⁶km² occurring in 2007 (US National Snow and Ice Data Center, http://nsidc.org/arcticseaicenews/)).

In the Antarctic, total sea-ice extent has not changed significantly in recent decades. Most ice that forms in the winter melts in the summer; the total extent ranges from about 18 × 10⁶km² at the end of winter to 3 × 10⁶km² at the end of summer. Much of the Antarctic has been sheltered from climate warming as a result of increased winds (which isolate the Antarctic from warming regions to the north) and strong vertical ocean mixing (which efficiently removes heat from surface waters). An important difference between the two hemispheres is geographically derived: the Arctic Ocean provides a maritime polar climate, while Southern Hemisphere sea ice is greatly influenced by the massive Antarctic ice sheet and its polar continental climate. Sea ice has retreated, however, around the Antarctic Peninsula, the one part of Antarctica that is rapidly warming.

Sea ice in both hemispheres is sensitive to climate warming, in part because of ice-albedo feedback. When highly reflective ice is removed, solar radiation absorbed in the darker ocean can melt the remaining ice from beneath, accelerating the retreat. The feedback works in the other direction as well: cooling leads to more extensive ice cover and more cooling.

Thus, sea ice is important for climate because it insulates the cold atmosphere from the warmer ocean in winter, and it reflects much of the incoming solar radiation in spring and summer. The ice also alters the density structure of the ocean, increasing the salinity of the upper ocean as it grows and freshening the ocean as it melts. As sea ice drifts on the surface of the ocean, it carries pollutants and dust that can affect the ocean ecosystem. Also, sea ice contains a great deal of non-ice material, mostly taken up from the ocean as the ice freezes, such as salt and algae. These inclusions not only influence the ocean, but often affect the sea ice itself by modifying radiative transport, heat conduction, and the solid/brine matrix within the ice.

Physical processes acting on sea ice are commonly divided into two categories: thermodynamic processes, which involve the transfer of heat or radiation, and dynamic processes, which move and deform the ice. Modern sea-ice models treat both kinds of processes. The roots of modern sea-ice modeling go back to the International Geophysical Year (IGY), 1957−58. Based on IGY observations, Untersteiner (1964) developed a detailed model of heat conduction and surface growth and melting of sea ice. A few years later, Maykut and Untersteiner (1971) created the first sophisticated numerical model of sea-ice thermodynamics, providing a foundation for today’s models. During the 1970s, field observations from the Arctic Ice Dynamics Joint Experiment (AIDJEX) spawned major advances in the understanding and modeling of sea-ice dynamics (e.g. Parmerter and Coon, 1972; Coon and others, 1974; Hibler, 1979; Hopkins and others, 1991). Several seminal papers from this period form the basis for current model components, as discussed below.

During the 1980s, sea-ice models began to be incorporated into large-scale climate models. Hibler and Bryan (1987) coupled a sea-ice model to an ocean model and opened up many new considerations for sea-ice models. The first sea-ice models in GCMs were necessarily simple because of limited computing power. Sea ice was typically modeled as a motionless or freely drifting material with zero heat capacity and uniform thickness. During the past two decades, however, the models have become steadily more complex. The authors of this review have been closely involved in the development of the Los Alamos Sea Ice Model, CICE, which is the sea-ice component of the Community Earth System Model (CESM) based at the US National Center for Atmospheric Research. Other well- known models include the University of Victoria (Bitz and others, 2001), Thickness and Enthalpy Distribution (TED; Zhang and Rothrock, 2001) and Louvain-la-Neuve (LIM; Vancoppenolle and others, 2009a) ice models.

Sea-ice models continue to evolve, as smaller-scale and high-order processes are added. The following sections describe the current state of the art in sea-ice modeling for GCMs, noting areas that need to be improved and current efforts to do so.

2. Ice Thickness Distribution

The Arctic and Antarctic sea-ice packs are mixtures of open water, thin first-year ice, thicker multi-year ice, and thick pressure ridges. Sea ice varies in thickness from a few centimeters in newly opened cracks, or leads, to 10−20 m in ridges. (In comparison, icebergs calved from glacial ice can be tens to hundreds of meters thick.) First-year sea ice, which dominates the ice pack in the peripheral Arctic seas and most of the Antarctic, grows up to 2 m thick. Perennial ice – that which survives at least one melt season – typically grows to a thickness of 2−3 m, but the extent and thickness of perennial ice have decreased significantly during the past few decades (Kwok and Rothrock, 2009), due in part to increased export (Rigor and Wallace, 2004; Hunke and Bitz, 2009).

An essential aspect of sea-ice thermodynamics is the variation of growth and melting rates for different ice thicknesses. Because conduction is proportional to the vertical derivative of the temperature, thin ice grows and melts more quickly than does thicker ice. Similarly, thinner ice is more likely to undergo mechanical deformation than thicker ice. Most early sea-ice models neglected thickness variations, but modeling studies have confirmed that both thermodynamic and dynamic properties of the ice pack depend on how much ice lies in a given thickness range (Hibler, 1980; Maykut, 1982).

A basic goal of sea-ice models is to predict the evolution of the ice thickness distribution (ITD) in time and space. In global climate models, gridcells are large (10−100 km) compared with typical ice floes (1 m−1 km), so we try to capture the variation in thermodynamic growth and melt rates by dividing the ice in each gridcell into thickness categories. This is a distribution in the mathematical sense, but in practice only a few categories are tracked, typically 5−20, with more categories assigned to thinner ice to better resolve growth rates. Both thickness and ice area fraction or concentration are tracked for each category, and thus the ITD essentially increases the thermodynamic resolution.

In mathematical terms, the ice thickness distribution function dh is the fractional area covered by ice in the thickness range (h, h + dh) at a given time and location. Thus, integrating the ice thickness distribution over all thicknesses, one must cover the entire area. The fundamental equation describing evolution of the ITD ties together the various processes in the sea-ice model:

(1)

where is the horizontal ice velocity, is the rate of thermodynamic ice growth, ψ is a ridging redistribution function (Thorndike and others, 1975) and L is lateral melting (Hibler, 1980). The four terms on the right-hand side describe different kinds of sea-ice transport: (1) horizontal transport in (x, y) space, (2) transport in thickness space due to ridging and other mechanical processes, (3) thickness transport due to thermodynamic growth and melting, and (4) replacement of ice with open water by lateral melting. To compute thermodynamic transport in thickness space, we need the ice growth rate f in each thickness category, and to solve the horizontal transport and ridging equations we must know the ice velocity Each of these pieces is discussed in the following sections.

3. Thermodynamics

Early sea-ice models were one-dimensional (1−D) thermodynamic descriptions of the melting and freezing of snow and ice, without reference to ice motion or deformation. The first sea-ice model was a mathematical description of ice growth (Stefan, 1891), work that is now widely referenced in sciences where moving boundaries are considered (‘Stefan problems’). We still use 1−D thermodynamic parameter- izations in GCMs, because the thermodynamic processes are still primarily vertical in nature, but the models are becoming more sophisticated as we learn how to better treat salinity and hydrological processes in the ice.

The first detailed thermodynamic sea-ice model was developed by Maykut and Untersteiner (1971). This model treats the ice as a homogeneous slab with energy fluxes at both surfaces and heat conduction in the interior. If the upper surface temperature is at the melting point (0°C), the net surface energy flux is used to melt snow or ice:

where q is the energy per unit volume required to melt the top surface material (either snow or ice), h is thickness, F _s is the sensible heat flux, F _l is the latent heat flux, F _L↓ is the incoming longwave flux, F _L↑ is the outgoing longwave flux, F _sw is the incoming shortwave flux, α is the shortwave albedo, and i ₀ is the fraction of absorbed shortwave flux that penetrates into the ice. A similar relation holds at the bottom of the ice.

The rate of temperature change in the ice interior is given by

(2)

where ρ_i is the sea-ice density, c_i(T, S) is the specific heat of sea ice, K_i(T, S) is the thermal conductivity of sea ice, l _pen is the flux of penetrating solar radiation at depth z, and z is the vertical coordinate, defined to be positive downward. The heat capacity and conductivity depend on both salinity and temperature of the bulk sea-ice material.

In the 1970s the Maykut and Untersteiner (1971) model was too complex for use in global climate models. Semtner (1976) therefore proposed two simpler versions. His three- layer model uses essentially the same surface flux and heat conduction equations as Maykut and Untersteiner (1971), but with pure-ice values of heat capacity and thermal conductivity, and with just two layers of ice and one of snow. Penetrating radiation is stored during the summer in a fictitious energy reservoir, and autumn cooling is delayed until the reservoir drains. Semtner showed that this model agrees closely with the more detailed model for a wide range of parameters, although its seasonal amplitude is somewhat greater. He also presented a zero-layer model, which assumes a linear temperature profile in the ice and hence neglects internal energy storage. This model can produce a reasonable mean thickness, but with substantial errors in the phase and amplitude.

Semtner (1984) pointed out that the zero-layer model is more sensitive than the three-layer model to increases in downward longwave radiation and cautioned against the use of such simple models in climate simulations. Nevertheless, many large-scale sea-ice models have used some variation of the zero-layer model (e.g. Parkinson and Washington, 1979; Hibler, 1980; Flato and Hibler, 1995), including most early global climate models (e.g. Washington and Meehl, 1984;Manabe and others, 1991). An innovation by Parkinson and Washington (1979) was to add a fractional area of open water to each model gridcell. When the energy flux to the lead is negative, ice grows in the lead and is merged with the ice already present. When the energy flux is positive, part of the energy is used to melt ice laterally and the remainder is used to warm the lead. Subsequent versions of this parameterization (e.g. Ebert and Curry, 1993) have applied some of the energy to bottom melting as well. Ebert and Curry (1993) also introduced parameterizations of brine pockets and melt ponds, as well as a sophisticated albedo parameterization with five surface types and four different spectral bands for each type.

Bitz and Lipscomb (1999) developed an energy-conserving treatment of internal heating and surface melting for thermodynamic models with brine-pocket parameteriz- ations. As sea ice nears its melting point, the energy absorbed internally is used not only to raise its temperature, but also to expand brine pockets. The porosity of warm ice can reach 20−30% in the upper 30 cm (Eicken and others, 1995) and even 40−50% under old melt ponds (Maykut and others, 1992). The energy required to melt the remaining ice is less than it would be for pure, nonporous ice. Bitz and Lipscomb (1999) therefore adjusted the latent heat of melting at the upper surface.

Thorndike (1992) introduced a ‘toy model’ of sea-ice thermodynamics that treats the ice very simply (the ice has a linear temperature profile, a single fixed albedo, and no snow, brine or leads) but includes an important longwave radiation feedback. With reasonable parameter choices it can support perennial ice and an ice-free ocean as stable states, but not seasonal ice.

The following subsections highlight several important aspects of sea-ice thermodynamics, including radiative forcing and surface conditions (snow or ponds), internal changes to the ice structure associated with heat and salt, formation of new ice, and finally, updating the ice thickness distribution.

3.1. Radiation and albedo

Sea-ice albedo has long been recognized as a critical aspect of the global heat balance. Sea ice is highly sensitive to atmospheric radiation (e.g. Maykut and Perovich, 1987; Perovich, 2003) because net radiative forcing on sea ice is a difference of two large terms, net downward shortwave and net upwelling longwave. A small change in eitherterm results in a comparatively large sea-ice response. Clouds, which represent a large uncertainty in GCMs, are critical for sea-ice simulation because they contribute to both terms, shading solar radiation and increasing downwelling longwave.

Many sea-ice models use a relatively simple albedo parameterization that specifies four albedo values: cold snow; warm, melting snow; cold, bare ice; and warm, melting ice. Other models use more complex formulations that take into account the ice and snow thickness, spectral band and other parameters. Barry (1996) provides an excellent compilation of modeling efforts and observations that indicate the complexity of the pack-ice albedo. Sea-ice albedos are often tuned to produce a realistic simulation of sea-ice extent, compensating for deficiencies in surface forcing (e.g. Losch and others, 2010).

Taylor and Feltham (2004) developed a two-stream radiation model that determines albedo from the seasonal evolution of bulk optical properties in three layers: sea ice, ponded meltwater and refrozen ponds. In CICE, solar radiation may be distributed within the ice column assuming exponential decay (Beer’s law) or via a multiple-scattering delta-Eddington radiative-transfer scheme (Briegleb and Light, 2007). In the latter case, ‘inherent’ optical properties for snow and sea ice, such as extinction coefficient and single-scattering albedo, are prescribed based on physical measurements. Reflected, absorbed and transmitted shortwave radiation (‘apparent’ optical properties, including albedo) are then computed for each snow and ice layer in a self-consistent way. Absorptive effects of inclusions such as dust and algae can also be simulated, along with radiative effects of melt ponds and other changes in physical properties (e.g. granularization associated with snow aging). As more of these higher-order processes are explicitly included in sea-ice models, their physical feedbacks on the ice and other constituents (such as biology) need to be accounted for through more realistic treatments of radiation.

3.2. Melt ponds and snow

There are several schemes for modeling the ponds that form on the surface of melting sea ice. The simplest and most widely used method does not involve tracking meltwater, but rather decreases the ice surface albedo under warm, melting conditions (e.g. Holland and others, 2006; Vancoppenolle and others, 2009a). This method simulates the critical radiative effect of melt ponds, but not their hydrological influence (e.g. the delay of internal ice cooling as ponds refreeze in the fall).

A second method, used in version 4 of the Community Climate System Model (CCSM4, http://www.cesm.ucar.edu/models/ccsm4.0), tracks melt-pond area and volume for each ice thickness category, but applies those quantities only for radiative calculations using the delta-Eddington multiplescattering scheme. The ponds are ‘virtual’ in the sense that all meltwater runs off into the ocean.

A third method involves more detailed modeling of melt- pond physics, including melting and freezing within the ponds (Ebert and Curry, 1993; Taylor and Feltham, 2004; Flocco and Feltham, 2007). In these studies, more complex radiative schemes are used to account for pond-induced albedo modifications. A portion of the meltwater runs off, while the rest is withheld until the melt pond melts fully through the ice or refreezes in the fall. The main difficulty with melt-pond modeling is lack of knowledge of the ice and snow topography.

Snow properties also need more attention in GCM sea- ice components. In the natural world, snow density varies; snow compacts in the wind, and snow grains change size as they age and are exposed to various environmental conditions. Flooded snow may refreeze to form snow ice (Jeffries and others, 1997) or may turn into a gap layer (Haas and others, 2001) in which biological communities thrive. Comprehensive snow models that include metamorphic processes are available (e.g. Jordan and others, 1999), but snow is usually treated quite simply in sea-ice models – for example, with constant density. Some snow processes such as flooding are included, but lack of observational data for validation hinders snow model development in sea-ice models. Also, more comprehensive snow models are generally formulated in terms of detailed representations of various properties (e.g. density and moisture content) within distinct snow layers. Although this is a sensible approach for snow on a fixed substrate, a fresh approach is needed to prevent the likely diffusion of layer properties when the underlying sea ice is advected.

3.3. Sea-ice microstructure

The processes by which ice forms contribute greatly to its final composition and structure. The thermodynamic model of Maykut and Untersteiner (1971) addresses the accretion of new ice onto the bottom of existing ice, resulting in columnar crystal structures. In turbulent, open-ocean conditions, frazil ice crystals collect at the surface to form a slush called grease ice, because of the appearance it gives the water’s surface. Further freezing under the influence of waves converts the grease ice into ‘pancakes’ up to ˜80cm thick (Ackley and Sullivan, 1994). This process results in a granular texture that can be differentiated from columnar ice in ice cores taken in the field. Although climate-scale sea-ice models incorporate both kinds of ice formation, the processes themselves are highly parameterized, and the resulting columnar or granular ice is combined into one uniform ice type in the simulated ice column. More detailed parameterizations are under development for use in climate models (e.g. Dai and others, 2004; Skogseth and others, 2009).

In general, sea ice is quite heterogeneous, in large part because of its salt content. As ice crystals grow in sea water, brine is captured among the crystals, altering the heat capacity and conductivity of the ice (Feltham and others, 2006). During much of the year the ice is cold at the top and warm at the bottom, which can induce convection. In addition, as pure ice crystals form from brine, they leave behind even saltier brine, which is dense. Thus the brine tends to flow out of the bottom of growing sea ice, while less saline sea water flows into it, a process referred to as gravity drainage (Notz and Worster, 2009). When the ice becomes cold, brine left in the ice structure collects in brine pockets, and when the ice warms again, the brine pockets expand and interconnect, creating brine channels. Melt ponds form in depressions on the surface of the ice and can drain through the brine channels when the ice becomes warm and permeable. This flushing can effectively clean the ice of salt, nutrients and other inclusions (Eicken and others, 2004). These changes within the ice affect the albedo, conductivity and biogeochemical processes, thereby playing a role in climate change. The effort to learn how these processes interact, and determine their importance in the polar environment, constitutes a vibrant area of sea-ice research.

In many sea-ice models used for climate study, the sea- ice salinity assumes a fixed value that is used for freshwater and salt exchanges at the ice-ocean interface, ensuring salt conservation in the coupled modeling context. The salinity used internally to the sea-ice model may vary independently of the coupling value, however. Some models allow the internal salinity to vary in time (e.g. Schramm and others, 1997), while others (e.g. Bitz and Lipscomb, 1999) assume a variable salinity profile that is constant in time and representative of multi-year sea ice. This profile is less appropriate for younger ice types found in the Antarctic and, increasingly, in the Arctic, spurring interest in improving the models.

Cox and Weeks (1988) introduced an empirical formulation of winter sea-ice desalination in a simple model. Using this approach, Vancoppenolle and others (2005) developed a 1-D model that successfully simulates the desalination of Arctic sea ice as it first grows and then transitions from first- year to multi-year ice. A simplification of this model was used in a full sea-ice model, LIM3, to establish the importance of changes in sea-ice microstructure for its salinity and evolution, including feedbacks associated with changes in ocean stratification (Vancoppenolle and others, 2009a, b).

3.4. Transport in thickness space

Changes to the ice thickness through all of these processes must be applied to determine the new ice thickness distribution. Early sea-ice models made the simple but unrealistic assumption of a single, uniform ice thickness in each gridcell. With the introduction of an ice thickness distribution, a method must be provided for transferring sea ice among thickness categories, that is, for solving

(3)

which is obtained from Equation (1) by neglecting the ridging, horizontal transport and lateral melt terms. Once the surface growth and melt rates f are determined for a given thickness category, the ice thickness can be adjusted. However, the new thickness may lie in a different category.

Early methods for solving Equation (3) fixed the mean thickness in each category while varying the category areas (Hibler, 1980). These schemes were highly diffusive, requiring many categories of both ridged and undeformed ice to improve accuracy (Flato and Hibler, 1995). Other schemes allowed the thicknesses to vary, adding categories when needed and later merging them (Bjork, 1992; Schramm and others, 1997). A similar method allows the ice thickness to vary within category boundaries, occasionally transferring all the ice in a category to a neighboring category when the thickness changes sufficiently (Bitz and others, 2001). This ‘delta-function’ scheme can cause sudden changes in properties such as ice strength if the ITD is not well resolved.

CICE uses the remapping scheme of Lipscomb (2001), which can be thought of as a 1-D analog of the two-dimensional (2-D) incremental remapping scheme described in section 4.3 below. In the 1-D scheme, the ‘velocity’ is the thermodynamic growth or melt rate, and the ice ‘advects’ in thickness space. Thickness categories are represented as Lagrangian gridcells whose boundaries are projected forward in time as the ice grows or melts. The thickness distribution function g is approximated as a linear

remapped onto the original categories. This method is numerically smooth and not too diffusive, requiring only a few (approximately five) thickness categories to simulate large-scale properties of the ice pack.

3.5. Lateral melting

Rather than changing the actual ice thickness, lateral melting affects the ice thickness distribution g directly by changing the ice area (Equation (1)). Solar radiation warms leads and other open water within the ice pack, and some of this heat melts the edges of ice floes. The fraction of available ocean heat used for lateral melting will depend on floe size and geometry, which are not generally known in large-scale sea-ice models used in GCMs. In sea-ice models with ice thickness distributions, melting through of ice in thinner categories can partially account for this process (Bitz and others, 2001; Vancoppenolle and others, 2009b). In other models (e.g. Hunke and Lipscomb, 2010), a fraction of the total ocean heat available for melting ice is apportioned to lateral melting following Maykut and Perovich (1987) and Steele (1992).

4. Dynamics

The first two terms on the right-hand side of Equation (1) are due to ice dynamics. In the Arctic, sea ice moves from Siberia toward northern Greenland in the Transpolar Drift Stream. Some ice piles up near the Canadian coast, some circulates in the Beaufort Gyre, and each year about 10% of the pack flows through Fram Strait and the Canadian Arctic Archipelago to the North Atlantic, where it melts. Sea ice is relatively fresh (up to 25 psu; Untersteiner, 1968) compared to sea water, and so the melt dilutes the dense, salty water brought northward by the Gulf Stream, potentially disrupting the density-driven movement of global ocean currents. In the Antarctic, winds flowing northward from the continental ice sheet drive a freshwater conveyor in which sea ice formed in the south is carried away from the coast, leaving saline water (produced by brine rejection of the freezing ice) to descend along the continental shelf, and transporting fresh meltwater northward. This process creates some of the densest sea water on Earth (Aagaard and others, 1985; Foldvik and others, 2004).

In order to simulate this ice movement, sea-ice models generally include equations for the momentum, rheology (a constitutive law that describes the material properties of the ice), transport, and mechanical deformation, often referred to as ‘ridging’. These equations must take into account some basic properties of ice motion. Because sea ice is highly fractured, it can diverge easily. However, it is still an essentially rigid material and therefore resists convergent forcing. If the forcing is strong enough, the ice will break, forming ridges along floe edges or other weak points in the pack. These considerations lead to a momentum equation and constitutive law, used to determine the ice velocity. The deformation or rates of strain, which are spatial derivatives of the velocity components, are also important. They control the amount of ridging, which in turn determines the amount of open water that is exposed to the atmosphere.

4.1. Momentum

The force balance per unit area in the ice pack is given by a 2-D momentum equation, obtained by integrating the three-dimensional equation with respect to the thickness of the ice in the vertical direction (Hibler, 1979):

(4)

where m is the combined mass of ice and snow per unit area and and are wind and ocean stresses, respectively. The strength of the ice is represented by the internal stress tensor σ_ij , and the last two terms on the right-hand side are stresses due to Coriolis effects and the sea surface slope.

Near the ice edge, where floe concentrations are low, the primary balance in Equation (4) is among the inertial (left- hand side), wind stress and ocean stress terms. Early dynamic models used a free-drift approach, in which the Coriolis term was included and inertia was taken to be zero (Fel’zenbaum, 1961; Bryan and others, 1975; Manabe and others, 1979; Parkinson and Washington, 1979), but lack of internal ice stress allowed too much motion in regions where the ice pack should be relatively stiff.

In the ice interior, the primary force balance is given by the internal ice stress, wind stress, and ocean stress. The internal stress, σ, is given by a constitutive equation that describes the rheology of the material. While the ice motion itself is critical for transport of ice properties and tracers, the rheology is strongly linked with deformational processes (Geiger and others, 1998).

4.2. Rheology

Early sea-ice rheologies treated the ice as a Newtonian viscous fluid (Campbell, 1965), a linear viscous fluid (Hibler, 1974; Hibler and Tucker, 1979) or a plastic material. AIDJEX researchers proposed an elastic-plastic rheology (Coon and others, 1974; Pritchard and others, 1977), and several other nonlinear plastic rheologies have since been studied (see Feltham, 2008, for a review). A nonlinear viscous-plastic (VP) constitutive law proposed by Hibler (1979) has become the standard rheological model. Most plastic rheologies describe a material that resists deformation (often through an elastic response) until the forcing stress becomes large enough to cause it to deform, or ‘yield’. The essential assumption in the VP approach is that sea ice is in a continual state of yielding, and only when the strain rates approach zero (and effective viscosities become infinite) is a viscous ‘creep’ behavior imposed.

Simulations with more complicated rheologies than the standard elliptical yield curve (Hibler, 1979) – such as teardrop (Coon and others, 1974), sine wave lens (Bratchie, 1984), Mohr-Coulomb (Tremblay and Mysak, 1997) and square (Ip and others, 1991) shapes – show that the rheology can have a significant effect on long-term simulations of ice drift and deformation (Ip and others, 1991; Geiger and others, 1998). Changing the aspect ratio of the elliptical yield curve has also been shown to significantly alter the pattern of sea-ice thickness (Miller and others, 2005).

The enormous range of effective VP viscosities severely limits the time-step for stability of an explicit numerical scheme, especially in regions where the ice is relatively rigid. When run on parallel computers, implicit solution methods such as successive over-relaxation (Hibler, 1979) and line relaxation (Holland and others, 1993; Oberhuber, 1993; Zhang and Hibler, 1997) entail a great deal of communication between processors and are difficult to parallelize efficiently, unlike explicit methods. Simpler versions of the constitutive model, such as free-drift descriptions with no ice interaction and cavitating-fluid models in which the ice has no resistance to shear forces (Nikiforov and others, 1967; Flato and Hibler, 1990, 1992), are more tractable numerically, but the model behavior is sensitive to these simplifications (Holland and others, 1993). Hunke and Dukowicz (1997, 2002) modified the VP model by adding elastic behavior, thus realizing large gains in numerical efficiency by permitting a fully explicit implementation with an acceptably long time-step. This elastic- viscous-plastic (EVP) model reduces to the original VP model on long timescales and was shown to be more accurate for transient behavior (Hunke and Zhang, 1999). The recent Jacobian-free Newton-Krylov approach (Lemieux and others, 2010) promises computationally efficient VP solutions for GCMs.

At moderately high resolution, the VP model captures an essential property of the large-scale ice pack: the presence of long, linear features. Sea ice tends to move as large plates with long faults related to shearing. The coarse-resolution (1°) grids typically used in GCMs cannot resolve these features, but at higher resolutions the plastic aspect of the model causes the deformation to be concentrated in long, narrow bands of high shear. Closer examination using modern satellite data, however, shows that the VP model does not reproduce these features very accurately (Lindsay and others, 2003; Coon and others, 2007).

Many efforts to move beyond the classic VP model toward more physical descriptions of sea-ice momentum and material properties have been undertaken or are under way. For instance, the assumption of isotropy was deemed reasonable when gridcells were on the order of 100 km, but as mesh resolution increases, anisotropic rheology is more realistic (Wilchinsky and Feltham, 2004, 2006; Coon and others, 2007). Newer sea-ice dynamics models include a particle-in-cell method (Flato, 1993), a granular, Lagrangian approach (Hopkins, 2004; Hopkins and Thorndike, 2006), smoothed-particle hydrodynamics (Gutfraind and Savage, 1997; Lindsay and Stern, 2004), a continuum elastic- decohesive rheology (Schreyer and others, 2006; Sulsky and others, 2007) and an elasto-brittle progressive damage rheology (Girard and others, 2010a, b). Because the dynamics component often uses a majority of the computational time in sea-ice models, an important goal is to develop an algorithm that captures the relevant physical processes at the desired scales while still remaining computationally feasible for use in GCMs.

4.3. Transport

Given the ice velocity field, the next step is to move the ice and overlying snow from one grid location to another while conserving properties such as mass and internal energy. The densities of ice and snow are often assumed to be constant to enable numerical conservation of volume. We therefore solve the continuity or transport equation,

(5)

for the fractional ice area, a _i, in each thickness category n. Equation (5) describes the conservation of ice area under horizontal transport. It is obtained from Equation (1) by discretizing g and neglecting the second, third and fourth terms on the right-hand side, which are treated separately.

There are similar conservation equations for ice volume, snow volume, ice energy and snow energy. In addition, there may be equations describing the transport of tracers T_n .

These equations typically have one of the following three forms

(6)

(7)

(8)

where v _in and v _sn are the ice and snow volume, respectively, in category n. Note that the conservation equation for ice volume can be cast as Equation (6), considering ice thickness as a tracer following ice area.

These generic equations are solved in many fluid- dynamical applications. A desirable numerical method for solving transport equations would be conservative, stable, accurate, sign- and monotonicity-preserving, and computationally efficient. In addition, tracers at a given grid location should remain bounded by their values in neighboring gridcells at the previous time-step, a property referred to as compatibility. The simplest method commonly used for sea- ice transport is the first-order upwind, or donor-cell, scheme. This method has most of the characteristics just mentioned but is also very diffusive, reducing accuracy. Nevertheless, because of their efficiency, upwind schemes often are used for sea-ice transport in climate models. The simplest second- order-accurate methods are centered-in-space schemes, but these can be highly oscillatory and are generally not sign-preserving. Oscillations can be suppressed by adding artificial diffusion terms, at the cost of reduced accuracy.

The second-order moment (SOM) scheme of Prather (1986) is used for sea-ice transport (Merryfield and Holloway, 2003; Vancoppenolle and others, 2009a). SOM is a modified upwind scheme that reduces diffusion by transporting not only the mean fields, but also their first- and second-order moments: a total of six fields in two dimensions. This scheme nearly preserves monotonicity for conserved fields (ice area and volume), but not for tracers. It is third-order accurate in space except where limited to prevent overshoots and undershoots. SOM is relatively inexpensive when used to transport just a few fields but becomes costly when applied to large numbers of tracers. Russell and Lerner (1981) developed a similar method, the linear-upstream or slopes scheme, which transports first- order but not second-order moments and thus is second- order accurate in space.

CICE uses a 2-D incremental remapping scheme (Lipscomb and Hunke, 2004), in which the transport equation is solved by projecting model gridcells backward in time along Lagrangian trajectories. Scalar fields are reconstructed and integrated over the Lagrangian departure regions, then remapped onto the grid. When the conserved fields are constructed appropriately, remapping has all the desirable features listed earlier. The spatial accuracy depends on the accuracy of the reconstruction: if the fields are constant within each gridcell, remapping is a first-order scheme, whereas linear fields produce a second-order method. Remapping preserves monotonicity when the field gradients are limited appropriately, although this limiting may reduce the accuracy locally. In terms of efficiency, the method is relatively expensive for the basic state variables. Additional tracers are inexpensive, however, because geometrical quantities for departure cells are computed once and then are reused for tracers with little extra computational expense. With many new biological and diagnostic tracer fields, incremental remapping becomes more efficient relative to other transport schemes in ice models.

4.4. Mechanical redistribution and ice strength

When densely packed ice converges within a gridcell, the solution of Equation (5) can result in an ice fractional area greater than 1. Early sea-ice models simply limited the fractional area to be ≤1, redistributing the extra ice within a gridcell by increasing the thickness. Modern models use a mathematical, continuum description of mechanical redistribution, or ridging, as represented by the second term on the right-hand side of Equation (1). This term converts thinner ice to thicker ice under convergence and shear, ensuring that the total ice area does not exceed the gridcell area.

Mechanical redistribution occurs in two steps. First, a participation function defines which portion of the ice thickness distribution participates in ridging. In the original formulation of Thorndike and others (1975), only the thinnest 15% of the ice pack would participate, but discontinuous derivatives caused numerical instabilities in high-resolution models. Treating the participation function as an exponential has mitigated this problem (Lipscomb and others, 2007).

The second piece is a redistribution function, which describes how the ice participating in ridging is redistributed among other thickness categories, based on the initial thickness and the desired height-to-area ratio of the ridge. Hibler (1980) specified that ridges are triangular, with a maximum thickness proportional to the square root of the thickness of the ridging ice. This parameterization agreed better with observations than the single-thickness ridges of Thorndike and others (1975). However, Amundrud and others (2004) found that this scheme produced too much deep ridged ice and suggested modifications to the assumed keel shape to improve the distribution of ice draft. CICE assumes an exponential distribution of ridge thicknesses, consistent with observations of thick ice. The correct form of the redistribution function remains an open question recognized since AIDJEX (Untersteiner, 1980).

In some models the ridging scheme is closely related to the ice strength, or pressure P, which may be defined in terms of the frictional energy loss and potential energy production needed to obtain the new ridge distribution (Rothrock, 1975). The ice strength in turn affects the ice rheology and momentum through the appearance of horizontal gradients of pressure in the constitutive law. Here the ice strength is the compressive stress below which the ice is considered rigid and at which it fails, or yields. Many climate models, however, use an empirical description for P given as a function of ice area fraction a and thickness h (Hibler, 1979),

where P* and a* are constants. This form for the ice strength captures the essential characteristics of weak, low-concentration ice and very strong, high-concentration ice, with a linear dependence on ice thickness.

Although they have not been used in climate models because of their complexity and expense, more detailed models of ridging are available. For example, Parmerter and Coon (1972) constructed a kinematic model that simulates the displacement of rubble in a lead as the surrounding thick ice floes converge. Rubble can accumulate above and below the thick floes, causing them to bend and eventually break. Parmerter and Coon found that ridges grow to a limiting height, with additional ridging extending laterally. Hopkins and others (1991) adopted a similar conceptual approach but used a particle simulation to explicitly calculate the dynamical interactions of the ice floes and rubble. Although these approaches provide a more complete understanding of the energetics, such detailed simulations of ridging are not yet feasible for GCMs.

4.5. Dynamic coupling

Wind stress is arguably the primary external forcing mechanism for the ice motion (Equation (4)), although the ice- ocean stress, Coriolis force, and slope of the ocean surface are also important (Steele and others, 1997). Coupling between sea-ice models and atmospheric models or data generally employs a quadratic form for the wind stress,

where ρ_a and denote atmosphere density and wind, respectively, and D_a is a drag coefficient, usually taken to be constant. This can be written in the equivalent form

in which the friction velocity u* is defined via MoninObukov similarity theory for turbulent exchange (Andreas and others, 2004). In this case, u ^* is computed simultaneously and consistently with analogous coefficients for the turbulent fluxes of sensible and latent heat.

Likewise, a quadratic ice-ocean stress term is typically applied at the bottom of the ice,

Again, the exchange coefficient D _w may be constant or defined in terms of a friction velocity. Here, and are the ocean and ice velocities, respectively, θ is the turning angle between geostrophic and surface currents, ρ_w is the density of sea water, and is the vertical unit vector. A turning angle θ is not necessary if the top ocean model layers are able to resolve the Ekman spiral in the boundary layer. For a sea-ice undersurface roughness of 0.03 m, θ takes values up to about 25° for ice speeds greater than 0.05 m s^-1 (McPhee, 2008).

The parameterization for the wind and ice-ocean stress terms must contain the ice concentration as a multiplicative factor to be consistent with the theory of free drift in regions of low ice concentration (Connolley and others, 2004).

In some models, a stress equal and opposite to is applied to the ocean surface, thus conserving momentum. In this view, sea ice does not displace the underlying ocean but instead ‘levitates’ on the ocean surface. In reality, however, sea ice floats in its own liquid phase, and representing the ocean and ice as separate models is not at all natural. Thus, while is usually used for forcing sea-ice models, other forms for the ice-ocean stress are often applied in ocean models. These alternative approaches generally strive to describe the combined ice and ocean mass transport (e.g. Hibler and Bryan, 1987; Heil and Hibler, 2002). In an intermediate approach, the surface pressure felt by the ocean due to the ice may be included by using a z* coordinate system in the ocean model (Campin and others, 2008). Buoyancy affects the ice-ocean stress as well (McPhee, 1994); the interface of sea-ice and ocean modeling will benefit from further investigation.

5. Beyond The Physical Model

In addition to the core dynamic and thermodynamic components, sea-ice models include many other parameterizations, such as the surface fluxes of heat and moisture, flooding of ice floes, and so on. A current thrust in sea-ice modeling is the addition of new tracers to help elucidate model processes or to compare with new measurements. For example, estimates of sea-ice age are available from remotesensing data and are used to estimate changes in sea-ice volume (Maslanik and others, 2007). Age tracers are being added to sea-ice models, but there are different definitions of ‘age’, both for observations and for models, which complicate the comparisons (Hunke and Bitz, 2009; Vancoppenolle and others, 2009a). Likewise, it is difficult to compare in situ observations with model output for other sea-ice properties, such as the ratio of level to deformed ice. Model-data comparisons will continue to be a rich area for research.

Another new thrust for climate modeling is biogeochemistry, which in sea ice depends critically on the ice microstructure and salinity profile. Several groups are working toward coupling sea-ice biological and chemical descriptions into GCMs. For the Arctic, the effort focuses on bottom ice communities (Lavoie and others, 2005; Jin and others, 2006), but for the Southern Ocean the full sea-ice column must be included (Vancoppenolle and others, 2009b). Related work involves the vertical transport and cycling of quantities such as aerosols (e.g. Bailey and others, 2010) and gases (Nomura and others, 2010) that pass gradually through the ice and can modify oceanic or atmospheric chemistry. These inclusions interact with sea ice through radiative transfer and other physical processes such as flushing by meltwater.

Finally, every physical model rests on an infrastructure that includes input and output files, a grid (structured, unstructured or a combination), an approach for using parallel computer architectures, and so on. For climate models, computationally efficient algorithms and coding practices are critical for enabling multiple-century simulations at fine grid resolutions desired for accuracy.

6. Summary

Sea ice has long been recognized as an important element of the global climate system. The threat of high-latitude climate change (e.g. Holland and Bitz, 2003) provides impetus to improve existing sea-ice models. Most large-scale, physical sea-ice processes are well understood and well represented in models; for example, the basic thermodynamic description has been known for 40 years (Maykut and Untersteiner, 1971), although it was implemented in GCMs only within the past decade. Likewise, a relatively simple approach for sea-ice dynamics is 30 years old (Hibler, 1979, 1980) and finally appearing (modified) in most modern GCMs. These thermodynamic and dynamic models capture the first-order behavior of sea ice in the climate system.

Sea-ice model development now follows two paths, both arguably addressing higher-order effects: (1) more precise descriptions of physical processes and characteristics, and (2) extensions of the model for ‘Earth system’ simulations with biogeochemistry.

For example, multiple-scattering radiation schemes that take into account the effects of melt ponds and sea-ice inclusions provide better estimates of reflected and absorbed radiation, and of temperature profiles in the ice. New approaches for determining the evolution of salinity and, more generally, the sea-ice microstructure, are important for modeling biological and chemical species in sea ice. Efforts continue to improve the representation of other processes that influence the pack ice evolution, such as the development of frazil ice into pancakes and eventually a solid ice cover, and melt ponds. Various anisotropic descriptions of the rheology promise more accurate representation of linear kinematic features and ice deformation. In addition to the research development currently under way, mechanical deformation and snow processes both need closer examination and development for sea-ice components of GCMs.

Stronger interactions and feedbacks among the many model elements will become important as climate modeling progresses. For instance, inclusions of dust, aerosols and biology affect solar absorption and the sea-ice microstructure, and can thus contribute to faster melting and weakening of the ice pack.

Likewise, as climate component models (e.g. atmosphere, land ice, sea ice and ocean) become more intimately coupled, careful coupling strategies must be developed, and sea-ice components will need to be fully validated within the coupled context. Fresh numerical approaches and algorithm improvements play a critical role in the development process, as climate models continue to push the limits of computational power.