Abstract

A stand-alone sea ice model is tuned and validated using satellite-derived, basinwide observations of sea ice thickness, extent, and velocity from the years 1993 to 2001. This is the first time that basin-scale measurements of sea ice thickness have been used for this purpose. The model is based on the CICE sea ice model code developed at the Los Alamos National Laboratory, with some minor modifications, and forcing consists of 40-yr ECMWF Re-Analysis (ERA-40) and Polar Exchange at the Sea Surface (POLES) data. Three parameters are varied in the tuning process: Ca, the air–ice drag coefficient; P*, the ice strength parameter; and α, the broadband albedo of cold bare ice, with the aim being to determine the subset of this three-dimensional parameter space that gives the best simultaneous agreement with observations with this forcing set. It is found that observations of sea ice extent and velocity alone are not sufficient to unambiguously tune the model, and that sea ice thickness measurements are necessary to locate a unique subset of parameter space in which simultaneous agreement is achieved with all three observational datasets.

1. Introduction

Sea ice is important to the climate system because it insulates the cold polar atmosphere from the warmer ocean in winter, reflects a high proportion of incoming shortwave radiation in summer, and both stores and transports freshwater and negative latent heat. Brine released during ice formation is also thought to play a role in the formation of oceanic deep water.

It is therefore important that sea ice is modeled as accurately as possible in the global circulation models (GCMs) that are used to make predictions of how the earth's climate will respond to various greenhouse gas emission scenarios (Houghton et al. 2001). Possible future changes in the ice–albedo feedback due to a changing Arctic sea ice cover represent a large uncertainty in the prediction of future temperature rise (Carson 1999; Rind et al. 1997).

However, uncertainty still surrounds some aspects of contemporary sea ice models. This was recently illustrated by Rothrock et al. (2003), who compared predictions of Arctic ice thickness from various published studies to highlight the substantial disagreement among the models used (their Fig. 12). They conclude that the wide range of thickness predictions must have been a result of differences in the representation of sea ice physics and/or the forcing data used. Unfortunately, no single clear reason for the differences emerged. To increase confidence in simulations of the ice cover, progress must therefore be made in assessing the quality of forcing data (Fischer and Lemke 1994; Curry et al. 2002) and in improving the representation of physical processes in sea ice models.

In recent decades, considerable progress has been made in improving the representation of sea ice in models (Hibler 2004). An evolution equation for the ice thickness distribution (ITD) was introduced by Thorndike et al. (1975), and the widely used viscous–plastic (VP) rheology was introduced by Hibler (1979). Bitz and Lipscomb (1999) introduced an energy-conserving sea ice model that explicitly accounts for the effect of internal brine pocket melting on surface melt, and observations made during the yearlong SHEBA experiment (Perovich et al. 1999, 2003) have resulted in better parameterizations of surface albedo (Curry et al. 2001; Perovich et al. 2002) as well as more detailed assessments of previously used parameterizations (Eicken et al. 2004).

Notwithstanding these developments, even models with identical parameterizations of sea ice physical processes often differ in the values of parameters they use. For example, a quadratic air drag parameterization (McPhee 1975) [see Eq. (4) below] is used in many models, but the air drag coefficient, Ca, often differs. For geostrophic wind forcing, Hibler (1979) and Lindsay and Stern (2004) both use Ca = 0.0012, and Rothrock et al. (2003) vary Ca sinusoidally between 0.0006 on 1 January and 0.0012 on 1 July. However, Zhang and Hunke (2001) also use Ca = 0.0012 for forcing by 10-m surface winds, which are generally weaker than geostrophic winds (Walter et al. 1984).

Similarly, Hibler (1979) introduced a widely used, empirical ice strength parameterization [see Eq. (7) below] with the strength parameter P* = 5 kN m−2, whereas Hibler and Walsh (1982) used a larger value of P* = 27.5 kN m−2, and Holland et al. (1993) used the intermediate value, P* = 10 kN m−2.

Harder and Fischer (1999) discuss the uncertainties surrounding both Ca and P* and use observations of buoy drift trajectories to tune them in their sea ice model, achieving closest agreement with observations when Ca = 0.0016 and P* = 20 kN m−2. In this paper we also address the uncertainty surrounding Ca and P* as well as the broadband albedo of cold bare ice, α, using a stand-alone sea ice model with prescribed atmospheric and oceanic forcing.

However, our study advances previous work in two respects. First, the model we use, CICE, the Los Alamos Sea Ice Model developed at the Los Alamos National Laboratory (Hunke and Lipscomb 2004), incorporates all of the recent improvements to model physics mentioned above. Second, to tune the model we now use an extended, satellite-derived, basin-scale dataset comprising sea ice extent, velocity, and, for the first time, thickness (Laxon et al. 2003).

Section 2 contains a description of the observations we use to tune the model. This is followed by a brief description of the model and forcing in section 3. In section 4 we describe the metrics that we use to quantify the agreement between the model and observations. Much of the detailed discussion of the comparison of model predictions and observations is placed in the appendix, but we summarize our findings in section 5 and describe the method used to determine the subset of parameter space that gives the best simultaneous agreement with all three sets of observations. Our conclusions are given in section 6.

2. Observations used to tune and calibrate the model

We use basin-scale, satellite-derived estimates of sea ice thickness, extent, and velocity to tune the model, which was run using 1980 to 2001 forcing data derived from the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40) (see www.ecmwf.int) and the Polar Exchange at the Sea Surface (POLES; see http://psc.apl.washington.edu/POLES/) sea ice model forcing set. This section describes these observations in detail, as well as their estimated errors.

a. Sea ice thickness

Laxon et al. (2003) describe new techniques to obtain ice thickness from satellite estimates of ice freeboard using data from the 13.8-GHz radar altimeters carried by the ERS-1 and ERS-2 satellites, which have been in orbit since July 1991 and April 1995, respectively, and provide coverage to 81.5°N. The mean winter (1993–2001) regional distribution of ice thickness determined by Laxon et al. (2003, their Fig. 1) confirms and extends previously published submarine climatologies (Bourke and McLaren 1992). Furthermore, Laxon et al. (2003) found a considerable interannual variability in winter ice thickness that is highly correlated with the length of the intervening melt season.

To validate the sea ice thickness measurements, Laxon et al. (2003, their Fig. 2) compared near-coincident satellite altimeter- and submarine-derived ice thicknesses in the Beaufort Sea in the 1990s. Taking into account the estimated measurement errors in both datasets, they found that the correlation between the altimeter and the submarine estimates was significant at the 99.9% level using a χ2 goodness-of-fit test.

Difficulties in determining the origin of the echoes under melt conditions mean that ice thickness cannot be determined between May and August. Furthermore, owing to difficulties in distinguishing thin ice (h < 0.5 m) from open water, the thickness estimates exclude areas including leads, open water, and new thin ice. We use the Laxon et al. ice thickness estimates in this study, but reduce problems with thin ice by using measurements made during the six winter months from November until April only, when thin ice and open water cover only a small percentage of the total area.

All ice thickness measurements, h, made during the 50 winter months between November 1993 and December 2001 were regridded onto the model grid (section 3c) using a Gaussian weighting scheme with a maximum search radius of 250 km. Grid cell means, hcell, of fewer than 100 weighted observations were excluded from our analysis to reduce the effects of measurement noise. This left 12 978 monthly grid cell means in this period, with a mean thickness hobs = 2.71 m and standard deviation σobs = 0.89 m.

We restricted thickness comparisons to the central Arctic ocean below the 81.5°N observational limit of the altimeter. This is due to the known inaccuracy (Hibler and Walsh 1982; Hibler and Bryan 1984) of stand-alone sea ice models in the Barents and Greenland–Iceland–Norwegian (GIN) Seas, where simulated winter ice tends to be too extensive due to the lack of lateral heat transported by the northward Norwegian and West Spitzbergen Currents (Hibler and Bryan 1987), which help to keep these seas free of ice throughout the year.

1) Sea ice thickness errors

Laxon et al. (2003) discuss in detail the errors in the method of calculating sea ice thickness using ERS altimeter data. We find that monthly grid cell means have an estimated error, Δhcell, of 4–6 cm. Errors in monthly, regional averages of ice thickness, hobst, comprising ncells grid cells, are given by

 
formula

where Δh2ρ = 0.0121 m2 is due to uncertainty in ice and water densities, and Δh2s = 0.0081 m2 is taken from the Warren et al. (1999) estimate of the interannual variability of snow depth for each month.

b. Sea ice extent

Sea ice extent, defined here as the total area with a concentration in excess of 15%, varies in the Arctic from a minimum of ∼7 × 106 km2 at the end of the summer melt season in September to a maximum of ∼15 × 106 km2 in March (Parkinson et al. 1999).

