## Abstract

The granular sea ice model (GRAN) from Tremblay and Mysak is converted from Cartesian to spherical coordinates. In this conversion, the metric terms in the divergence of the deviatoric stress and in the strain rates are included. As an application, the GRAN is coupled to the global Earth System Climate Model from the University of Victoria. The sea ice model is validated against standard datasets. The sea ice volume and area exported through Fram Strait agree well with values obtained from in situ and satellite-derived estimates. The sea ice velocity in the interior Arctic agrees well with buoy drift data. The thermodynamic behavior of the sea ice model over a seasonal cycle at one location in the Beaufort Sea is validated against the Surface Heat Budget of the Arctic Ocean (SHEBA) datasets. The thermodynamic growth rate in the model is almost twice as large as the observed growth rate, and the melt rate is 25% lower than observed. The larger growth rate is due to thinner ice at the beginning of the SHEBA period and the absence of internal heat storage in the ice layer in the model. The simulated lower summer melt is due to the smaller-than-observed surface melt.

## 1. Introduction

Sea ice dynamics plays an important role in shaping the sea ice cover in polar regions. In the last few decades, several models have been developed to represent the dynamics of sea ice. Crucial to the model representation of dynamics is the formulation of the rheology, that is, the relationship between applied stresses and the resulting deformations.

In these models, sea ice has been modeled either as a continuum (e.g., Coon et al. 1974) or as a collection of discrete particles (e.g., Hopkins 1996). Observations during the Arctic Ice Dynamics Joint Experiment (AIDJEX) campaign suggested that on a scale of about 100 km, there are enough randomly distributed cracks, ridges, and leads to allow sea ice to be treated as an isotropic material and as a continuous medium (Coon et al. 1974). Until recently, most sea ice models were based on the continuum approach and assumed that sea ice is an isotropic material. Recent studies, however, are questioning these assumptions (e.g., Coon et al. 2006), and anisotropic/discontinuous models are being developed (e.g., Hibler and Schulson 2000; Schreyer et al. 2007).

The models assuming isotropy and a continuous medium can be divided into two categories: (i) the elastic–plastic models (e.g., Coon et al. 1974; Pritchard 1975) and (ii) the viscous–plastic models (e.g., Hibler 1979; Hunke and Dukowicz 1997; Tremblay and Mysak 1997; Zhang and Rothrock 2005).

In the elastic–plastic approach, sea ice undergoes small elastic strains under low stress conditions, while plastic deformations occur when the stresses reach critical values defined by a yield curve (Coon et al. 1974). Since the elastic part of the deformation is reversible, the history of the strain has to be known. This usually requires a Lagrangian description of the motion. Furthermore, the time step has to be small in order to resolve elastic waves. In the viscous–plastic approach, the elastic strains are neglected (Hibler 1977). When the deformations are small, sea ice behaves as a very viscous fluid (i.e., exhibits creeping flow). Since the viscous–plastic approach can be used with an Eulerian description and is numerically easier to implement, these models are now widely used for large-scale sea ice simulations. Nevertheless, recent work shows promising developments with an elastic–decohesive model (Schreyer et al. 2006) using an efficient numerical scheme in a Lagrangian framework (Sulsky et al. 2007), which is argued to be as numerically efficient as the Hunke and Dukowicz (1997) model.

To complete the formulation of the viscous–plastic rheology, a yield curve, which marks the transition between the viscous phase and the plastic phase, and a flow rule, which describes the relation between the strain rates and the stresses on the yield curve, have to be specified. Owing to the difficulties in measuring in situ large-scale stresses, the yield curve and the flow rule that are the most realistic on geophysical scales remain unclear. For this reason, various yield curves have been used in sea ice modeling. For instance, current models use an ellipse (Hibler 1979; Hunke and Dukowicz 1997), a lens (Zhang and Rothrock 2005), a square (Ip et al. 1991), a teardrop (Zhang and Rothrock 2005), a line (Flato and Hibler 1992), or a triangle (Ip et al. 1991; Flato and Hibler 1992; Tremblay and Mysak 1997).

For an elastic–plastic rheology, the principle of plastic irreversibility states that the yield curve should be convex and that the resulting deformation should be normal to the yield curve (Goodier and Hodge 1958). Although it is not required for a viscous–plastic rheology to ensure proper energy dissipation (Hibler and Schulson 2000), most of the viscous–plastic models use a normal flow rule (e.g., Hibler 1979). But others, like the ones based on a Mohr–Coulomb failure criterion (triangle), use a nonnormal flow rule (e.g., Ip et al. 1991; Tremblay and Mysak 1997).

