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, Reference UntersteinerUntersteiner (1964) developed a detailed model of heat conduction and surface growth and melting of sea ice. A few years later, Reference Maykut and UntersteinerMaykut 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. Reference Parmerter and CoonParmerter and Coon, 1972; Reference Coon, Maykut, Pritchard, Rothrock and ThorndikeCoon and others, 1974; Reference HiblerHibler, 1979; Reference Hopkins, Hibler,III and FlatoHopkins 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. Reference Hibler and BryanHibler 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 (Reference Bitz, Holland, Weaver and Eby.Bitz and others, 2001), Thickness and Enthalpy Distribution (TED; Reference Zhang and RothrockZhang and Rothrock, 2001) and Louvain-la-Neuve (LIM; Reference Vancoppenolle, Fichefet, Goosse, Bouillon, Madec and Morales MaquedaVancoppenolle 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 (Reference Kwok and RothrockKwok and Rothrock, 2009), due in part to increased export (Reference Rigor and WallaceRigor and Wallace, 2004; Reference Hunke and BitzHunke 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 (Reference HiblerHibler, 1980; Reference MaykutMaykut, 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:

where is the horizontal ice velocity, is the rate of thermodynamic ice growth, ψ is a ridging redistribution function (Reference Thorndike, Rothrock, Maykut and ColonyThorndike and others, 1975) and L is lateral melting (Reference HiblerHibler, 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 (Reference StefanStefan, 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 Reference Maykut and UntersteinerMaykut 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

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 Reference Maykut and UntersteinerMaykut and Untersteiner (1971) model was too complex for use in global climate models. Reference SemtnerSemtner (1976) therefore proposed two simpler versions. His three- layer model uses essentially the same surface flux and heat conduction equations as Reference Maykut and UntersteinerMaykut 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.

Reference SemtnerSemtner (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. Reference Parkinson and WashingtonParkinson and Washington, 1979; Reference HiblerHibler, 1980; Reference Flato and Hibler,IIIFlato and Hibler, 1995), including most early global climate models (e.g. Reference Washington and MeehlWashington and Meehl, 1984;Reference Manabe, Stouffer, Spelman and BryanManabe and others, 1991). An innovation by Reference Parkinson and WashingtonParkinson 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. Reference Ebert and CurryEbert and Curry, 1993) have applied some of the energy to bottom melting as well. Reference Ebert and CurryEbert 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.

Reference Bitz and LipscombBitz 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 (Reference Eicken, Lensu, Leppäranta, Tucker, Gow and SalmelaEicken and others, 1995) and even 40−50% under old melt ponds (Reference Maykut, Grenfell and WeeksMaykut and others, 1992). The energy required to melt the remaining ice is less than it would be for pure, nonporous ice. Reference Bitz and LipscombBitz and Lipscomb (1999) therefore adjusted the latent heat of melting at the upper surface.

Reference ThorndikeThorndike (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.

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; Reference UntersteinerUntersteiner, 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 (Reference Aagaard, Swift and CarmackAagaard and others, 1985; Reference FoldvikFoldvik 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 (Reference HiblerHibler, 1979):

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 (Reference Fel’zenbaumFel’zenbaum, 1961; Reference Bryan, Manabe and PacanowskiBryan and others, 1975; Reference Manabe, Bryan and SpelmanManabe and others, 1979; Reference Parkinson and WashingtonParkinson 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 (Reference Geiger, Hibler,III and AckleyGeiger and others, 1998).

4.2. Rheology

Early sea-ice rheologies treated the ice as a Newtonian viscous fluid (Reference CampbellCampbell, 1965), a linear viscous fluid (Reference HiblerHibler, 1974; Reference HiblerHibler and Tucker, 1979) or a plastic material. AIDJEX researchers proposed an elastic-plastic rheology (Reference Coon, Maykut, Pritchard, Rothrock and ThorndikeCoon and others, 1974; Reference Pritchard, Coon and McPheePritchard and others, 1977), and several other nonlinear plastic rheologies have since been studied (see Reference FelthamFeltham, 2008, for a review). A nonlinear viscous-plastic (VP) constitutive law proposed by Reference HiblerHibler (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 (Reference HiblerHibler, 1979) – such as teardrop (Reference Coon, Maykut, Pritchard, Rothrock and ThorndikeCoon and others, 1974), sine wave lens (Reference BratchieBratchie, 1984), Mohr-Coulomb (Reference Tremblay and MysakTremblay and Mysak, 1997) and square (Reference Ip, Hibler,III and FlatoIp and others, 1991) shapes – show that the rheology can have a significant effect on long-term simulations of ice drift and deformation (Reference Ip, Hibler,III and FlatoIp and others, 1991; Reference Geiger, Hibler,III and AckleyGeiger 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 (Reference Miller, Laxon and FelthamMiller 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 (Reference HiblerHibler, 1979) and line relaxation (Reference Holland, Mysak, Manak and OberhuberHolland and others, 1993; Reference OberhuberOberhuber, 1993; Reference Zhang and Hibler,IIIZhang 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 (Reference Nikiforov, Gudkovich, Yefimov and RomanovNikiforov and others, 1967; Reference Flato and Hibler,IIIFlato and Hibler, 1990, Reference Flato and Hibler,III1992), are more tractable numerically, but the model behavior is sensitive to these simplifications (Reference Holland, Mysak, Manak and OberhuberHolland and others, 1993). Reference Hunke and DukowiczHunke and Dukowicz (1997, Reference Hunke and Dukowicz2002) 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 (Reference Hunke and ZhangHunke and Zhang, 1999). The recent Jacobian-free Newton-Krylov approach (Reference LemieuxLemieux 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 (Reference Lindsay, Zhang and RothrockLindsay and others, 2003; Reference Coon, Kwok, Levy, Pruis, Schreyer and SulskyCoon 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 (Reference Wilchinsky and FelthamWilchinsky and Feltham, 2004, 2006; Reference Coon, Kwok, Levy, Pruis, Schreyer and SulskyCoon and others, 2007). Newer sea-ice dynamics models include a particle-in-cell method (Reference FlatoFlato, 1993), a granular, Lagrangian approach (Reference HopkinsHopkins, 2004; Reference Hopkins and ThorndikeHopkins and Thorndike, 2006), smoothed-particle hydrodynamics (Reference R and SavageGutfraind and Savage, 1997; Reference Lindsay and SternLindsay and Stern, 2004), a continuum elastic- decohesive rheology (Reference Schreyer, Sulsky, Munday, Coon and KwokSchreyer and others, 2006; Reference Sulsky, Schreyer, Peterson, Kwok and CoonSulsky and others, 2007) and an elasto-brittle progressive damage rheology (Reference Girard, Amitrano and WeissGirard and others, 2010a, Reference Girard, Bouillon, Weiss, Amitrano, Fichefet and Legatb). 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,

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

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 Reference PratherPrather (1986) is used for sea-ice transport (Reference Merryfield and HollowayMerryfield and Holloway, 2003; Reference Vancoppenolle, Fichefet, Goosse, Bouillon, Madec and Morales MaquedaVancoppenolle 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. Reference Russell and LernerRussell 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 (Reference Lipscomb and HunkeLipscomb 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 Reference Thorndike, Rothrock, Maykut and ColonyThorndike 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 (Reference Lipscomb, Hunke, Maslowski and JakackiLipscomb 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. Reference HiblerHibler (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 Reference Thorndike, Rothrock, Maykut and ColonyThorndike and others (1975). However, Reference Amundrud, Melling and IngramAmundrud 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 (Reference UntersteinerUntersteiner, 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 (Reference RothrockRothrock, 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 (Reference HiblerHibler, 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, Reference Parmerter and CoonParmerter 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. Reference Hopkins, Hibler,III and FlatoHopkins 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 (Reference Steele, Zhang, Rothrock and SternSteele 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 (Reference Andreas, Jordan and MakshtasAndreas 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 (Reference McPheeMcPhee, 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 (Reference Connolley, Gregory, Hunke and McLarenConnolley 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. Reference Hibler and BryanHibler and Bryan, 1987; Reference P and HiblerHeil 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 (Reference Campin, Marshall and FerreiraCampin and others, 2008). Buoyancy affects the ice-ocean stress as well (Reference McPheeMcPhee, 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 (Reference Maslanik, Fowler, Stroeve, Drobot and ZwallyMaslanik 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 (Reference Hunke and BitzHunke and Bitz, 2009; Reference Vancoppenolle, Fichefet, Goosse, Bouillon, Madec and Morales MaquedaVancoppenolle 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 (Reference Lavoie, Denman and MichelLavoie and others, 2005; Reference JinJin and others, 2006), but for the Southern Ocean the full sea-ice column must be included (Reference Vancoppenolle, Fichefet, Goosse, Bouillon, Madec and Morales MaquedaVancoppenolle and others, 2009b). Related work involves the vertical transport and cycling of quantities such as aerosols (e.g. Reference BaileyBailey and others, 2010) and gases (Reference Nomura, Yoshikawa-Inoue, Toyota and ShirasawaNomura 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. Reference Holland and BitzHolland 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 (Reference Maykut and UntersteinerMaykut 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 (Reference HiblerHibler, 1979, Reference Hibler1980) 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.