We calculate observed ice extent from 1994 to 2001 using the ice concentration data from the Special Sensor Microwave Imager (SSM/I) passive microwave radiometer (PMR), which is available from the National Snow and Ice Data Center (NSIDC: see http://nsidc.org) and has been calculated using the NASA team algorithm (Cavalieri et al. 2002). We assume concentrations in excess of 15% above the SSM/I latitudinal limit of 87.6°N.

Monthly ice concentration fields were regridded onto the model grid using a Gaussian weighting scheme with a maximum search radius of 250 km and then used to calculate the observed monthly ice extent.

1) Sea ice extent errors

Because of melt ponds and other effects, PMR algorithms underestimate summer ice concentration (Comiso and Kwok 1996; Serreze et al. 2003) and possibly sea ice extent. Partington et al. (2003) compared manual U.S. National Ice Center (NIC) ice charts to PMR sea ice concentration derived using the NASA team algorithm and found differences that peaked in summer, the PMR concentrations underestimating sea ice concentration by up to 23% in early August.

To account for the possibility that observations of summer ice concentration that lie between 0% and 15% may have actual concentrations in excess of 15%, we estimate the error in summer ice extent as the area of the grid cells in which the observed ice concentration lies between 0% and 15%.

c. Sea ice velocity

In this study we use 1994 to 2001 gridded monthly mean ice motion vector fields from the Polar Pathfinder Daily 25-km Equal-Area Scalable Earth Grid (EASE-grid) Sea Ice Motion Vectors dataset, computed by Fowler (2003) and available from the NSIDC.

Fowler used optimal interpolation to create daily gridded ice motion fields from Scanning Multichannel Microwave Radiometer (SMMR), SSM/I, Advanced Very High Resolution Radiometer (AVHRR), and International Arctic Buoy Program (IABP) buoy data. The daily gridded fields combine data from all sensors and are available for each month from November 1978 until March 2003. Monthly gridded fields were calculated from means of the daily dataset.

We regridded the monthly ice velocity fields onto the model grid using a Gaussian weighting scheme with a maximum search radius of 100 km. Grid cell means of less than ten weighted observations were excluded from our analysis to reduce the effects of measurement noise. This procedure resulted in 48 523 monthly grid cell means in the Arctic basin between 1994 to 2001, with an overall mean speed υobs = 2.8 cm s−1, and a standard deviation συ = 0.8 cm s−1.

1) Sea ice velocity errors

To gauge the accuracy of the blended dataset, Fowler (2003) interpolated several years of daily SMMR, SSM/I, and AVHRR vectors {u, υ} to the same grid but without using buoy data. The mean difference between the interpolated EASE-grid u components and the buoy vectors was 0.1 cm s−1 with a root-mean-square (rms) error of 3.36 cm s−1. For the υ component, the mean difference was 0.4 cm s−1 and the rms error 3.40 cm s−1. He also compared AVHRR and SSM/I data separately to daily buoy data and found typical errors of 3–4 and 4–5 cm s−1, respectively.

If we therefore assume an error in each daily velocity component Δdaily ≈ 5 cm s−1, then the error in each monthly velocity component will decrease as Δmonthly ≈ 5/d cm s−1, where d is the number of days in the month. However, we regridded each monthly component onto the model grid, so this further reduces the gridded errors. For each grid cell we demanded a minimum of 10 observations within 100 km, so the maximum error in monthly grid cell velocity components is Δmonthly/10 ≈ 5/300 ≈ 0.3 cm s−1, where we have assumed a 30-day month.

3. Model and forcing

For this study we have used, in stand-alone mode, version 3.0.1 of CICE, the Los Alamos Sea Ice Model, developed at the Los Alamos National Laboratory and released in 2001. Full details of the latest version (3.1) of the model can be found in Hunke and Lipscomb (2004), but we summarize here CICE's main features and describe in detail our modifications.

a. CICE dynamics

CICE characterizes the state of the ice cover in each grid cell at a time t and location x using an ITD function g(x, h, t) (Thorndike et al. 1975), where g(x, h, t) is the fractional grid cell area covered by ice in the thickness range (h, h + dh). [We abbreviate g(x, h, t) to g(h) in the following discussion.] CICE then solves Thorndike et al.'s evolution equation for g(h):

 
formula

where u is the horizontal ice velocity, f = f (h) is the thermodynamic ice growth/melt rate, and ψ = ψ(g(h)) is the ridging redistribution function, by discretizing g(h) in each grid cell into five ice thickness categories. The CICE ridging parameterization, described by ψ in Eq. (2), is identical to the formulation of Flato and Hibler (1995).

The ice velocity, u, is determined from a two-dimensional sea ice force balance per unit area (Hibler 1979):

 
formula

where m is the mass of snow and ice per unit area, σ is the ice stress tensor, k is the vertical unit vector, f is the Coriolis parameter, g is the acceleration due to gravity, and Hsea is the sea surface height. The terms on the right-hand side are the air stress, water stress, internal ice force per unit area, Coriolis force per unit area, and the force per unit area due to sea surface tilt, respectively.

The standard CICE model uses an air drag coefficient that depends on atmospheric boundary layer stability. However, to facilitate comparisons with previous studies we use a standard quadratic expression (McPhee 1975; Hibler 1979; Kreyscher et al. 2000) for air drag:

 
formula

where ua is the 10-m wind field, and Ca is the air drag coefficient that we wish to determine (see section 3d). We have assumed here that wind speeds are significantly higher than the maximum ice speeds and that the turning angle is negligible.

Drag at the sea ice–ocean interface is given by

 
formula

where uw are the POLES mean 1979–93 ocean currents at 40-m depth from the coupled ice–ocean model of Zhang et al. (1998), θ is the turning angle, ρw is the ocean density, and Cw is the ocean drag coefficient (see Table 1).

Table 1.

Fixed parameter values of the CICE sea ice model during the optimization process.

Fixed parameter values of the CICE sea ice model during the optimization process.
Fixed parameter values of the CICE sea ice model during the optimization process.

To calculate the internal ice force, the stress tensor σij is related to the deformation of the ice cover through a constitutive law or rheology. The CICE model uses the elastic–viscous–plastic (EVP) constitutive law of Hunke and Dukowicz (1997, 2002), which is based on the VP rheology introduced by Hibler (1979).

The EVP rheology relates σij to the strain rate tensor ε̇ij, an internal ice strength P, nonlinear bulk ζ(ε̇ij, P), and shear η(ε̇ij, P) viscosities:

 
formula

where E is a Young modulus. The VP rheology in Eq. (6) is recovered when E → ∞ and/or in the steady-state limit; in this case, the principal components of the stress tensor lie on an elliptical yield curve. One advantage of using the EVP over the VP rheology is that it allows larger time steps to be used (Hunke and Zhang 1999).

Central to both VP and EVP rheologies is the grid cell ice strength P. We have chosen the parameterization of Hibler (1979):

 
formula

where h is the mean ice thickness, P* (kN m−2) is the ice strength parameter that we wish to determine (see section 3d), C is a positive constant (see Table 1), and A is the ice area fraction. Thus, the ice exponentially weakens (strengthens) as the open water fraction increases (decreases) and it linearly strengthens as h increases.

CICE can also calculate sea ice strength P in its ridging parameterization as being proportional to the change in ice potential energy per unit area of compressional deformation, following Rothrock (1975) and Flato and Hibler (1995). We chose this simpler parameterization on the basis of its simplicity, but this does make it easier to compare our results to earlier model studies that used the same parameterization, such as Hibler (1979), Hibler and Walsh (1982), Holland et al. (1993), and Harder and Fischer (1999). Moreover, the parameterization is still in current use. For example, the sea ice component of the Hadley Centre's latest general circulation model (HadGEM1) uses the EVP rheology and the same, simple, strength parameterization used here (Johns et al. 2004).

b. CICE thermodynamics

To calculate the thermodynamic growth rate, f, in Eq. (2) and vertical heat transfer CICE uses the multilayer, energy-conserving thermodynamic model of Bitz and Lipscomb (1999) with four layers of ice and one layer of snow in each of the five ice thickness categories and assumes a constant, but nonlinear, vertical salinity profile. This model includes a temperature-dependent heat capacity to take into account the presence of brine pockets.

1) Albedo

CICE uses a parameterization of albedo similar to the National Center for Atmospheric Research (NCAR) Climate System Model version 1.3 (CSM1.3) parameterization described by Curry et al. (2001) and chosen to fit observations made during the Surface Heat Budget of the Arctic (SHEBA) field experiments (Hunke and Lipscomb 2004). Broadband grid-cell albedo depends on surface type (snow or bare ice), surface temperature, and the distributions of ice and snow thickness.

Two bare ice albedos are calculated for each ice thickness category N, one each for both the visible wavelength band (<700 nm), αice,NV, and the infrared wavelength band (>700 nm), αice,NIR. If the ice thickness in category N is given by hNmod, then visible albedo is first scaled with ice thickness as

 
formula

where αw is the albedo of open water given in Table 1, αV is the visible ice albedo parameter with the SHEBA-derived value, αV = 0.78, and fh is the (nondimensional) thickness scaling factor fh = min(2hNmod, 1). Similarly,

 
formula