The most widely used viscous–plastic model is the one developed by Hibler (1979), and modifications thereof. This model uses an elliptical yield curve with a normal flow rule. Hunke and Dukowicz (1997) added a simplified elastic component to the viscous–plastic approach. This elastic part is not intended to represent the physical elastic strains but is introduced to allow an explicit numerical scheme with a large sub–time step. This explicit numerical scheme is therefore suitable for parallelization. In the cavitating fluid approach (Flato and Hibler 1992), sea ice is assumed to have no shear strength. This way of representing ice interactions is mathematically simple and easy to implement. However, the lack of shear strength induces too-large sea ice drifts and led some authors to question the validity of this model (Steele et al. 1997; Arbetter et al. 1999). Flato and Hibler also described how shear stresses can be added to the cavitating fluid model using a Mohr–Coulomb failure criterion. As an extension to the work of Flato and Hibler (1992), Tremblay and Mysak (1997) developed a viscous–plastic sea ice model based on the physics of a slowly deforming granular material. The failure in shear is represented by a Mohr–Coulomb failure criterion. Reorganization of the granules (ice floes) during shear deformation can lead to densification or dilatation of the material. This phenomenon is referred to as the dilatancy effect. With the continuum approach, the microphysical dilatancy effect is taken into account by allowing convergence (densification) or divergence (dilatation) during shear deformation. Hence, the dilatancy effect modifies the flow rule of a yield curve using a Mohr–Coulomb failure criterion. In this model, the macroscopic angle of friction is expressed in terms of the angle of dilatancy and the microphysical angle of friction. This clearly indicates the link between the deformations and the resistance to deformations.

There is currently a tendency among the different global climate modeling groups to adopt the Hibler (1979) rheology or modified versions thereof. Of the 23 models participating in the Intergovernmental Panel on Climate Change (IPCC) fourth assessment, 4 have no dynamics or a very crude representation of the rheology, 4 use the cavitating fluid model (i.e., without shear strength), and the remaining 15 models use the Hibler (1979) model or modified versions (e.g., Hunke and Dukowicz 1997). Many authors have shown that the choice of yield curve and flow rule has a significant effect on the simulated thickness and velocity fields (Ip et al. 1991; Arbetter et al. 1999; Zhang and Rothrock 2005). The use of different rheologies would raise different questions and further our understanding of the role of sea ice rheology on the sea ice drift and its impact on the climate system in general.

The main purpose of this work is to present a spherical coordinate version of the granular sea ice model (GRAN) of Tremblay and Mysak (1997). The yield curve in this model is based on a Mohr–Coulomb failure criterion as compared to an ellipse in the Hibler (1979) model. Recent laboratory measurements and field experiments suggest the use of a Coulombic yield curve (e.g., Overland et al. 1998; Hibler and Schulson 2000). In GRAN, a nonnormal flow rule with the dilatancy effect is used.

Starting from the Cartesian coordinate version of the GRAN (Tremblay and Mysak 1997), we derive the governing equations in spherical coordinates. Due to the curvature of the grid, extra terms^{1} arise in the formulation of the divergence of the stress tensor and of the strain rates in spherical coordinates owing to the differentiation of the scale factors. The metric terms are important near the poles, and can have a significant contribution to the simulated thickness and velocity fields (Hunke and Dukowicz 2002). The metric terms are incorporated into the model presented here.

As an application, the spherical coordinates version of GRAN is coupled to the Earth System Climate Model of the University of Victoria (UVic ESCM; Weaver et al. 2001). Currently, this model includes the elastic–viscous–plastic (EVP) dynamics of Hunke and Dukowicz (1997). However, it is important to note that the GRAN presented in this paper could be coupled to any global or regional atmosphere–ocean model. The GRAN will be part of a future release of the UVic ESCM (M. Eby 2006, personal communication). This will provide a platform to perform climate simulations using two different rheologies.

Various aspects of the GRAN simulations of sea ice in the Arctic are validated with and compared to the standard datasets of recent sea ice observations such as ice thickness distribution, sea ice drift, and Fram Strait ice volume and ice area fluxes. The dynamics of sea ice is strongly coupled to its thermodynamics through the sea ice thickness and concentration-dependent compressive and tensile strengths. Accordingly, we also validate the thermodynamical part of the model against the year-long datasets from the Surface Heat Budget of the Arctic Ocean (SHEBA) field experiment (see Uttal et al. 2002 for a summary). While these datasets are for single point measurements only, they are useful verification tools due to their high temporal resolution.

This paper is structured as follows. In section 2, we give a general overview of the GRAN and describe the UVic ESCM. The transformation of the sea ice model equations to spherical coordinates and the coupling of the GRAN to the UVic ESCM are presented in section 3. At the end of this section we compare the convergence and the numerical efficiency of the GRAN with that of the EVP. Various simulation results of the new coupled model and their comparison with data and simulations of the EVP model are shown in the first three parts of section 4. In the fourth part of section 4, the simulation results for the sea ice seasonal heat budget are compared with data from the SHEBA campaign. The main conclusions are summarized in section 5.

## 2. Models

### a. The granular sea ice model

In GRAN both the advection and the acceleration terms in the momentum equation are neglected. These assumptions are valid for length scales larger than a couple of kilometers and for forcing time scales of about half a day and longer. The remaining terms in the momentum equation (which is now diagnostic) are the Coriolis force, wind stress, water drag, divergence of the stress tensor, and sea surface tilt, namely,

where *ρ _{i}* is the density of the ice,

*h*the sea ice thickness,

*f*the Coriolis parameter,

**u**

*the sea ice velocity,*

_{i}*A*the sea ice concentration,

*τ*the wind stress,

_{a}*τ*the water drag,

_{w}*σ*the internal ice stress tensor,

*g*the gravity, and

*H*is the sea surface height.

_{d}In GRAN, sea ice is considered to be a slowly deforming granular material (Balendran and Nemmat-Nasser 1993). The continuum approach is adopted; individual ice floes are not resolved but the microphysical dilatancy effect is represented in the macrophysical constitutive equation. Hence, the present formulation allows for convergence (densification) or divergence (dilatation) of the pack ice during shear deformation. The dilatancy effect is determined by *δ*, the angle of dilatancy, which is the angle between the microscopic plane of motion of the individual floes and the macroscopic sliding plane (see Tremblay and Mysak 1997, their Fig. 6). A positive angle of dilatancy leads to dilatation, while a negative *δ* characterizes densification. Observations show that, on average, shear deformations are associated with dilatation (Stern et al. 1995).

The yield curve is represented in stress invariant space by a triangle (Fig. 1). The two stress invariants are the pressure *p*, which is the (negative) average of the normal stresses, and *q*, the maximum shear stress at a point. The triangular yield curve is based on a Mohr–Coulomb failure criterion for which the maximum allowable *q* (*q*_{max}) is proportional to the ice pressure *p*. The Mohr–Coulomb failure criterion is defined by the following equation:

where *ϕ _{f}* is the angle of friction and

*q*

_{max}defines the largest shear stress that sea ice can sustain as a function of pressure. When

*q*=

*q*

_{max}, plastic deformation occurs, consisting of shearing along sliding lines plus some convergence/divergence determined by

*δ*and the rate of shearing.

Failure in compression occurs when *p* reaches the maximum allowable pressure denoted by *P*_{max}. This pressure is parameterized, according to Hibler (1979), by the following equation:

where *P** is the ice strength per meter, *A* is the ice concentration, and *C* is the ice concentration parameter, an empirical constant characterizing the dependence of the compressive strength of sea ice on the ice concentration. When the pressure reaches *P*_{max}, sea ice can no longer sustain the compressive load; in this case, convergence (ridging) occurs.

GRAN also allows for cohesion or tensile stresses. Failure in tension occurs when *p* reaches the tensile strength *P*_{min}. In this case, the ice drifts freely.

To close the system of equations, a flow rule must be specified. We assume that the plane of maximum (minimum) normal strain rate is aligned with the plane of maximum (minimum) normal stress, that is, the principal strain rates are aligned with the principal stresses. In Flato and Hibler (1992), the sea ice pressure is calculated assuming that the flow is nondivergent when the pressure is between *P*_{max} and *P*_{min}. In GRAN, the divergence (*ε̇*_{kk}) is related to the angle of dilatancy and the rate of shearing according to the equation

where *γ̇* is the rate of shearing.

From these considerations and upon neglecting the much smaller elastic deformations (Coon et al. 1974), the internal ice stresses can be written as

where *δ _{ij}* is the Kronecker delta and

*ε̇*

_{ij}are the strain rates. See Tremblay and Mysak (1997) for a detailed derivation of Eq. (5). The shear viscosity

*η*is given by

When the strain rates are small, the denominator in Eq. (6) tends to zero and the shear viscosity is capped at *η*_{max}. When the shear viscosity is equal to *η*_{max}, the state of stress lies inside the yield curve and sea ice behaves as a very viscous fluid. When the shear viscosity is smaller than *η*_{max}, the state of stress is on one of the two limbs of the triangle and plastic deformations occur (shear plus divergence/convergence).

### b. The UVic ESCM

The UVic ESCM belongs to the family of Earth System Models of Intermediate Complexity (EMIC) (see Claussen et al. 2002). This means that some components of the earth system are described in a less sophisticated way than in coupled general circulation models (GCMs); that is, they include fewer processes or use a coarser resolution. The reduced complexity allows transient model runs over long periods with reasonable computational time. The first release of the UVic ESCM is described in Weaver et al. (2001). The UVic ESCM used in this study is version 2.6, which consists of oceanic, atmospheric, and sea ice components. The resolution of the model is 3.6° zonally and 1.8° meridionally, and the time step used is 1.25 days for the ocean and 0.625 days for the atmosphere; the coupling of the two components is done every 2.5 days.

The ocean model is the 3D primitive equation GCM Modular Ocean Model version 2.2 (MOM 2.2; Pacanowski 1995) with 19 unequally spaced layers in the vertical. The near-surface layer is 50 m thick and the bottom layer is 518 m thick. The atmospheric model is an energy moisture balance model (EMBM) (Fanning and Weaver 1996) with heat and moisture transport parameterized by a diffusion process. In this model, the moisture can also be advected. The wind stress forcing and the atmospheric CO_{2} concentration are prescribed. The sea ice component is a dynamic–thermodynamic model. The thermodynamics is calculated with a zero-layer Semtner model (Bitz et al. 2001) and the original dynamical part of the UVic ESCM is the EVP model (Hunke and Dukowicz 1997). The ice advection is calculated with a first-order upstream advection scheme.

In the UVic ESCM, the computational grid can be rotated to avoid the singularity at the poles due to the convergence of the meridians of longitude. In the rotated system, the North Pole is in Greenland and the South Pole is over East Antarctica. Iceland is removed to allow for a faster North Atlantic drift, and the Bering Strait is closed owing to the coarse resolution.

## 3. Transformation and coupling