where αIR is the infrared ice albedo parameter with the SHEBA-derived value, αIR = 0.36. Both category N bare ice albedos are then reduced linearly by up to 0.075 as each category's surface temperature TNsfc rises from −5°C to the fresh ice melting temperature of 0°C. These thickness and temperature scalings are summarized in Fig. 1.

Fig. 1.

The CICE broadband albedo, α, for bare ice. The albedo decreases when the bare ice thins below 0.5 m and/or when its surface temperature increases above −5°C and approaches the fresh ice melting temperature of 0°C.

Fig. 1.

The CICE broadband albedo, α, for bare ice. The albedo decreases when the bare ice thins below 0.5 m and/or when its surface temperature increases above −5°C and approaches the fresh ice melting temperature of 0°C.

To determine both the final visible and infrared albedos for category N, an area-weighted average of both ice and snow albedos is taken:

 
formula

and

 
formula

where αsnowV and αsnowIR are the visible and infrared snow albedo parameters with the SHEBA-derived values αsnowV = 0.98 and αsnowIR = 0.70, respectively (slightly adjusted for a temperature dependence), and fs is fractional area of snow on ice given by fs = hNsnow/(hNsnow + 0.04), where hNs is the snow thickness for category N.

We assume that 52% of the incoming (ERA-40) shortwave flux, Fsw, is in the visible wavelength band and 48% is in the near-infrared wavelength band. Thus, the final, broadband albedo for category N is given by

 
formula

Albedo values are scaled with neither temperature nor thickness for ice in a category that is cold (<−5°C) and thick (>0.5 m) and, if it does not have a snow cover, the SHEBA-derived value for the albedo is given by

 
formula

using Eqs. (8)(11). Similarly, the SHEBA-derived value for the broadband snow albedo, αs, as given in Table 1 is given by

 
formula

using Eqs. (10) and (11).

In this study we have chosen to adjust the broadband albedo for bare ice by changing one parameter, αV, from its SHEBA value, while keeping the original relationship αV = αIR + 0.42 intact. Thus, by Eq. (13), αNαV − 0.2. All other parameters are left unchanged.

2) Heat fluxes

We have replaced the standard CICE model's boundary layer stability-dependent parameterization of sensible and latent heat fluxes with bulk parameterizations. Sensible heat flux for ice thickness category N is given by

 
formula

where ρa is the air density, cma is the specific heat of moist air, Cs is the exchange coefficient (see Table 1), Tair is the 2-m air temperature, and TNsfc is the surface temperature for category N. The second term in the first factor takes into account sensible heat transfer in windless conditions (Jordan et al. 1999), and all heat fluxes are positive downward.

Similarly, latent heat flux for ice thickness category N is given by

 
formula

where Lsub is the latent heat of sublimation of freshwater, Qair is the surface specific humidity, Cl is the exchange coefficient (see Table 1), and QNsfc is the surface saturation specific humidity for category N.

We also changed the 50-m-deep oceanic mixed layer model to mimic the treatment of Ebert and Curry (1993). The mixed layer temperature, Tw, changes as shortwave radiation is exchanged between air and ocean through leads and a proportion that passes through the ice, and in response to the heat transfer into the bottom of ice due to molecular and turbulent diffusion, Fb, given by (Maykut and McPhee 1995)

 
formula

where Tf is the (salinity dependent) mixed layer freezing temperature, cw is the specific heat of seawater, u* = |τw|/ρw is a friction velocity, and Cb is an empirical heat transfer coefficient (see Table 1). We also assume that there is no vertical heat transfer from the deeper ocean. This is likely to be a good approximation in the highly stratified central Arctic Ocean, where studies performed by Hibler and Bryan (1987) and Zhang and Rothrock (2003) suggest that the average annual deep ocean heat flux is small, approximately 2 W m−2.

c. Model grid, spinup, and forcing

The model grid is a curvilinear orthogonal grid with a resolution of 1° covering the Arctic Ocean and its peripheral seas above 56°N (see Fig. 2). It was generated by rotating a latitude–longitude grid by 90°, keeping the Greenwich meridian invariant. This makes the grid almost regular over the Arctic region and avoids problems with the North Pole in global grids. In the CICE model, metric terms associated with the curvature of the grid are explicitly incorporated into the discretization of the EVP dynamics (Hunke and Lipscomb 2004; Hunke and Dukowicz 2002). The boundaries were set outside the region of maximum sea ice extent, but Bering Strait is closed off.

Fig. 2.

Model domain (black and gray regions combined). The gray region is where the Arctic basin and ERS observations overlap, and the latitudinal limit of 81.5°N is due to the orbit of the ERS-1 and ERS-2 satellites. We compare modeled and observed sea ice thickness in this gray region but include the region above 81.5°N in speed and extent comparisons, and call this larger region the “Arctic basin” in the paper.

Fig. 2.

Model domain (black and gray regions combined). The gray region is where the Arctic basin and ERS observations overlap, and the latitudinal limit of 81.5°N is due to the orbit of the ERS-1 and ERS-2 satellites. We compare modeled and observed sea ice thickness in this gray region but include the region above 81.5°N in speed and extent comparisons, and call this larger region the “Arctic basin” in the paper.

Atmospheric forcing data are derived from various sources and regridded onto the model grid. For the years 1980–2001, 6-hourly precipitation (snowfall), incoming longwave and shortwave radiative fluxes, and 10-m winds are taken from the ERA-40 reanalysis. However, daily surface (2 m) air temperatures are taken from the optimally interpolated IABP/POLES set (Rigor et al. 2000) since the ERA-40 temperatures were found to regularly exceed the freezing point over ice-covered grid cells in summer (Flato 1995; Vowinckel and Orvig 1970). An annual climatology of daily surface specific humidity from the POLES sea ice model forcing set is used to calculate latent heat fluxes.

Spacially variable but temporally constant mean 1979–93 40-m ocean currents from the coupled ice–ocean model of Zhang et al. (1998) are also taken from the POLES forcing set, and a spacially varying but temporally constant mixed layer salinity is taken from the Polar Science Center Hydrographic Climatology (PHC) climatology (Steele et al. 2001).

The ice in the model is started from rest with a uniform thickness of 1.85 m everywhere above 70°N, spun up for 11 years using repeated 1980 forcing data, and then run from 1980 until 2001 with each year's atmospheric and oceanic forcing.

d. Optimization of CICE model parameters

As discussed in section 1, considerable uncertainty surrounds the correct values of Ca [see Eq. (4)], P* [see Eq. (7)], and α (see Fig. 1) to use in sea ice models. Yet, these dynamic and thermodynamic parameters exert considerable influence on the modeled thickness, extent, and velocity (Ebert and Curry 1993; Chapman et al. 1994; Fischer and Lemke 1994; Steele et al. 1997). It is therefore important to reduce this uncertainty if we are to be confident in the predictions of similarly parameterized sea ice models.

Following the approach taken by Chapman et al. (1994), Harder and Fischer (1999), and Kreyscher et al. (2000), we tune these three parameters in the modified CICE model, forced with ERA-40 and POLES data, to locate an optimal subset of parameter space in which the discrepancies between model predictions and observations are minimized. All other parameters are fixed at commonly used values and are given in Table 1.

4. Comparison of model and observations

To illustrate the metrics that we have used to quantify the discrepancies between model predictions and observations, we have used a standard model run with parameters fixed at the values shown in Table 1 except that P* = 10 kN m−2, Ca = 0.0011, and α = 0.58 (see Figs. 3 –7) and compared its output to the observations described in section 2. These values lie roughly in the middle of the region of parameter space we searched to find the optimum parameter set (see section 5).

Fig. 3.

The observed and standard model run Arctic basin sea ice extent from 1994 to 2001 with {α, Ca, P*} = {0.58, 0.0011, 10 kN m−2}; RMSDsep = 0.25 × 106 km2.

Fig. 3.

The observed and standard model run Arctic basin sea ice extent from 1994 to 2001 with {α, Ca, P*} = {0.58, 0.0011, 10 kN m−2}; RMSDsep = 0.25 × 106 km2.

Fig. 7.

Scatterplot of the observed and the standard model run basin-averaged sea ice thicknesses for all 50 winter months between 1994 and 2001. As in Fig. 6, {α, Ca, P*} = {0.58, 0.0011, 10 kN m−2}. The dashed line is the least squares fit, and χ2thk = 175.

Fig. 7.

Scatterplot of the observed and the standard model run basin-averaged sea ice thicknesses for all 50 winter months between 1994 and 2001. As in Fig. 6, {α, Ca, P*} = {0.58, 0.0011, 10 kN m−2}. The dashed line is the least squares fit, and χ2thk = 175.

a. Ice extent comparisons

We first compared the standard model run's Arctic basin ice extent to the value calculated from SSM/I ice concentrations, described in section 2b.