### a. Transformation

The steady-state sea ice momentum equation without advection contains two terms with derivatives: the rheology and the sea surface tilt terms [see Eq. (1)]. The latter term can be rewritten using the geostrophic approximation for the ocean currents. The divergence of the stress tensor can be split into two parts: the gradient of the pressure *p* and the divergence of the deviatoric stress tensor (**∇** · *σ*′). The gradient of the pressure in spherical coordinates is given by

where *λ* is the longitude, *ϕ* is the latitude, *R _{e}* is the radius of the earth, and

**and**

*λ***are unit normal vectors. The divergence of a two-dimensional tensor in spherical coordinates can be found in Schade (1997). Using Eq. (5), the divergence of the deviatoric stress tensor can be written as follows:**

*ϕ*where the *ε̇*_{ij} (*i*, *j* = 1, 2) are the strain rates. The top line in Eq. (8) is the zonal component and the bottom line is the meridional component. Owing to the curvature of the grid, metric terms arise in the right-hand side of Eq. (8) upon carrying out the differentiation with respect to *ϕ*.

The strain rates *ε̇*_{11}, *ε̇*_{22}, and *ε̇*_{12} in spherical coordinates are given by

The strain rate term is symmetric, that is, *ε̇*_{12} = *ε̇*_{21} (Tremblay and Mysak 1997). In spherical coordinates, the strain rates also include metric terms [i.e., the second term on the right-hand side of Eqs. (9) and (11)].

### b. Coupling

The GRAN sea ice dynamics is coupled to the atmospheric component of the UVic ESCM and uses the same resolution and time step as the latter. The thermodynamics and advection are calculated with the original UVic ESCM model. The dynamical part can now be calculated using the original EVP sea ice dynamics of the UVic ESCM or the GRAN model.

The numerical scheme of the latter is the same as in the Cartesian version of the model (Tremblay and Mysak 1997). Step 1 calculates the free-drift solution as a first approximation for the velocity field. In step 2, the pressure term is added and the initial velocities and pressure (assumed equal to zero in step 1) are corrected. During this correction, the off-diagonal terms (such as, e.g., the Coriolis term) are considered constant. In step 3, the shear viscosities are calculated from the strain rates and ice pressure calculated in step 2. This linearizes the momentum equation. Finally, in step 4, the velocities are updated taking into account the complete rheology; however, the shear viscosities are kept constant.

Because the momentum equation is linearized, steps 2, 3, and 4 are repeated several times. We refer to this as the “superloop” iteration. The number of superloops is analogous to the number of pseudo time steps in Zhang and Rothrock (2000) and the sub–time step in Hunke and Dukowicz (1997). To ensure sufficient convergence of the velocity field, the spherical version of GRAN needs at least six superloop iterations. With sufficient convergence we mean that the total kinetic energy (KE) of the Arctic sea ice is within 10% of the final value if an infinite number of superloops were used.

The pressure *p* in step 2 is calculated in a similar way as described in Flato and Hibler (1992). A Reynolds decomposition for the velocity and pressure fields is inserted into the momentum equation and the basic state is subtracted. The shear stresses, the deviatoric stresses, and the off-diagonal terms are assumed to be constant. The pressure perturbation is then calculated as a function of the divergence of the sea ice velocity field.

The free-drift velocities are computed on an Arakawa B grid, whereas the velocity updates are computed on an Arakawa C grid. At the end of the dynamics calculation the ice velocities are interpolated again onto an Arakawa B grid and passed to the UVic ESCM. The other variables such as *p*, *ε̇*, and *η* are calculated at the tracer point, in the center of the grid cell.

The analysis of the numerical efficiency and convergence depends highly on the convergence criteria. Inspired by Zhang and Rothrock (2000), we use the total Arctic KE as a measure of the convergence. A simulation of 50 days is performed. We consider the GRAN to have converged when the KE has reach 10% of its final value. Similarly, the necessary length of the sub–time steps for the EVP is determined. Using this criterion, GRAN needs six superloops and the EVP sub–time step has to be approximately 90 s, on average. The corresponding mean CPU times are of similar magnitude with 2400 ms per time step for the GRAN and 2200 ms per time step for the EVP.

## 4. Simulation results

The spinup of the model was done in two phases. We started from an existing equilibrium run provided by the UVic ESCM with an atmospheric CO_{2} concentration of 280 ppm. The first 200 years of spinup was done with a constant CO_{2} concentration. The second spinup phase between 1800 and 1948 was done with a linearly increasing CO_{2} concentration going from 280 to 316 ppm in 1959. During this phase the daily wind stress fields are specified from a random permutation of the years from 1950 to 2000 of the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis. This was done in order to reduce any biases due to different modes in the atmospheric circulation. The actual run from 1948 onward was done with a CO_{2} concentration following the Mauna Loa record (Keeling and Whorf 2005) and the actual daily reanalysis wind stress fields from NCEP–NCAR.

In Table 1 the physical parameters and constants used for this simulation are listed. The diagnostic variables calculated by the sea ice dynamics model are the sea ice velocities and pressure. The sea ice thermodynamics and sea ice advection are calculated with the original UVic model. In the following three sections (4a–c), we compare the GRAN output with some standard observations and with an EVP run. In section 4b, the simulation results with the GRAN are compared with data from the SHEBA campaign.