In Fig. 3 we have plotted the time series of monthly Arctic basin ice extent from 1994 to 2001. Monthly ice extent has a maximum error of Δhe ≈ 0.5 × 106 km2 in September. There is also considerable interannual variability, with the lowest ice extent seen during the summers of 1995, 1998, and 1999.

Because both modeled and observed ice fill the Arctic basin region completely during winter, we have concentrated on the summer months between 1994 and 2001 in this study, and chosen as our sea ice extent metric the root-mean-square difference between modeled and observed September ice extents, RMSDsep. For the standard model run, we found RMSDsep = 0.25 × 106 km2 (see Fig. 3).

b. Ice speed comparisons

In Fig. 4 we plot the observed monthly ice speeds in the Arctic basin from 1994 to 2001. As with sea ice extent, there is considerable variability in the time series with a maximum monthly mean ice speed of ∼5.5 cm s−1 in July 1994 and a minimum of ∼1.4 cm s−1 in July 1997. Since each monthly basin-averaged speed is typically the mean of approximately 500 grid cell speeds, the error associated with points in Fig. 4 is too small to plot. Similarly, the estimated error in the overall mean, υobs, is tiny on account of the large (48 523) number of individual grid cell means available.

Fig. 4.

The observed and standard model run Arctic basin sea ice speeds from 1994 to 2001. As in Fig. 3, {α, Ca, P*} = {0.58, 0.0011, 10 kN m−2}. There is a high correlation between the time series, Rspeed = 0.81, but the modeled mean speed, υmod = 3.8 cm s−1, is considerably higher than the observed mean, υobs = 2.8 cm s−1.

Fig. 4.

The observed and standard model run Arctic basin sea ice speeds from 1994 to 2001. As in Fig. 3, {α, Ca, P*} = {0.58, 0.0011, 10 kN m−2}. There is a high correlation between the time series, Rspeed = 0.81, but the modeled mean speed, υmod = 3.8 cm s−1, is considerably higher than the observed mean, υobs = 2.8 cm s−1.

We have chosen three metrics with which to assess the model's ability to reproduce observed ice speeds. First we compared υmod, the modeled, and υobs, the observed mean of all 48 523 monthly, mean grid cell speeds between 1994 and 2001. The standard model run's mean ice speed υmod = 3.8 cm s−1 overestimates the observed mean speed, υobs = 2.8 cm s−1, by 1.0 cm s−1 or ∼36%.

For our second metric we chose to assess the standard model run's correlation, Rspeed, with the observed month-to-month variability seen in Fig. 4 and found Rspeed = 0.81.

The third metric is a measure of how closely the modeled distribution of ice speeds reproduces the observed distribution of all 48 523 monthly mean grid cell speeds, χ2speed. Arranging all modeled and observed ice speeds separately into 0.5 cm s−1 bins from 0 to 20 cm s−1 (see Fig. 5), we first calculate the normalized observed distribution, Oi, where i = 1, . . . , 40 and Σ40i=1Oi = 1.χ2speed is then given by

 
formula

where Nm is the total number of modeled speeds that have been binned, and Mi is the number of modeled speeds in bin i, such that Σ40i=1Mi = Nm. Thus, lower χ2speed. values indicate better agreement, and we found that χ2speed=16 359 for the standard model.

Fig. 5.

Distribution of the observed and the standard model run's (dashed line) Arctic basin sea ice speeds from 1994 to 2001. As in Fig. 3, {α, Ca, P*} = {0.58, 0.0011, 10 kN m−2}. The distance between the distributions is χ2speed = 16 359.

Fig. 5.

Distribution of the observed and the standard model run's (dashed line) Arctic basin sea ice speeds from 1994 to 2001. As in Fig. 3, {α, Ca, P*} = {0.58, 0.0011, 10 kN m−2}. The distance between the distributions is χ2speed = 16 359.

c. Ice thickness comparisons

Before we could compare modeled and observed ice thickness, we first had to address the fact that the observations exclude areas of open water and ice of thickness h < 0.5 m (see section 2a). As already stated, we first minimized problems with thin ice by using only the ice thickness measurements made during the six winter months from November until April. However, we also addressed the problem explicitly in the model by treating ice thinner than 0.5 m as open water in each of the model's monthly grid cell diagnostics, specifically those of the model's first thickness category, which represents ice thinner than 0.64 m in our model. We then recalculated the monthly mean ice thickness for each grid cell before comparing the model predictions to the observed data.

In Fig. 6 we plot the time series of monthly Arctic basin (below 81.5°N) ice thickness from November 1993 until December 2001 to demonstrate the considerable interannual variability. These regional averages have a typical errors, Δhobst, as given by Eq. (1), of 14–15 cm.

Fig. 6.

The observed and standard model run's Arctic basin (below 81.5°N) sea ice thickness for all 50 winter months between 1994 and 2001 when {α, Ca, P*} = {0.58, 0.0011, 10 kN m−2}. For this model run we find hmod = 2.78 m and σmod = 0.91 m.

Fig. 6.

The observed and standard model run's Arctic basin (below 81.5°N) sea ice thickness for all 50 winter months between 1994 and 2001 when {α, Ca, P*} = {0.58, 0.0011, 10 kN m−2}. For this model run we find hmod = 2.78 m and σmod = 0.91 m.

To assess ice thickness predictions we compared hmod, the modeled, and hobs = 2.71 m, the observed mean of all 12 978 winter grid cell ice thicknesses measured between November 1993 and December 2001, as well as the corresponding standard deviations, σmod and σobs = 0.89 m. These metrics represent global measures of the accuracy of each model run's ice thickness and variability in the Arctic basin below 81.5°N.

The standard model run with ice thinner than 0.5 m excluded yielded values of hmod = 2.78 m and σmod = 0.91 m, whereas including the thin ice reduces the overall mean slightly to hmod = 2.66 m, but leaves the variability unchanged at σmod = 0.91 m.

Our third metric measures a model run's ability to accurately reproduce the observed monthly variability in area-averaged ice thickness seen in Figs. 6 and 7. We calculated

 
formula

where hmodt is the modeled area-averaged ice thickness in month t and Δhobs2t is its error variance given in Eq. (1), and for the standard model run found χ2thk = 175.

5. Model optimization

We now present the results of our search for a subset of parameter space in which the closest agreement with the observed extent, speed/velocity, and thickness data is found. For the reasons outlined in section 3, we have restricted our attention to three parameters in this study: Ca, the air–ice drag coefficient; P*, the ice strength parameter; and α, the broadband albedo of cold, bare ice. Furthermore, we restricted our search to the range of parameter values given in Table 2, chosen because they cover values used in previous studies. This resulted in a total of 168 model runs.

Table 2.

Parameter ranges explored in the optimization process.

Parameter ranges explored in the optimization process.
Parameter ranges explored in the optimization process.

We used the metrics introduced in section 4 to examine the predictions of ice extent, speed, and thickness for model runs using each of the 168 parameter combinations. Each dataset allowed us to progressively restrict the range of an optimal parameter set giving simultaneous agreement with all three datasets. We summarize the optimization process here, with full details provided in the appendix.

a. Summary of optimization procedure

We found (see appendix section a) that the value of α at which the lowest RMSDsep is reached depends upon the ice strength parameter. In particular, when P* ≤ 10 kN m−2, the lowest values of RMSDsep are found when α = 0.54–0.56 (Tables A1a and A1b); when 10 kN m−2 < P* < 100 kN m−2, the lowest RMSDsep values are found when α = 0.56–0.58 (Tables A1b and A1c); and when P* = 100 kN m−2, the lowest values are found when α = 0.60 (Table A1d). Figure 8 illustrates the variation in RMSDsep with P* and Ca, when α = 0.56. (Although the closest agreement is found where Ca = 0.0003, models run with this air drag coefficient predict ice speeds that are much lower than observed; see appendix section b.)

Table A1. Rms difference between observed and modeled September sea ice extent minima, RMSDsep (106 km2). The air drag parameter, Ca, varies vertically, and the ice strength parameter, P* (kN m−2), varies horizontally.

Table A1. Rms difference between observed and modeled September sea ice extent minima, RMSDsep (106 km2). The air drag parameter, Ca, varies vertically, and the ice strength parameter, P* (kN m−2), varies horizontally.
Table A1. Rms difference between observed and modeled September sea ice extent minima, RMSDsep (106 km2). The air drag parameter, Ca, varies vertically, and the ice strength parameter, P* (kN m−2), varies horizontally.
Fig. 8.

Sea ice extent optimization using values in Table A1b, where α = 0.56. The + symbol marks the optimized model run and the × symbol marks the standard model run. Solid lines: RMSDsep (106 km2); contour interval: 0.02 × 106 km2.

Fig. 8.

Sea ice extent optimization using values in Table A1b, where α = 0.56. The + symbol marks the optimized model run and the × symbol marks the standard model run. Solid lines: RMSDsep (106 km2); contour interval: 0.02 × 106 km2.