### a. Fram Strait sea ice export

The major export site for sea ice into the North Atlantic is Fram Strait. The freshwater flux associated with this export can affect the meridional overturning circulation (Mauritzen and Häkkinen 1997; Holland et al. 2001; Mysak et al. 2005) and hence plays an important role in determining the local and global climate. This flux through Fram Strait involves the interplay between the thermodynamics (ice thickness) and the dynamics (velocity) of sea ice and, hence, its simulation provides a measure of the overall performance of the sea ice model.

Vinje et al. (1998) computed the ice volume flux through Fram Strait from August 1990 to July 1996 at around 79°N. A velocity profile across Fram Strait was determined from buoys and satellite measurements. To get a continuous velocity record, the drift was related to the pressure gradient across Fram Strait and a climatological ocean velocity. The ice thickness profile was obtained from upward-looking sonars deployed in Fram Strait. From observations at four draft measurement stations between September 1992 and March 1993, a relationship was derived in order to get a thickness profile from the single-point measurement time series from the other years. Kwok and Rothrock (1999) provide an ice area flux record derived from satellite measurements at 81°N. The ice motion record spans the period from October 1979 up to May 1996. The summer velocities obtained from the satellite observations are less reliable; thus these velocities are parameterized in terms of the sea level pressure gradient. In addition, Kwok and Rothrock (1999) used the same thickness profiles as in Vinje et al. to compute the ice volume flux through Fram Strait at 79°N from October 1990 to May 1996.

Figure 2 shows the simulated ice volume flux at 79.2°N (solid), along with the estimates from Vinje et al. (1998) (dashed), Kwok and Rothrock (1999) (dash–dotted), and the EVP simulation (dotted). The GRAN flux is in good agreement with the estimations. The correlation coefficient between the simulated ice volume flux and the Vinje et al. data is 0.58. For the Kwok and Rothrock estimations the correlation coefficient is 0.65. The EVP has a correlation coefficient of 0.60 with Vinje et al.’s estimates and 0.70 with Kwok and Rothrock’s data. The EVP and GRAN curves are remarkably similar, and indeed have a correlation coefficient of 0.97. All correlation coefficients are well above the 99% significance level. Köberle and Gerdes (2003), who employed a higher-resolution Hibler-type sea ice model forced by NCEP–NCAR reanalysis data, obtained a correlation coefficient of 0.68 with Vinje et al.’s estimates.

The winter (October–May) sea ice area flux is also compared with that of Kwok and Rothrock (1999). These fluxes are considered more reliable as they are based on satellite measurements only. The amplitude of the simulated total area flux is generally in good agreement with the observations from Kwok and Rothrock (1999; Fig. 3). Owing to the model’s land–ocean configuration, the simulated area flux is calculated at 79.2°N in comparison with 81°N for the observations. The simulated mean winter area of ice exported over the whole period is 574 × 10^{3} km^{2}, which is about 15% smaller than the observed mean area of 669 × 10^{3} km^{2}. The correlation coefficient between the two time series is 0.79, significant at the 99% level. The correlation with the EVP run is 0.76. Köberle and Gerdes (2003) reported a correlation coefficient of 0.84 over the same period, and Hilmer and Jung (2000) obtained a correlation coefficient of 0.88.

### b. Sea ice thickness

In 1991 and 1995 the *European Remote Sensing Satellites ERS-1* and *ERS-2* were launched. Laxon et al. (2003) used the measurements from the radar altimeter to estimate the ice freeboard from 1993 to 2001 in the latitudinal band from 65° to 81.5°N. The ice freeboard, which is the thickness above sea level, was converted to ice thickness using fixed ice and seawater densities in addition to a climatological snow cover. The resulting uncertainty is of the order of 40 cm for the mean thickness. Their altimeter estimates excluded thicknesses below 1 m.

Figure 4 compares the simulated thickness distribution (Fig. 4a) with the winter thickness climatology for 1993–2001 derived from satellite observations (Fig. 4b). The winter here is defined as the months between October and March, inclusive. Qualitatively, the model agrees well with the observations. The thickest simulated ice of 3–5 m is north of Greenland and north of the Canadian Archipelago. In these regions the observed ice thickness ranges from 3.5 to 4.5 m. Overall, the thickness distribution in the Arctic Ocean is well captured. Further details on the simulated seasonal ice thickness in the Beaufort Sea will be given in section 4d.

Figure 4a also reveals that the simulated ice edge extends well into the Norwegian and Barents Seas. This is a common problem in coarse-resolution ocean models. This discrepancy occurs because the simulated North Atlantic Drift and the associated ocean heat transport are too weak. This causes a nonrealistic sea ice edge position, which is also the case for the UVic ESCM with the EVP model (Mysak et al. 2005).

### c. Sea ice velocity

The Arctic sea ice atlas from the Environmental Working Group (EWG) provides sea ice motion data. The data were compiled from measurements taken from the International Arctic Buoy Program (IABP), the Russian Drifting Automatic Radiometeorological Stations (DARMS) and manned research stations. In this dataset, the climatological ice motion was computed for the years 1950 up to 1996. Then, the monthly mean velocity distribution was calculated as a deviation from the climatological field. For each particular month, a certain region around an actual trajectory measurement was selected. A weighting function was used to compute the deviatoric velocity around the measurement site, whereas for the remaining velocity field the climatological ice motion was used. The motion data are available on a CD-ROM (Arctic Climatology Project 2000).