We also found (see appendix section b) that the ice speed metrics were largely insensitive to albedo (Tables A2a–d). When Ca = 0.0006, we required low P* values between 5 and 10 kN m−2 to obtain modeled speeds near the observed value, υobs. However, when we increased Ca to 0.000 85, we had to increase P* to approximately 27.5 kN m−2 to maintain the agreement. Increasing Ca still further, to 0.0011 and 0.0016, required us to raise P* to 55 and 100 kN m−2, respectively.

Table A2. (top entry) Modeled means of every monthly grid cell speed, υmod (cm s−1), in the Arctic basin from 1994 to 2001. (second entry) Correlation of observed and modeled monthly, basin-averaged ice speeds, Rspeed, from 1994 to 2001. (bottom entry) Measure of the difference between the modeled and observed ice speed distributions, χ2speed.

Table A2. (top entry) Modeled means of every monthly grid cell speed, υmod (cm s−1), in the Arctic basin from 1994 to 2001. (second entry) Correlation of observed and modeled monthly, basin-averaged ice speeds, Rspeed, from 1994 to 2001. (bottom entry) Measure of the difference between the modeled and observed ice speed distributions, χ2speed.
Table A2. (top entry) Modeled means of every monthly grid cell speed, υmod (cm s−1), in the Arctic basin from 1994 to 2001. (second entry) Correlation of observed and modeled monthly, basin-averaged ice speeds, Rspeed, from 1994 to 2001. (bottom entry) Measure of the difference between the modeled and observed ice speed distributions, χ2speed.

In addition, the monthly correlation, Rspeed, was found to be high (∼0.80) for all values of Ca when P* ≤ 10 kN m−2, but decreased with increasing P*. Furthermore, the lowest values of χ2speed were found where modeled speeds were near the observed υobs.

Figure 9 illustrates the variation in all three speed metrics with P* and Ca for α = 0.56, but the results of the sea ice speed optimization were found to be insensitive to α in the ranges we explored. Values of |υmod − υobs| ≤ 0.2 cm s−1 and χ2speed ≤ 5, 000 (near the + symbol) clearly coincide when both Ca and P* are raised simultaneously. However, Rspeed values (dotted lines) decrease as P* is raised and/or as Ca is lowered.

Fig. 9.

Sea ice speed optimization using values in Table A2b, where α = 0.56. The + symbol marks the optimized model run and the × symbol marks the standard model run. Solid lines: υmodυobs (cm s−1); contour interval: 0.2 cm s−1; dotted lines: Rspeed; dashed lines: χ2speed, variable contour interval.

Fig. 9.

Sea ice speed optimization using values in Table A2b, where α = 0.56. The + symbol marks the optimized model run and the × symbol marks the standard model run. Solid lines: υmodυobs (cm s−1); contour interval: 0.2 cm s−1; dotted lines: Rspeed; dashed lines: χ2speed, variable contour interval.

Combining the findings described in appendix sections a and b, we could conclude that simultaneous agreement with observed sea ice extent and speed was obtained in model runs with the following range of model parameter choices: {α, Ca, P*} = {0.54–0.56, 0.0006, 5–10 kN m−2}, {0.58, 0.000 85, 27.5 kN m−2}, {0.58, 0.0011, 55 kN m−2}, and {0.60, 0.0016, 100 kN m−2}. However, despite each of these model runs predicting ice extents and speeds in agreement with observations, their predictions of sea ice thickness and its variability differed widely (see appendix section c).

In particular, each such model run with P* ≥ 10 kN m−2 predicted a mean ice thickness that was too low and showed less variability than the observed data, but the agreement with observed ice thickness did improve when P* was reduced, at points {α, Ca, P*} = {0.54–0.56, 0.0006, 5–7.5 kN m−2}. As well as predicting mean ice thickness and its standard deviation close to the observed values of hmod = 2.78 m and σmod = 0.91 m, model runs with these parameter values also showed some of the lowest χ2thk values found anywhere in the parameter space we explored, that is, χ2thk ≈ 147–623, demonstrating that, not only was the absolute ice thickness being accurately modeled, but the variability in thickness was also.

Figure 10 illustrates the variation of ice thickness predictions with changing P* and Ca, for α = 0.56. Best agreement, that is, |hmodhobs| ≤ 0.2 m (solid lines), |σmodσobs| ≤ 0.2 m (dotted lines), and χ2thk ≤ 400 (dashed lines), is indicated in the figure by the + symbol and clearly found in an overlapping region of parameter space for all three measures when P* ≤ 10 kN m−2.

Fig. 10.

Sea ice thickness optimization using values in Table A3b, where α = 0.56. The + symbol marks the optimized model run and the × symbol marks the standard model run. Solid lines: hmodhobs (m); contour interval: 0.2 m; dotted lines: σmodσobs (m); contour interval: 0.2 m; dashed lines: χ2thk, variable contour interval.

Fig. 10.

Sea ice thickness optimization using values in Table A3b, where α = 0.56. The + symbol marks the optimized model run and the × symbol marks the standard model run. Solid lines: hmodhobs (m); contour interval: 0.2 m; dotted lines: σmodσobs (m); contour interval: 0.2 m; dashed lines: χ2thk, variable contour interval.

The optimum region in parameter space was restricted further by observing that ice thickness and its variability were still too low when P* = 7.5 kN m−2. Finally, using Tables A1a, A1b, A2a, A2b, A3a and A3b to compare the two remaining parameter choices, namely {α, Ca, P*} = {0.54 and 0.56, 0.0006, 5 kN m−2}, the model run with albedo α = 0.54 was found to agree slightly better with observed September sea ice extent minima (by 0.03 × 106 km2), but the model run with albedo α = 0.56 performed better on most measures of sea ice speed and thickness agreement. The point {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2} was therefore chosen as the final, optimal parameter choice, which in the following we describe as our “optimized model.”

Table A3. (top entry) hmod (m), (second entry) σmod (m), and (third entry) χ2thk, for various values of the air drag parameter, Ca, and ice strength, P* (kN m−2).

Table A3. (top entry) hmod (m), (second entry) σmod (m), and (third entry) χ2thk, for various values of the air drag parameter, Ca, and ice strength, P* (kN m−2).
Table A3. (top entry) hmod (m), (second entry) σmod (m), and (third entry) χ2thk, for various values of the air drag parameter, Ca, and ice strength, P* (kN m−2).

In Fig. 11 we have used one measure each of the model run predictions of thickness, speed, and extent to illustrate the simultaneity of agreement for the optimized model, indicated in the figure by the + symbol. The standard model run, on the other hand, indicated in the figure by the × symbol, gives satisfactory predictions of sea ice thickness and extent but overestimates mean sea ice speed by 1 cm s−1.

Fig. 11.

Model optimization, illustrated using one metric to gauge thickness, speed, and extent predictions, where α = 0.56. Red lines: hmodhobs (m); contour interval: 0.2 m; green lines: υmodυobs (cm s−1); contour interval: 0.2 cm s−1; blue lines: RMSDsep (106 km2); contour interval: 0.02 × 106 km2. Optimal values are indicated with bold lines. The + symbol marks the optimized model run and the × symbol marks the standard model run. We also indicate parameter values used in previous studies that use the same air drag and ice strength parameterizations: □, Hibler (1979); ▵, Holland et al. (1993); ⋄, Hibler and Walsh (1982); *, Harder and Fischer (1999). Note, however, that the albedo parameterizations and values used in these studies differ.

Fig. 11.

Model optimization, illustrated using one metric to gauge thickness, speed, and extent predictions, where α = 0.56. Red lines: hmodhobs (m); contour interval: 0.2 m; green lines: υmodυobs (cm s−1); contour interval: 0.2 cm s−1; blue lines: RMSDsep (106 km2); contour interval: 0.02 × 106 km2. Optimal values are indicated with bold lines. The + symbol marks the optimized model run and the × symbol marks the standard model run. We also indicate parameter values used in previous studies that use the same air drag and ice strength parameterizations: □, Hibler (1979); ▵, Holland et al. (1993); ⋄, Hibler and Walsh (1982); *, Harder and Fischer (1999). Note, however, that the albedo parameterizations and values used in these studies differ.

b. Output from the optimized model

In Figs. 12 to 17 we have plotted the optimized model's sea ice thickness, extent, and speed predictions against the observed data.

Fig. 12.

Observed Arctic basin (below 81.5°N) sea ice thickness for all 50 winter months between 1994 and 2001 and the predictions of the optimized model run, where {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2}. The modeled thickness hmod = 2.69 m and its variability σmod = 0.96 m.

Fig. 12.

Observed Arctic basin (below 81.5°N) sea ice thickness for all 50 winter months between 1994 and 2001 and the predictions of the optimized model run, where {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2}. The modeled thickness hmod = 2.69 m and its variability σmod = 0.96 m.

Agreement with observed thickness in Fig. 12 is close before 1999 after which the monthly modeled ice is thicker than observed by up to 35 cm but, in general, the observed inter- and intra-annual variability is captured. The quality of the spatial agreement can be gauged by inspecting Fig. 13, where we have plotted the observed mean Arctic basin (below 81.5°N) ice thickness field for winters 1993–2001 and the output of both the optimized model and the standard model run. Both model runs have similar mean ice thickness fields (though, as stated in section 5a, the standard model run's ice speeds are too high) and both overestimate observed ice thickness in the western Arctic and underestimate it in the eastern Arctic, particularly in the Laptev Sea. This may be due to the simple ice strength and air drag parameterizations used and also the absence of a coupled ocean model in these shallow waters where tides and river runoff would be expected to play an important role.

Fig. 13.

Mean Arctic basin (below 81.5°N) sea ice thickness for winters 1993–2001 (top) as observed, (middle) as predicted by the optimized model, and (bottom) as predicted by the standard model.

Fig. 13.

Mean Arctic basin (below 81.5°N) sea ice thickness for winters 1993–2001 (top) as observed, (middle) as predicted by the optimized model, and (bottom) as predicted by the standard model.

Figure 14 demonstrates that the optimized model is largely able to capture the considerable interannual variability in summer ice extent, and in Fig. 15 we have plotted the observed and modeled September ice extent anomalies to illustrate this in more detail. We find that the time series are highly correlated (R = 0.98), with the largest discrepancies occurring in the summer of 1994 when the model underestimates the observed decrease in summer ice extent.

Fig. 14.

Observed Arctic basin sea ice extent from 1994 to 2001, and the predictions of the optimized model, where {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2}. We find that RMSDsep = 0.15 × 106 km2.

Fig. 14.

Observed Arctic basin sea ice extent from 1994 to 2001, and the predictions of the optimized model, where {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2}. We find that RMSDsep = 0.15 × 106 km2.

Fig. 15.

Observed Arctic basin 1994–2001 September sea ice extent anomalies and the predictions of the optimized model, where {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2}. The correlation R = 0.98.

Fig. 15.

Observed Arctic basin 1994–2001 September sea ice extent anomalies and the predictions of the optimized model, where {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2}. The correlation R = 0.98.

Figure 16 shows the high correlation (Rspeed = 0.81) between the optimized model's predictions of mean Arctic Ocean ice speed and the observed Polar Pathfinder data (Fowler 2003). The greatest discrepancies generally arise during the summer months. For example, the model overestimates summer ice speeds in 1997 and 1999 but underestimates them in 1994.

Fig. 16.

Observed Arctic basin sea ice speeds from 1994 to 2001 and the predictions of the optimized model, where {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2}. The correlation Rspeed = 0.81 and modeled mean speed υmod = 3.01 cm s−1.

Fig. 16.

Observed Arctic basin sea ice speeds from 1994 to 2001 and the predictions of the optimized model, where {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2}. The correlation Rspeed = 0.81 and modeled mean speed υmod = 3.01 cm s−1.

Finally, in Fig. 17 we use the modeled and observed ice speed distributions to show that the model predicts too little ice with speeds between 0 and 2 cm s−1 and too much with speeds between 2 and 5 cm s−1. However, the agreement for speeds above 5 cm s−1 is close.

Fig. 17.

Distribution of all observed Arctic basin sea ice speeds from 1994 to 2001, and the predictions of the optimized model, where {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2}. The difference between the distributions χ2speed = 4969.

Fig. 17.

Distribution of all observed Arctic basin sea ice speeds from 1994 to 2001, and the predictions of the optimized model, where {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2}. The difference between the distributions χ2speed = 4969.

c. Optimal parameter sensitivity

The model tuning process just described is valid for this particular model forced with the data described in section 3c. It has highlighted the importance of tuning and validating models, however complex, with multiple datasets. It is of course possible that an identical model forced with alternative datasets, such as the National Centers for Environmental Prediction (NCEP) reanalyses, would require different parameter values to achieve the same level of agreement with observations.

To investigate the sensitivity of the optimal parameters to alternative forcing data, we have rerun the optimized model with its ocean current field set to zero, that is, with the ocean acting as a motionless drag on the ice. We found a slight increase in the model's thickness predictions (hmod = 2.80 m, σmod = 0.88 m, χ2thk = 178), generally lower speed predictions (υmod = 2.33 cm s−1, Rspeed = 0.79, χ2speed = 5140), and a significant increase in the September sea ice extent, with RMSDsep = 0.35 × 106 km2. However, changing the model parameters to {α, Ca, P*} = {0.52, 0.007 25, 5 kN m−2} improved the model's thickness predictions (hmod = 2.64 m, σmod = 0.90 m, χ2thk = 201), its speed predictions (υmod = 2.66 cm s−1, Rspeed = 0.79, χ2speed = 4766), and its predictions of September sea ice extent (RMSDsep = 0.20 × 106 km2).

Thus, although the POLES ocean current field used (see section 3) is clearly idealistic and lacks temporal variability, treating the ocean as a source of passive drag instead does not require a large change in the model's optimal parameters.

6. Summary and conclusions

We have used basin-scale observations of Arctic sea ice thickness, extent, and velocity to tune and validate a modified, stand-alone, version of the CICE sea ice model (Hunke and Lipscomb 2004) forced with POLES and ERA-40 reanalysis data. This is the first time that basin-scale observations of sea ice thickness (Laxon et al. 2003) have been used for this purpose.

Three parameters were varied in the tuning process: Ca, the air–ice drag coefficient; P*, the ice strength parameter; and α, the broadband albedo of cold bare ice—parameters chosen because of their large influence on sea ice thickness, extent, and velocity, and because there remains considerable uncertainty regarding the correct values to use in basin-scale sea ice models (Hibler 1979; Chapman et al. 1994). After running the model with 168 parameter choices covering a subset of the full three-dimensional parameter space and comparing the model run predictions of ice extent, speeds, and sea ice thickness to observations, we were able to specify the point {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2} as the optimal parameter choice for this particular forcing set (see Figs. 12 –17). However, as described in section 5b, there remained differences between the optimized model run and all three sets of observations.

One possible reason for these discrepancies is the existence of errors in the forcing data. Furthermore, some sea ice processes are either inadequately parameterized or omitted entirely from the model. CICE has an ice thickness distribution (Thorndike et al. 1975), an EVP rheology (Hibler 1979; Hunke and Dukowicz 1997, 2002), and multilayer, energy-conserving thermodynamics (Bitz and Lipscomb 1999), but we changed its original, sophisticated parameterizations of air–ice drag and ice strength (Flato and Hibler 1995) to those of Eqs. (4) and (7), respectively. Future work will involve restoring these complex parameterizations to examine whether it is possible to improve upon the predictions of the optimized model.

Still greater improvements should be possible by coupling the sea ice model to an ocean model and by increasing its resolution. Future enhancements of the model should also involve explicit parameterizations of melt ponds (Eicken et al. 2004; Taylor and Feltham 2004) and removal of the implicit assumption of lead isotropy (Wilchinsky and Feltham 2004).

Notwithstanding these probable sources of error, the availability of basinwide observational datasets of sea ice concentration, velocity (NSIDC), and thickness (Laxon et al. 2003) has enabled us to tune and validate a sea ice model and to reduce ambiguity in commonly used parameterizations of air–ice drag, ice strength, and albedo in a model forced with a combination of ERA-40 and POLES data. We located an optimal region in parameter space where a model run's predictions agreed satisfactorily and simultaneously with the observed data. Different forcing can lead to different optimal parameters (see section 5c), but we have demonstrated the importance of tuning and validating models, however complex, with multiple datasets.

Acknowledgments

We are grateful to the model developers and to the Los Alamos National Laboratory for making CICE freely available to the scientific community. We thank the European Space Agency for the provision of the ERS data, and Neil Peacock and Rob Potter for processing it. We would also like to thank the National Snow and Ice Data Center for providing sea ice concentration and velocity data, and the IABP/POLES teams for providing additional forcing data. The authors are grateful to two anonymous referees for their helpful comments. P. M. was funded by the NERC COAPEC thematic program.

REFERENCES

REFERENCES
Bitz
,
C. M.
, and
W. H.
Lipscomb
,
1999
:
An energy-conserving thermodynamic sea ice model for climate study.
J. Geophys. Res
,
104
,
15669
15677
.
Bourke
,
R. H.
, and
A. S.
McLaren
,
1992
:
Contour mapping of Arctic basin ice draft and roughness parameters.
J. Geophys. Res
,
97
,
4533
4544
.
Carson
,
D. J.
,
1999
:
Climate modelling: Achievements and prospects.
Quart. J. Roy. Meteor. Soc
,
125
,
1
27
.
Cavalieri
,
D. J.
,
C. L.
Parkinson
,
P.
Gloersen
, and
H. J.
Zwally
,
2002
:
Sea ice concentrations from Nimbus-7 SMMR and DMSP SSM/I passive microwave data. National Snow and Ice Data Center, Boulder, CO, CD-ROM
.
Chapman
,
W. L.
,
W. J.
Welch
,
K. P.
Bowman
,
J.
Sacks
, and
J. E.
Walsh
,
1994
:
Arctic sea ice variability: Model sensitivities and a multidecadal simulation.
J. Geophys. Res
,
99
,
919
935
.
Comiso
,
J. C.
, and
R.
Kwok
,
1996
:
Surface and radiative characteristics of the summer Arctic sea ice cover from multi-sensor satellite observations.
J. Geophys. Res
,
101
,
28397
28416
.
Curry
,
J. C.
,
J. L.
Schramm
, and
D. K.
Perovich
,
2001
:
Applications of SHEBA/FIRE data to evaluation of snow/ice albedo parameterizations.
J. Geophys. Res
,
106
,
15345
15355
.
Curry
,
J. C.
,
J. L.
Schramm
,
A.
Alam
,
R.
Reeder
,
T. E.
Arbetter
, and
P.
Guest
,
2002
:
Evaluation of data sets used to force sea ice models in the Arctic Ocean.
J. Geophys. Res
,
107
.
8027, doi:10.1029/2000JC000466
.
Ebert
,
E. E.
, and
J. A.
Curry
,
1993
:
An intermediate one-dimensional thermodynamic model for investigating ice-atmosphere interactions.
J. Geophys. Res
,
98
,
10085
10109
.
Eicken
,
H.
,
T. C.
Grenfell
,
D. K.
Perovich
,
J. A.
Richter-Menge
, and
K.
Frey
,
2004
:
Hydraulic controls of summer Arctic pack ice albedo.
J. Geophys. Res
,
109
.
C08007, doi:10.1029/2003JC001989
.
Fischer
,
H.
, and
P.
Lemke
,
1994
:
On the required accuracy of atmospheric forcing fields for driving dynamic-thermodynamic sea ice models. The Polar Oceans and Their Role in Shaping the Global Environment, Geophys. Monogr., No. 85, Amer. Geophys. Union, 373–381
.
Flato
,
G. M.
,
1995
:
Spatial and temporal variability of Arctic ice thickness.
Ann. Glaciol
,
21
,
323
329
.
Flato
,
G. M.
, and
W. D.
Hibler
III
,
1995
:
Ridging and strength in modeling the thickness distribution of Arctic sea ice.
J. Geophys. Res
,
100
,
18611
18626
.
Fowler
,
C.
,
2003
:
Polar pathfinder daily 25km EASE-grid sea ice motion vectors. National Snow and Ice Data Center, Boulder, CO, digital media. [Available online at http://nside.org/data/nside-0166.html.]
.
Harder
,
M.
, and
H.
Fischer
,
1999
:
Sea ice dynamics in the Weddell Sea simulated with an optimized model.
J. Geophys. Res
,
104
,
11151
11162
.
Hibler
III,
W. D.
,
1979
:
A dynamic–thermodynamic sea ice model.
J. Phys. Oceanogr
,
9
,
815
846
.
Hibler
III,
W. D.
,
2004
:
Modelling the dynamic response of sea ice.
Mass Balance of the Cryosphere, J. L. Bamber and A. J. Payne, Eds., Cambridge University Press, 227–334
.
Hibler
III,
W. D.
, and
J. E.
Walsh
,
1982
:
On modeling the seasonal and interannual fluctuations of Arctic sea ice.
J. Phys. Oceanogr
,
12
,
1514
1523
.
Hibler
III,
W. D.
, and
K.
Bryan
,
1984
:
Ocean circulation: Its effects on seasonal sea-ice simulations.
Science
,
224
,
489
492
.
Hibler
III,
W. D.
, and
K.
Bryan
,
1987
:
A diagnostic ice–ocean model.
J. Phys. Oceanogr
,
7
,
987
1015
.
Holland
,
D. M.
,
L. A.
Mysak
, and
D. K.
Manak
,
1993
:
Sensitivity study of a dynamic thermodynamic sea ice model.
J. Geophys. Res
,
98
,
2561
2583
.
Houghton
,
J. T.
,
Y.
Ding
,
D. J.
Griggs
,
M.
Noguer
,
P. J.
van der Linden
,
X.
Dai
,
K.
Maskell
, and
C. A.
Johsnon
,
2001
:
Climate Change 2001: The Scientific Basis.
Cambridge University Press, 881 pp
.
Hunke
,
E. C.
, and
J. K.
Dukowicz
,
1997
:
An elastic–viscous–plastic model for sea ice dynamics.
J. Phys. Oceanogr
,
27
,
1849
1867
.
Hunke
,
E. C.
, and
Y.
Zhang
,
1999
:
A comparison of sea ice dynamics models at high resolution.
Mon. Wea. Rev
,
127
,
396
408
.
Hunke
,
E. C.
, and
J. K.
Dukowicz
,
2002
:
The elastic–viscous–plastic sea ice dynamics model in general orthogonal curvilinear coordinates on a sphere—Incorporation of metric terms.
Mon. Wea. Rev
,
130
,
1848
1865
.
Hunke
,
E. C.
, and
W. H.
Lipscomb
,
2004
:
CICE: The Los Alamos Sea Ice Model, documentation and software user's manual. Los Alamos National Laboratory. [Available online at http://climate.lanl.gov/Models/CICE/index.htm.]
.
Johns
,
T.
, and
Coauthors
,
2004
:
HadGEM1—Model description and analysis of preliminary experiments for the IPCC Fourth Assessment Report. Hadley Centre Tech. Note 55, Met Office, Exeter, United Kingdom, 75 pp. [Available online at http://www.metoffice.gov.uk/research/hadleycentre/pubs/HCTN/index.html.]
.
Jordan
,
R. E.
,
E. L.
Andreas
, and
A. P.
Makshtas
,
1999
:
Heat budget of snow-covered sea ice at North Pole 4.
J. Geophys. Res
,
104
,
7785
7806
.
Kreyscher
,
M.
,
M.
Harder
,
P.
Lemke
, and
G. M.
Flato
,
2000
:
Results of the Sea Ice Model Intercomparison Project: Evaluation of sea ice rheology schemes for use in climate simulations.
J. Geophys. Res
,
105
,
11299
11320
.
Laxon
,
S.
,
N.
Peacock
, and
D.
Smith
,
2003
:
High interannual variability of sea ice thickness in the Arctic region.
Nature
,
425
,
947
950
.
Lindsay
,
R. W.
, and
H. L.
Stern
,
2004
:
A new Lagrangian model of Arctic sea ice.
J. Phys. Oceanogr
,
34
,
272
283
.
Maykut
,
G. A.
, and
M. G.
McPhee
,
1995
:
Solar heating of the Arctic mixed layer.
J. Geophys. Res
,
100
,
24691
24703
.
McPhee
,
M. G.
,
1975
:
Ice–ocean momentum transfer for the AIDJEX ice model.
AIDJEX Bull
,
29
,
93
111
.
Parkinson
,
C. L.
,
D. J.
Cavalieri
,
P.
Gloersen
,
H. J.
Zwally
, and
J. C.
Comiso
,
1999
:
Arctic sea ice extents, areas, and trends, 1978–1996.
J. Geophys. Res
,
104
,
20837
20856
.
Partington
,
K.
,
T.
Flynn
,
D.
Lamb
,
C.
Bertoia
, and
K.
Dedrick
,
2003
:
Late twentieth century Northern Hemisphere sea-ice record from U.S. National Ice Center ice charts.
J. Geophys. Res
,
108
.
3343, doi:10.1029/2002JC001623
.
Perovich
,
D. K.
, and
Coauthors
,
1999
:
Year on ice gives climate insights.
Eos, Trans. Amer. Geophys. Union
,
80
,
481
485–486
.
Perovich
,
D. K.
,
T. C.
Grenfell
,
B.
Light
, and
P. V.
Hobbs
,
2002
:
Seasonal evolution of the albedo of multiyear Arctic sea ice.
J. Geophys. Res
,
107
.
8044, doi:10.1029/2000JC000438
.
Perovich
,
D. K.
,
T. C.
Grenfell
,
J. A.
Richter-Menge
,
B.
Light
,
W. B.
Tucker
III
, and
H.
Eicken
,
2003
:
Thin and thinner: Sea ice mass balance measurements during SHEBA.
J. Geophys. Res
,
108
.
8050, doi:10.1029/2001JC001079
.
Rigor
,
I. G.
,
R. L.
Colony
, and
S.
Martin
,
2000
:
Variations in surface air temperature observations in the Arctic 1979–97.
J. Climate
,
13
,
896
914
.
Rind
,
D.
,
R.
Healy
,
C.
Parkinson
, and
D.
Martinson
,
1997
:
The role of sea ice in 2 × CO2 climate model sensitivity: Part II: Hemispheric dependencies.
Geophys. Res. Lett
,
24
,
1491
1494
.
Rothrock
,
D. A.
,
1975
:
The energetics of the plastic deformation of pack ice by ridging.
J. Geophys. Res
,
80
,
4514
4519
.
Rothrock
,
D. A.
,
J.
Zhang
, and
Y.
Yu
,
2003
:
The Arctic ice thickness of the 1990s: A consistent view from observations and models.
J. Geophys. Res
,
108
.
3083, doi:10.1029/2001JC001208
.
Serreze
,
M. C.
, and
Coauthors
,
2003
:
A record minimum Arctic sea ice extent and area in 2002.
Geophys. Res. Lett
,
30
.
1110, doi:10.1029/2002GL016406
.
Steele
,
M.
,
J.
Zhang
,
D.
Rothrock
, and
H.
Stern
,
1997
:
The force balance of sea ice in a numerical model of the Arctic Ocean.
J. Geophys. Res
,
102C
,
21061
21079
.
Steele
,
M.
,
R.
Morley
, and
W.
Ermold
,
2001
:
PHC: Global ocean hydrography with a high-quality Arctic Ocean.
J. Climate
,
14
,
2079
2087
.
Taylor
,
P. D.
, and
D. L.
Feltham
,
2004
:
A model of melt pond evolution on sea ice.
J. Geophys. Res
,
109
.
C12007, doi:10.1029/2004JC002361
.
Thorndike
,
A. S.
,
D. A.
Rothrock
,
G. A.
Maykut
, and
R.
Colony
,
1975
:
The thickness distribution of sea ice.
J. Geophys. Res
,
80
,
4501
4513
.
Vowinckel
,
E.
, and
S.
Orvig
,
1970
:
The climate of the North Polar Basin.
Climates of the Polar Regions, S. Orvig, Ed., World Survey of Climatology, Vol. 14, Elsevier, 129–252
.
Walter
Jr.,
B. A.
,
J. E.
Overland
, and
R.
Gilmer
,
1984
:
Air–ice drag coefficients for first year sea ice derived from aircraft measurements.
J. Geophys. Res
,
89
,
3550
3560
.
Warren
,
S. G.
,
I. G.
Rigor
,
N.
Untersteiner
,
V. F.
Radionov
,
N. N.
Bryazgin
,
Y. I.
Aleksandrov
, and
R.
Colony
,
1999
:
Snow depth on Arctic sea ice.
J. Climate
,
12
,
1814
1829
.
Wilchinsky
,
A.
, and
D.
Feltham
,
2004
:
A continuum anisotropic model of sea ice dynamics.
Proc. Roy. Soc. London
,
460A
,
2105
2140
.
Zhang
,
J.
, and
D. A.
Rothrock
,
2003
:
Modeling global sea ice with a thickness and enthalpy distribution model in generalized curvilinear coordinates.
Mon. Wea. Rev
,
131
,
845
861
.
Zhang
,
J.
,
D. A.
Rothrock
, and
M.
Steele
,
1998
:
Warming of the Arctic Ocean by a strengthened Atlantic inflow: Model results.
Geophys. Res. Lett
,
25
,
1745
1748
.
Zhang
,
Y.
, and
E. C.
Hunke
,
2001
:
Recent Arctic change simulated with a coupled ice-ocean model.
J. Geophys. Res
,
106
,
4369
4390
.

APPENDIX

An Optimization of the Sea Ice Model Using Observations of Extent, Velocity/Speed, and Thickness

Sea ice extent optimization

Tables A1a–e show the values of the ice extent metric RMSDsep for α = 0.54, 0.56, 0.58, 0.60, and 0.62, respectively. Each subtable shows a “slice” through parameter space at a fixed α.

The eight blank entries in these tables indicate parameter values with which the model did not run to completion. In every case, this happened when using the highest air drag value in our parameter space, Ca = 0.0016, in combination with some of the lowest ice strength values, P* ≤ 10 kN m−2, and was a result of excessive unphysical ridging.

We found the highest values of RMSDsep when α = 0.62, irrespective of Ca and P*, so we will not consider this value of albedo in what follows. Similarly, we found the lowest values of RMSDsep where Ca = 0.0003 (see Table A1b). However, for these runs the ice speeds were too low (see appendix section b), so we did not run models with Ca = 0.0003 for albedo values other than α = 0.56.

When P* ≤ 10 kN m−2, the lowest values of RMSDsep were found when α = 0.54–0.56 (Tables A1a and A1b); when 10 kN m−2 < P* < 100 kN m−2, the lowest RMSDsep values were found when α = 0.56–0.58 (Tables A1b and A1c); and, when P* = 100 kN m−2, the lowest values were found when α = 0.60 (Table A1d).

Sea ice velocity/speed optimization

Tables A2a–d show the results of the ice speed assessment for α = 0.54, 0.56, 0.58, and 0.60, respectively. The first (top) entry is υmod, the second entry is Rspeed, and the third (bottom) entry is χ2speed.

We note that for each fixed value of α and P*, raising Cair increases υmod, and when α and Cair are fixed we find that υmod decreases as P* increases. In addition, albedo was found to have only a small and indirect effect on the ice speeds.

When Ca = 0.0006 we required low P* values of between 5 and 10 kN m−2 to obtain modeled speeds near the observed υobs. However, increasing Ca to 0.00085 requires a P* of approximately 27.5 kN m−2 to maintain the agreement. Increasing Ca still further, to 0.0011 and 0.0016, required us to raise P* to 55 and 100 kN m−2, respectively.

We found high correlations (Rspeed ≈ 0.80) for all values of Ca when P* ≤ 10 kN m−2, which decreased as P* increased. Also, the lowest values of χ2speed were found when the ice strength was high, with little variation with α. However, low values of χ2speed that is, χ2speed ≤ 5300, were found whenever modeled speeds, υmod, were close to the observed value υobs.

Taking all three measures into consideration, closest agreement was found at points {Ca, P*} = {0.0006, 5–10 kN m−2}, {0.000 85, 27.5 kN m−2}, {0.0011, 55 kN m−2}, and {0.0016, 100 kN m−2}. Furthermore, since the results of the sea ice speed optimization given above are insensitive to α in the ranges we have explored, we can we now combine these findings with those of appendix section a to conclude that simultaneous agreement with observed sea ice extent and speed is found at the following points: {α, Ca, P*} = {0.54–0.56, 0.0006, 5–10 kN m−2}, {0.58, 0.000 85, 27.5 kN m−2}, {0.58, 0.0011, 55 kN m−2}, and {0.60, 0.0016, 100 kN m−2}, which still leaves considerable ambiguity in the choice of an optimum parameter set. To reduce this uncertainty we need a third, independent dataset.

Sea ice thickness optimization

In Tables A3a–d we give the results of our ice thickness metrics for α = 0.54, 0.56, 0.58, and 0.60, respectively. The first (top) entry is hmod, the second entry is σmod, and the third (bottom) entry is χ2mod.

As we found in appendix section b, we were able to achieve satisfactory agreement with sea ice extent and velocity near the points {α, Ca, P*} = {0.54–0.56, 0.0006, 5–10 kN m−2}, {0.58, 0.000 85, 27.5 kN m−2}, {0.58, 0.0011, 55 kN m−2}, and {0.60, 0.0016, 100 kN m−2}. We have highlighted these points in boldface in Tables A3a–d.

Inspection of those points in Table A3c and Table A3d quickly reveals that despite the close agreement with observed extent and velocity, the agreement with observed sea ice thickness at higher values of P* is poor. Modeled ice at these points is too thin, and it shows less variability than the observed data. Similarly, inspection of highlighted points when P* = 10 kN m−2 in Table A3a and Table A3b shows that these runs also predict ice that is both too thin by more than 0.5 m and insufficiently variable.

However, the agreement with observed ice thickness improves near {α, Ca, P*} = {0.54–0.56, 0.0006, 5–7.5 kN m−2} (see Tables A3a and A3b). As well as predicting ice thickness and its standard deviation close to the observed values of hmod = 2.78 m, and σmod = 0.91 m, these regions also show some of the lowest χ2thk values found anywhere in the parameter space we have explored, χ2thk ≈ 147–623, demonstrating that not only is the absolute ice thickness being accurately modeled, but the intra-annual variability in thickness is also.

The optimum region in parameter space has therefore been reduced to the region near {α, Ca, P*} = {0.54–0.56, 0.0006, 5–7.5 kN m−2}. However, we can restrict this set further by first observing that ice thickness and its variability are still too low when P* = 7.5 kN m−2. Finally, using Tables A1a, A1b, A2a, A2b, A3a and A3b to compare the two remaining points, namely {α, Ca, P*} = {0.54–0.56, 0.0006, 5 kN m−2}, we see that the model run with albedo α = 0.54 agrees slightly better with observed September sea ice extent minima (by 0.03 × 106 km2), but the model run with albedo α = 0.56 performs better on most measures of sea ice speed and thickness agreement. We therefore choose the point {α, Ca, P*} = {0.56, 0.0006, 5 kN m−2} to be our optimal parameter set.

Footnotes

Corresponding author address: Dr. Seymour W. Laxon, Centre for Polar Observation and Modelling, Dept. of Earth Sciences, University College London, Gower Street, London WC1E 6BT, United Kingdom. Email: swl@cpom.ucl.ac.uk