To compare the simulated velocities with the EWG ice motion data we chose two periods. The mean velocities from December to March were first calculated for the years 1984–88. Then, the winter velocities, as above, were calculated for the period 1989–93. These two periods correspond to, respectively, the anticyclonic (1984–88) and cyclonic (1989–93) regimes described in Proshutinsky and Johnson (1997). The anticyclonic ice motion in the EWG data (Fig. 5b) shows a large Beaufort gyre (BG). The observed velocities are approximately 1 cm s^{−1} in the BG. The Transpolar Drift Stream (TDS) has a fairly straight trajectory from the Laptev and East Siberian Seas with a velocity of about 2 cm s^{−1}. The simulated anticyclonic regime (Fig. 5a) has a BG that is shifted too far west. However, the location of the BG agrees with the results from Proshutinsky and Johnson (1997). The simulated velocities are about 1 cm s^{−1} in the BG, which compare favorably with the observations. The simulated velocities north of Alaska also agree with the observations. The TDS starts mainly in the Laptev Sea and its speed is slightly higher than the one observed. As mentioned above, EWG data used climatological data if there were no observations for a certain period. Therefore, the average velocities were calculated only in the interior of the Arctic basin, where there are sufficient buoy trajectories. The EWG mean velocity is 1.9 cm s^{−1} and the simulated mean velocity is 2.1 cm s^{−1}. This difference is due to a faster simulated TDS.

The observed cyclonic regime (Fig. 6b) shows a smaller BG as compared to the anticyclonic regime, as described in Proshutinsky and Johnson (1997), and the center is closer to the northern Canadian coastline. The velocities are again about 1 cm s^{−1} in the BG. The TDS starts off moving from the Siberian coast toward the BG and then turns toward Fram Strait, producing the cyclonic motion. The observed velocities of the TDS are higher than in the anticyclonic regime with magnitudes around 3 cm s^{−1}. The simulated velocities for the cyclonic regime are shown in Fig. 6a. The BG is smaller but again shifted to the west. However, the simulated location also agrees with the results from Proshutinsky and Johnson (1997). The simulated magnitude of the ice motion is comparable to the observed speeds. The simulated TDS is similar to the observed trajectories from the Laptev Sea toward Fram Strait. Contrary to the observations, however, the simulated velocities of the TDS are lower in the cyclonic regime than in the anticyclonic regime. The average speed in the interior of the Arctic basin is 2.0 cm s^{−1} for the EWG data and 1.7 cm s^{−1} for the simulated ice motion. Again, the discrepancy is caused by a different TDS speed. Similar findings for the cyclonic regime were obtained by Mysak et al. (2005) with the UVic ESCM and using the EVP rheology.

### d. Comparison with SHEBA data

During the SHEBA field experiment, the Canadian icebreaker *Des Groseilliers* was frozen in the Beaufort Sea from October 1997 until October 1998. One of the observational programs of the field campaign included measuring the mass balance of sea ice (Perovich et al. 2003).

From autumn 1997 to autumn 1998 the SHEBA ice camp drifted from 75°N, 142°W to 80°N, 166°W. To compare the simulated data with the observations, we determined the time and position of the field station in relation to the model grid. Although the dataset is a point measurement, its temporal resolution is high, providing a good opportunity to validate the thermodynamic component of a sea ice model. A detailed study on the energy balance of the SHEBA dataset for model validation purposes is given by Huwald et al. (2005a,b).

#### 1) Sea ice thickness

The Pittsburgh site was located on a multiyear ice floe. We chose data from gauges 53, 69, and 71 to compare with our simulation. The data from these gauges were used by Huwald et al. (2005a) to produce a “best guess” ice and snow thickness time series. The two time series were adopted by the Sea Ice Model Intercomparison Project, Phase 2 (SIMIP2; violet line in Fig. 7 is the “best guess” thickness series). In the following we compare the simulated ice thickness evolution against the best-guess estimate of Huwald et al. The individual gauge time series provide an envelope of temporal evolutions of the ice thicknesses at the Pittsburgh site.

In the UVic ESCM, four processes contribute to seasonal ice thickness changes: 1) the melt/growth associated with the ocean and net atmospheric heat fluxes, 2) ice growth over open water, 3) sublimation/deposition, and 4) horizontal ice advection (i.e., ridging and lead opening). During the SHEBA campaign, the growth over open water is small in this region since the open water fraction is nearly zero during this period (Fig. 8). Furthermore, the sublimation/deposition is small. The simulated contribution to the thickness change due to advection is also small (Fig. 9). Since the SHEBA data are point measurements where ice only grows/melts thermodynamically and processes 2, 3, and 4 stated above are small, we examine the growth/melt of the sea ice due to thermodynamic processes only.

The initial simulated thickness in October 1997 (around Julian day 300) is 35 cm smaller than the observed best-guess initial thickness (Fig. 7). The simulated growth rate, however, is larger than that observed (as seen in the individual gauges and the best-guess estimate). As a consequence, the simulated maximum thickness is about 30 cm greater than the best-guess-observed thickness at around Julian day 500. Note that the sudden rises and falls in the simulated thickness are mainly due to the position change of the SHEBA camp from one model grid cell to another. The onset of the melting phase starting in June 1998 (around Julian day 540) is close to that observed. The simulated melt rate is 25% lower than the best-guess melt rate. Huwald et al. (2005b) estimated that about 60%–70% of the total melt at the Pittsburgh site was caused by surface melting. In the UVic ESCM the surface and ocean heat budgets contribute equally to the melt.

#### 2) Surface energy budget

The growth season can be separated into two parts: 1) the polar night, when the incident shortwave radiation is essentially zero (Julian days 300–410), and 2) the polar day, when shortwave radiation is present (Julian days 410–540). During the polar night, the net flux at the surface compares well with the observed mean flux (Fig. 10f). However, the simulated net flux does not exhibit the observed synoptic-scale variability associated with frontal activity, changes in cloud cover, and hence the net longwave radiation flux. The higher sea ice growth rate in early winter can be explained in part by the thinner sea ice, which grows faster, and the relatively low ocean heat flux (see Fig. 11). In addition, the ice thermodynamic model does not consider internal heat storage. This leads to a higher amplitude of the seasonal cycle and a faster growth rate in winter (Bitz and Lipscomb 1999). In the real world, the heat lost at the surface is used for brine pocket freeze up, which reduces the growth rate. Holland et al. (1993) observed a decrease in mean ice thickness with the inclusion of internal heat storage.

After the onset of the polar day (Julian day 410), the surface ice temperature is higher than observed (not shown), owing to the lower albedo value used in the UVic ESCM (Fig. 10a) and the fact that the thermodynamic model has no internal heat storage. As a consequence, the net heat flux emitted from the surface is higher than the observed flux (Fig. 10f). This also causes a much higher sensible heat flux after Julian day 450 when the downwelling shortwave radiation starts to become important (Fig. 10d). In spring, the air temperature begins to rise, but the northward transport of moisture in the model is smaller than in reality (not shown). This causes a larger simulated latent heat flux (Fig. 10e) as well. Note that the specific humidity in the model represents an average over the lower 1.8 km of the atmosphere (Weaver et al. 2001), whereas at SHEBA the specific humidity of the atmosphere is taken at 10 m above ground (Persson et al. 2002). In addition, the heat transfer coefficient in the model, which is a function of wind speed and temperature, is generally twice as large as the values typically used by other modeling groups. This issue has been resolved in version 2.8 of the UVic ESCM (M. Eby 2006, personal communication).

In the model, the ice layer starts to melt at Julian day 540 (25 June) (Fig. 7), a few days after the insolation reached its maximum value. From then on the incoming shortwave radiation is decreasing. This, together with the constant albedo value in the UVic ESCM, decreases the net shortwave radiation after the onset of melt (Fig. 10b). The sensible heat flux drops to almost zero 20 days later (Fig. 10d) as both the surface air and ice layer are around 0°C. However, the net longwave radiation is overestimated by about 20 W m^{−2} and the latent heat by about 35 W m^{−2} (Figs. 10c and 10e). The overly large simulated net longwave radiation is mainly caused by a low simulated downward longwave radiation. The higher latent heat loss is due to a continued drier atmosphere in the model. The maximum observed magnitude of the net surface heat flux is about 75 W m^{−2} as compared to about 30 W m^{−2} in the model (Fig. 10f).

#### 3) Basal energy budget

The observed ocean heat flux shows a large peak around Julian day 440 (Fig. 11). This peak is associated with the ice station drifting into shallow waters during a storm event and entraining warmer, deeper water (Perovich and Elder 2002). A second large change in measured ocean heat flux starting around Julian day 500 is again associated with a large storm event. Perovich and Elder (2002) calculated a mean ocean heat flux of 7.5 W m^{−2} over the entire period at the Pittsburgh site. This number is almost two to four times the values found in earlier studies (i.e., Maykut and McPhee 1995; Perovich et al. 1997). For instance, Perovich et al. (1997) found that in the Beaufort Sea the ocean heat flux in winter is close to zero and in summer it is around 5 W m^{−2}.

In the UVic ESCM the sensible heat flux of an ice-covered grid cell is calculated with the temperature difference between the ice base, set at the freezing point of water for a given salinity, and the surface ocean. In the following this sensible heat flux is referred to as the ocean heat flux. On the other hand, when the sea ice concentration is lower than 100%, heat that enters the surface waters from the atmosphere (radiation, turbulent) is used directly for basal ice melt. This is equivalent to assuming an infinite sensible heat transfer coefficient. Together with the ocean heat flux, this is referred to as the net ocean heat flux.

In the model, the ocean heat flux is on the order of 0.5 W m^{−2} (upward) during the snow/ice accumulation period (Fig. 11). After Julian day 490 the simulated ocean heat flux changes sign. At that point the snow layer starts to melt. In the model, the melted freshwater is released into the upper ocean and thus decreases its salinity. This increases the freezing point temperature (i.e., the ice base temperature), yet the actual surface ocean temperature remains nearly the same. This causes the ocean heat flux to change sign and basal ice to form (dashed line in Fig. 11). The second positive bulge of the ocean heat flux is associated with the melting of the ice layer. The positive ocean heat flux (from the ice to the ocean) in the model is associated with summer ice formation (“false bottoms”). These false bottoms were observed during the SHEBA campaign (Eicken et al. 2002; Perovich et al. 2003). The large simulated ocean heat flux pulse starting at Julian day 575 and ending at Julian day 605 is related to a sudden increase in ocean surface temperature. At this point the simulated SHEBA camp moves from shallow to deep bathymetry. The increase in temperature is presumably caused by diffusion or convection from below.

Around Julian day 520 the ice pack opens up (Fig. 8) and the sea ice concentration falls below 100%. From this date until Julian day 610 when the sea ice closes up again, heat that enters the open water is directly used for melting. The total available net ocean heat flux during the melting phase reaches a maximum of 22 W m^{−2} (solid line in Fig. 11).

## 5. Summary

The choice of rheology can lead to significant differences in the simulated sea ice thickness and velocity distributions (Ip et al. 1991; Arbetter et al. 1999; Zhang and Rothrock 2005). Which rheology is the most accurate to describe the sea ice dynamics on a large scale is still under debate. Therefore, it is useful to test different rheologies in global climate simulations. This can further our understanding of the role of sea ice rheology and its impact on the climate system and help raise new questions. Thus, in this paper we presented the transformation to spherical coordinates of the GRAN model (Tremblay and Mysak 1997). The GRAN is coupled to the UVic ESCM as an application and it will be part of a future release of this model (M. Eby 2006, personal communication). Its output is compared to different observational datasets and an EVP model simulation.

The simulated volume export through Fram Strait is in good agreement with the estimates from Vinje et al. (1998) and Kwok and Rothrock (1999). The comparison with the ice area flux estimated by Kwok and Rothrock (1999) shows that the area flux follows the estimated behavior, although the mean is 15% lower. The winter climatological thickness distribution corresponds well with the data provided by Laxon et al. (2003). The sea ice thickness distribution in the central Arctic Ocean is captured. The cyclonic and anticyclonic ice drift regimes found by Proshutinsky and Johnson (1997) are well captured and the simulated velocities are comparable to the buoy velocities (Arctic Climatology Project 2000). However, we note that the center of the simulated BG is located too far west in both regimes compared to the observations. Nonetheless, the location of the center in each regime agrees with that shown in Proshutinsky and Johnson (1997). The simulated mean velocity of the sea ice in the central Arctic is slightly too large for the anticyclonic regime and slightly too small in the cyclonic regime.

Finally, the thickness evolution in the Beaufort Sea was compared against the SHEBA datasets. The net growth of the ice layer compared to SHEBA data is more than twice as large. This is mainly caused by the absence of heat storage in the ice layer, the lower ocean heat flux, and the thinner simulated ice at the beginning of the SHEBA observation period. However, the simulated onset of melting coincides with the observed onset. The total melt in late summer is about 25% too low compared with SHEBA.

During the growing phase in winter, when the sun is below the horizon, the net surface heat flux is comparable with the SHEBA data. However, owing to the absence of internal heat storage in the ice layer and the initially thin ice, the growth rate is too large. During the polar day, the ice layer absorbs too much shortwave radiation due to the lower albedo value. Again, owing to the missing internal heat storage, the surface ice temperature increases too quickly, causing high sensible and latent heat fluxes. For the latent heat, this is exacerbated by the low atmospheric relative humidity simulated by the model. During the melting phase the net surface heat balance is lower than observed.

The simulated ocean heat flux is low during the growing season. At the onset of the snowmelt the model simulates the formation of false bottoms, similar to what was observed during SHEBA (Eicken et al. 2002; Perovich and Elder 2002). After breakup of the pack ice, an additional bottom melt effect appears in order to account for the heating of the mixed layer. In the model, the surface and bottom melt contribute equal amounts to the total simulated melt, but it is too low compared to observations.

## Acknowledgments

We thank M. Eby for his considerable help during the coupling of the GRAN to the UVic ESCM and H. Huwald for providing the skin temperatures at the SHEBA site. Furthermore, we thank the two anonymous reviewers for their comments and suggestions, which helped to improve the paper. This work was supported by research grants awarded to L. A. Mysak from NSERC, the Canadian Institute of Climate Studies (Arctic Node), and the Canadian Foundation for Climate and Atmospheric Sciences (Canadian CLIVAR Research Network). L. B. Tremblay was supported by an NSERC Discovery Grant and the National Science Foundation Grants OPP-0230264, OPP-1230325, and ARC-05-20496.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**.**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**.**

**,**

_{2}records from sites in the SIO air sampling network. Trends: A Compendium of Data on Global Change, Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**.**

**.**

**,**

**,**

**.**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**.**

## Footnotes

*Corresponding author address:* Jan Sedlacek, Earth System Modelling Group, Department of Atmospheric and Oceanic Sciences, McGill University, Montreal, QC H3A 2K6, Canada. Email: jan.sedlacek@mail.mcgill.ca

^{1}

Sometimes they are referred to as metric terms (e.g., Hunke and Dukowicz 2002).