## Abstract

Model simulations of Arctic sea ice and ocean systems are used to determine the major spatial and temporal modes of variability in the ice thickness. A coupled ice–ocean model is forced with daily NCEP–NCAR reanalysis surface air pressure and surface air temperature fields for the period 1951–2003 with the analysis of the results performed for the 51-yr period 1953–2003. Ice concentration data and ice velocity data (beginning in 1979) are assimilated to further constrain the simulations to match the observed conditions. The simulated ice thins over the study period with the area of greatest thinning in a band from the Laptev Sea across the Pole to Fram Strait. The thinning rate is greatest since 1988. The major spatial modes of variability were determined with empirical orthogonal functions (EOFs) for the ice thickness within the Arctic Ocean. The first three EOFs account for 30%, 18%, and 15%, respectively, of the annual mean ice thickness variance. The first EOF is a nearly basinwide pattern, and the next two are orthogonal lateral modes. Because of the nonstationary nature of the ice thickness time series, significant changes in the modes are found if a shorter period is analyzed. The second and third principal components are well correlated with the Arctic Oscillation. The model results are also used to simulate an observation system and to then determine optimal mooring locations to monitor the basinwide mean ice thickness as well as the spatial and temporal patterns represented in the EOF analysis. The nonstationary aspect of the ice thickness limits the strength of the conclusions that can be drawn.

## 1. Introduction

Arctic sea ice is in a period of rapid change, both in ice extent (e.g., Serreze et al. 2003) and in ice thickness (e.g., Rothrock et al. 1999). The ice extent and thickness are both important climate-state variables in the global climate system because of the strong dependence of the surface sensible, latent, and radiative heat fluxes on the presence and thickness of the ice. We wish to better understand where and by how much these parameters are changing in order to better determine causes of the changes and in order to design efficient monitoring systems to track them. One technique to determine the major modes of variability of the ice extent or thickness is through the use of empirical orthogonal functions (EOFs).

An EOF analysis of the weekly ice concentration as represented in the National Ice Center (NIC) ice charts from the period 1972–94 was made by Partington et al. (2003). They report that the major winter mode is a dipole of ice concentration between the Barents and Greenland Seas and the Labrador Sea. It is well correlated with the North Atlantic Oscillation from the previous winter. Another prominent mode is for the summer ice. It represents a dipole between the Kara Sea and the East Siberian Sea. It is weakly correlated with the North Atlantic Oscillation (NAO) from the previous winter. Singarayer and Bamber (2003) used EOF analysis of the ice concentration to compare the modes of variability from three ice concentration datasets: the NIC ice charts and two passive microwave algorithms. They found modes similar to those reported by Partington et al. (2003) and the major modes of variability in the three datasets were essentially the same.

Ice concentration variability is primarily located along the edges of the ice pack. As a consequence the modes of variability are concentrated in the edge regions and the marginal ice zones. The ice thickness, however, varies continuously throughout the pack. The modes of variability are not constrained to be near the ice edges and are very distinct from the ice concentration modes.

Direct observations of the ice thickness are few and far between and are far too sparse to permit an analysis of the modes of variability. Estimates of the ice thickness, however, can be made with a model. We use a coupled ice–ocean model that is forced with observed surface pressure and temperature fields and, in addition, uses data assimilation methods for ice concentration and ice velocity. Consequently the simulated ice thickness and concentration are thought to be reasonably faithful recreations of the actual parameters. The simulated ice thickness can then be used to determine the major modes of variability of the ice thickness and to determine how the weighting of the modes change over time.

Models alone are not sufficient because of uncertainties in their representations of the ice thickness. Direct observations of the ice thickness distribution are essential to document the true evolution of the air–ice–ocean system in a changing global climate. A long-term observing system could help to verify what the interannual variability of the total ice volume is and how the regional distribution of the ice mass changes from year to year. Answering these questions will both document the changing state of the Arctic system and provide information to validate and improve coupled ice–ocean models.

Current coupled ice–ocean models can aid in the design of an observing system. Model simulations of the ice thickness can help to estimate where to place moorings with upward-looking sonars to most efficiently monitor the thickness of the ice pack. In determining the best locations to place new moorings, it is prudent to consider ongoing measurement programs. Two considered here are the North Pole Environmental Observatory (Morison et al. 2002), which is located very close to the Pole and was first established in the spring of 2001, and a second one established in the Chukchi Sea in 2003. Other mooring programs might also be considered, such as those in the Fram Strait or in the eastern Beaufort Sea.

Here we study the variability of the ice thickness as represented in the model, focusing on the interannual variability of the annual mean thickness for the 51-yr period, 1953–2003. We determine the local trends in the ice thickness and the principal modes of variability. We use the model results as a simulated observation system to find how the local mean thickness is related to the basinwide mean thickness and hence determine the locations where measurements of the ice draft are most representative of the basinwide thickness changes. A different set of points are optimal for determining the temporal variation in the major spatial patterns of ice thickness. Here we focus on determining the modes of variability of the ice thickness; a companion paper (Lindsay and Zhang 2005) considers the dominant physical processes that caused the dramatic thinning of the ice from 1988 to 2003.

The model and data assimilation methods are briefly described in section 2. In section 3 the variability of the ice thickness over the 51-yr period is analyzed, including temporal trends and EOF analysis. The model results are also used in section 4 to determine the locations where the basinwide mean ice thickness is best monitored or, alternatively, where to best monitor the spatial and temporal variability of the annual mean thickness. Comments and conclusions are given in section 5.

## 2. Model description and data assimilation methods

We use a coupled ice–ocean model that has been used in a wide range of studies. The ice model is a multicategory ice thickness and enthalpy distribution model that consists of five main components: 1) a momentum equation that determines ice motion, 2) a viscous–plastic ice rheology with an elliptical yield curve that determines the relationship between ice deformation and internal stress, 3) a heat equation that determines ice temperature and ice growth or decay, 4) two ice thickness distribution equations for deformed and undeformed ice that conserve ice mass, and 5) an enthalpy distribution equation that conserves ice thermal energy. The first two components are described in detail by Hibler (1979). The ice momentum equation is solved using the Zhang and Hibler (1997) numerical model for ice dynamics. The heat equation is solved, over each category, using the Winton (2000) three-layer thermodynamic model, which divides the ice in each category into two layers of equal thickness beneath a layer of snow. The ice thickness distribution equations are described in detail by Flato and Hibler (1995). The ocean model is based on the Bryan–Cox model (Bryan 1969; Cox 1984) with an embedded mixed layer from Kraus and Turner (1967). Detailed information about the ocean model is in Zhang et al. (1998).

The model domain covers the Arctic Ocean, and the Barents, and Greenland–Iceland–Norwegian Seas. It has a horizontal resolution of 40 km × 40 km, 21 vertical ocean levels, and 12 thickness categories each for undeformed ice, ridged ice, ice enthalpy, and snow. The ice thickness categories and bottom topography are in Zhang et al. (2000). The model domain is illustrated in Fig. 1 and the region of primary interest, the Arctic Ocean, is marked.

The model is forced with surface wind and temperature records from the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis (Kalnay et al. 1996). The NCEP–NCAR reanalysis is the result of a global atmospheric weather prediction model that assimilates all available weather data to estimate the state of the global atmosphere. Here we use only the daily averaged sea level pressure (SLP) and the 2-m air temperature fields (T2m) for the 56-yr period, 1948–2003. The specific humidity and longwave and shortwave radiative fluxes are calculated following the method of Parkinson and Washington (1979) based on the SLP and T2m fields and spatially uniform, seasonally varying, climatological cloud fractions. Model input also includes river runoff and precipitation detailed in Hibler and Bryan (1987) and Zhang et al. (1998).

The ice concentration is assimilated from a dataset originally created by Chapman and Walsh (1993). The global ice concentration dataset (Gice) is obtained from the British Atmospheric Data Centre [a more recent version is the Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) dataset (Rayner et al. 1996) from the same center]. It consists of monthly averaged ice concentration on a 1° grid. In the satellite era it is based largely on various satellite measurements and in the presatellite era on ship reports and climatology. The ice concentration is assimilated in a method that emphasizes ice extent by heavily weighting the observations only when there is a large discrepancy between the model and the observed concentration (Lindsay and Zhang 2006). Each day the model estimate *C*_{mod} is nudged to a revised estimate *Ĉ*_{mod} with the relationship

and the gain (or weighting) function is

where *C*_{obs} is the observed concentration, *R*^{2} is the error variance of the observations, and the exponent *α* = 6. We use a fixed value of *R* = 0.05 that is consistent with the estimated errors of the Gice dataset. The observations are heavily weighted if |*C*_{obs} − *C*_{mod}| is greater than 0.5.

The ice velocity measurements are assimilated with an optimal interpolation scheme (Zhang et al. 2003). We use velocity measurements from both buoy- and Special Sensor Microwave Imager (SSM/I)-derived ice displacement measurements. The buoy velocities were obtained from the International Arctic Buoy Program (IABP), and SSM/I 85-Ghz ice displacement measurements were provided by the Jet Propulsion Laboratory Polar Remote Sensing Group. The buoy velocities are 24-h averages and the SSM/I velocities are based on 2-day displacements. The passive microwave displacement estimates are based on a maximum correlation method applied to sequential images of the ice cover. While the SSM/I daily velocity estimates have a substantially larger error standard deviation than the buoys (0.057 versus 0.007 m s^{−1}) their large number and excellent spatial coverage make them a valuable addition to our analysis.

Model spinup consists of an integration of 18 yr using 1948 forcing repeatedly. This allows the model thickness distribution to reach an approximate steady state from an initial uniform 3-m-thick ice cover. Sensitivity studies indicate that influence of the initial ice thickness on the thickness distribution becomes small after about 5 yr of spinup integration. To allow for a spinup time and to avoid the first 5 yr of the study period when the ice concentration data is largely based on climatological values, our analysis is confined to the 51-yr period, 1953–2003.

We use the 51 yr from our assimilation study (Lindsay and Zhang 2006) that includes assimilation of both ice concentration and ice velocity. Ice velocity measurements are sparse before 1979, so velocity assimilation begins with that year. The assimilation study compared the model ice draft with that measured by submarines in the period 1987–97 and found that the model-only simulation has a bias of −0.11 m and a correlation of *R* = 0.51 (*N* = 835). With assimilation of ice concentration the simulated mean thickness has a bias of −0.30 m and the correlation with the observed ice draft is *R* = 0.63. When ice velocity is also assimilated, the bias is reduced to –0.02 m and the correlation increases to *R* = 0.70. A strong spatial pattern in the bias of the model ice draft in comparison with the submarine ice draft was observed; the model shows ice averaging about 1 m too thick on the Pacific Ocean side of the basin and 1 m too thin on the Atlantic Ocean side. The spatial pattern of the model bias is puzzling because it persists even when ice velocity measurements are assimilated and hence the mean advective patterns are well estimated. This suggests that there may be some large-scale error in the model physics or in the forcing. The assumption of spatially uniform cloud cover may be an error and accounting for spatial variability in the clouds may reduce the spatial bias.

## 3. Variability and trends in ice thickness

### a. Mean and variability of the ice thickness for the Arctic Ocean

The annual mean ice thickness over the entire Arctic Ocean is shown in Fig. 2. The model-only simulation shows a mean thickness slightly greater than that seen with assimilation of ice concentration and ice velocity. The difference between the two is reduced to near zero shortly after the assimilation of ice velocity begins in 1979. The increase with assimilation of ice velocity may reflect increased model ice production due to increased ice deformation created by the velocity data assimilation. As mentioned above, the assimilated velocity data improves the comparisons with ice draft measurements. The thickness is the greatest in 1987 and there is a secondary maximum in the mid-1960s. There is little change in the mean thickness in the year after the 1987 maximum, but in the following year, 1989, the ice begins the rapid thinning trend that persists to the present. The decrease in the mean thickness over the Arctic Ocean from 1987 to 2003 is −1.33 m or 40%, and the ice in 2003 is at a minimum mean thickness for the 51-yr period. The maxima in 1966 and 1987 are similar to what Rothrock et al. (2003) found when looking at the published mean ice thickness from eight different sea ice models. Five of the eight models show a local maximum in 1966 and all eight show a local maximum in 1987.

The mean ice thickness and the standard deviation of the annual mean thickness for the 51-yr simulation are shown in Fig. 3. The mean ice thickness is greatest along the northern Canadian coast and is fairly uniform at near 3 m over much of the central part of the basin. The interannual variability is least in the central part of the basin, about 0.4 m. It is large both near the Canadian coast and in the East Siberian Sea with maxima near 1.0 m. The mean pattern represented here may be in error if the mean differences between the observed and modeled ice draft shown in Lindsay and Zhang (2006) can be extrapolated to the entire time series of ice thickness maps. The mean pattern would then show about 1-m-thinner ice on the Pacific side and 1-m-thicker ice on the Atlantic side of the basin. It is not known if such an extrapolation is appropriate, but such a bias could have important implications. A simple bias would not change the EOF analysis which focuses on variability, but a change in the mean ice thickness could change ice ridging and ice production rates, which could change the patterns of variability seen in the EOFs. However just how the variability would change is not known.

The trend in the annual mean ice thickness over the period 1953–2003 shows that the thinning rate is greatest near the New Siberian Islands where it amounts to about –0.037 m yr^{−1}, or 1.9 m over the 51 yr (Fig. 4). The thinning extends in a broad band toward Fram Strait. There is also thickening during this interval in a narrow band along the Canadian Islands. An analysis of the significance of the trend must include an allowance for the reduction in the degrees of freedom due to the temporal autocorrelation of the ice thickness. The autocorrelation (after removing the trend) indicates that on average there are 19 degrees of freedom in the 51-yr time series, although the value ranges from 6 to 49 depending on location. Much of the area with a positive trend is not significant, while more than half of the area with a negative trend is significant at the 90% confidence level. There is no significant trend in the basinwide mean ice thickness over this period.

In parts of the basin a substantial portion of the interannual variance is accounted for by the trend, in particular the region from the Laptev Sea to the North Pole, where about one-half the interannual variance is accounted for by the trend. In much of the rest of the basin the fraction is much lower, less than 10%. Looking at Fig. 2, clearly the magnitude of a trend in the mean ice thickness is highly dependent on the interval chosen for analysis. The trend in the ice thickness in the last 16 yr, from 1988 to 2003, and the kinematic and thermodynamic processes that contribute to it are considered in detail in Lindsay and Zhang (2005). This period is a time of greatly accelerated thinning of the ice; the trend of the basinwide annual mean thickness is 0.071 m yr^{−1}, significant at the 95% confidence level.

### b. Empirical orthogonal functions

The major modes of spatial variability of the simulated annual mean ice thickness are illustrated with the first three EOFs (Fig. 5). The spatial domain used to determine the EOF patterns included only the Arctic Ocean. The maps show the correlations of the annual mean ice thickness for all ice-covered locations with each of the first three principal component time series. The first mode, representing 30% of the Arctic Ocean variance, shows much of the basin increasing and decreasing in ice thickness in a coordinated fashion. The maximum of the pattern is near the center of the basin, in part because points near the edges have fewer neighbors with which to build high covariability. This bull's-eye pattern with diminished weighting near the edges of the domain is ubiquitous for fields in which a major mode of variability has a spatial scale commensurate with the size of the domain. This pattern likely reflects the basinwide interannual variability in the thermal forcing. The averages of submarine-measured ice draft for entire cruises shown by Rothrock et al. (2003) differ from year to year in a way that is consistent with the large-scale pattern of interannual variability as seen in the first EOF.

The second and third EOF patterns, explaining 18% and 15% of the variance, respectively, together represent modes of lateral variability, meaning points on opposite sides of the basin have an element of anticorrelation. The second EOF shows a mode with opposite signs centered near the Beaufort and Laptev Seas. It is similar in shape and location to the east–west Arctic anomaly pattern (EWAAP) identified by Zhang et al. (2000). The third pattern has centers in the East Siberian Sea and off Svalbard in an orientation orthogonal to the second EOF.

The interpretation of the EOFs is made more difficult because the time series of ice thickness fields is far from stationary (Fig. 2) and in recent years the nature of the thickness field is changing rapidly. The EOFs for the most recent 16 yr, 1988–2003, when the ice is thinning rapidly, are shown in Fig. 6. The first EOF, explaining 41% of the variance, is quite different from that of the 51-yr period. It is more broadly situated and the largest loading is in a band stretching from the Chukchi Sea to north of Canada to Fram Strait. This band is roughly coincident with the region of the greatest thinning rates for this same period seen in Lindsay and Zhang (2005). This is consistent with the fact that about 80% of the ice thickness variance during this period is accounted for by the trend as opposed to just 11% over the 51-yr period. The second EOF is still similar to the EWAAP, though the pattern is more accentuated, and it accounts for a similar amount of the variance as in the 51-yr period. The third EOF is very different from that seen in the longer period and accounts for less of the total variance.

The time series of the first three principal components (PCs) from the 51-yr analysis are shown in Fig. 7. The first component resembles the time series for the mean ice thickness (Fig. 2), reflecting the fact that it accounts for a substantial portion of the variance in the ice thickness fields. The two prominent maxima in the ice thickness correspond to maxima in the first PC and moderate values of the other two. In recent years it is in an extreme negative mode; much of the recent observed thinning of the basinwide mean ice thickness is accounted for by this mode. The second EOF was in a strong positive phase in the early 1990s, shortly after the 1987 maximum, and has since returned to near-zero values while the third had a strong positive phase shortly before the 1987 maximum. The autocorrelations (not shown) indicate that the dominant time scale of the variability of the PCs is greatest for the first and least for the third.

The EOFs and PCs are statistical constructs, but some physical interpretation can be obtained by seeing how well these components correlate with the dominant mode of atmospheric variability for the Northern Hemisphere. The winter (October–January) mean Arctic Oscillation index (AO; Thompson and Wallace 1998) is also shown in Fig. 7 as well as the lagged correlations of each PC with the AO. Only the second and third PCs show a significant correlation, greatest at a half-year lag, which is for the ice thickness in the calendar year following each winter AO value. Winter is the period when the AO has the most variability and the greatest hemispheric impact (Thompson and Wallace 1998). The October–January averaging period was chosen for the AO so as to maximize the correlation of the second PC, since the other two PCs show little correlation for any choice of the winter-averaging period. The maximum correlation is, however, not sensitive to the exact period selected for the winter averaging. A positive AO index implies enhanced offshore flow from the Asian sector, in particular near the Laptev Sea (Rigor et al. 2002), which diminishes the ice thickness on the Asian side of the basin. There is also less flow south into the Beaufort Sea from the thick-ice region to the northeast leading to thinner ice in the Beaufort. These patterns correspond to EOF-2 (thinner ice near the Laptev, thicker in the Beaufort) and gives rise to the positive correlation of the winter AO with the second PC. A correspondence of an enhanced EWAAP and the NAO was also noted by Zhang et al. (2000), which is consistent with the correlation found here between the second PC and the AO. A high AO winter is also associated with less ice moving into the Chukchi and East Siberian Seas and more rapid flow into Fram Strait (Rigor et al. 2002). This leads to thinner ice in the Chukchi and East Siberian Seas and thicker ice in Fram Strait, corresponding to the negative of EOF-3 and giving rise to the negative correlation of the third PC with the winter AO.

About one-half of the covariability between the second PC and the previous winter AO is contributed by the four or five years following the ice maximum in 1987. The fact that much of the covariability is contributed by one small interval of years suggests that a predictive scheme for this mode based on the AO is not likely to be robust.

## 4. Monitoring the basinwide mean ice thickness and the thickness spatial distribution

The model output may also be used as a simulated observing system to guide where to best locate fixed moorings to monitor the changing mean ice thickness of the basin as well as the changing spatial variability of the thickness. Beginning in the spring of 2001 a mooring has been deployed at the North Pole as part of the North Pole Environmental Observatory (Morison et al. 2002). A second mooring was deployed by the National Atmospheric and Oceanic Administration (NOAA) in the Chukchi Sea in the summer of 2003. How representative of the basinwide mean thickness are these locations? Figure 8a shows the squared correlation of the annual mean ice thickness at each model node with the basinwide mean value as seen in Fig. 2. The maximum variance explained (65%) is at a point near the Chukchi and East Siberian Seas; however, the North Pole is not a bad place for a first measurement, where the annual mean ice thickness explains 53% of the variance of the basinwide mean. Assuming a measurement is made at the Pole, where might a second point be placed to optimally predict the basinwide mean thickness? Figure 8b shows the additional variance explained by a second point given an observation at the Pole. The two points explain 90% of the interannual variance of the basinwide mean ice thickness. The location of the second maximum is not far from the location of the new Chukchi Sea mooring where the fraction of the variance explained is 0.86 if the measurement at the North Pole is included.

This analysis may help to determine the locations to best monitor the mean ice thickness of the basin but does not determine how to best monitor the temporal changes in the spatial patterns of the ice thickness. One way to find the locations that best represent both the temporal and spatial variability is to determine the points that maximize the explained variance for the most important PCs. Each PC represents the annual weighting (i.e., the temporal variability) of the spatial pattern of the corresponding EOF (i.e., the spatial variability). The PCs and EOFs are an efficient representation of the both spatial and temporal variability of the ice thickness estimated for all 4400 model grid cells that represent the Arctic Ocean. Here we use only the first 10 PCs, which embody 86% of the variance.

First, each of the PCs is normalized so that its variance *σ*_{i}^{2} is the same as the fraction of the total variance explained by the corresponding EOF, thus the first PC has a variance of 0.30, the second a variance of 0.18, etc. The mean annual thickness at each grid location *h*(*x, y, t*) is then correlated with the each of the first 10 PCs, *ϕ*_{i}(*t*), *i* = 1–10, where (*x, y*) is the location and *t* is the year. The correlation coefficient at each location is *R _{i}*(

*x, y*). The fraction of the total variance explained by each correlation is

*R*(

_{i}*x, y*)

^{2}

*σ*

_{i}

^{2}and the total variance explained by all 10 PCs is

Figure 9a shows the spatial pattern of the total variance explained by one point. The maximum (28%) is in the Siberian sector of the central basin at a location that shows strong loading on the first EOF (27%) and lesser loading on both of the next two (4% each). Figure 5 shows the spatial patterns associated with these three PCs.

Rather than find the optimum location for two points, which would require testing all pairs of grid points, we recognize the fact that there is an ongoing ice draft measurement program at the North Pole and investigate where a second point might be located assuming a measurement there. A multiple linear regression for each of the *ϕ*_{i}(*t*) with the ice thickness at the Pole *h*(0, 0, *t*) and with each of the other ice-covered grid points *h*(*x, y, t*) yields *R*_{2,i}(*x, y*)^{2}*σ*_{i}^{2} as the fraction of the total variance explained by each of these two-point correlations. Figure 9b shows the variance explained by the ice thickness measured at two locations one of which is the North Pole. The maximum total variance explained (44%) is located in the Beaufort Sea, but many locations in the Beaufort and Chukchi Seas do well. This location primarily measures the second PC because the North Pole location is well situated to monitor the first.

Last, we recognize that there is a new measurement site in the Chukchi Sea established by NOAA and we ask where a third location might be placed if measurements are available from both the North Pole and from the new mooring, CH01. Here a three-variable multiple regression is performed for each grid location. Figure 9c shows the variance explained by each location. Themaximum is in the East Siberian Sea and, with the two established moorings, it accounts for 57% of the total variance. This location primarily monitors the second PC (6% additional variance explained) and the third PC (7%).

## 5. Comments and conclusions

Our model simulations of the mean ice thickness in the Arctic Ocean show two prominent maxima: one in 1966 and one in 1987. The spatial pattern of the mean thickness exhibits the thickest ice in a thin band along the Canadian coast. The standard deviation of the annual mean thickness is least in the center of the basin and greatest in the East Siberian and Beaufort Seas. The mean thickness has shown a thinning trend in the 1953–2003 period, primarily on the Siberian coast and along the transpolar drift stream to Fram Strait. Since 1987 the simulated ice thinned considerably and there has been increasing amounts of summer open water. The details of the physical processes that lead to the dramatic thinning of the ice, particularly in the 16-yr period since 1987, are found in Lindsay and Zhang (2005).

The spatial patterns of the ice thickness variability were explored with EOFs. The first EOF pattern accounts for 30% of the interannual variability of the mean ice thickness in the Arctic Ocean. The pattern is nearly basin wide and probably is largely thermodynamically driven because of its broad shape. It is not strongly related to the winter AO index. The second EOF (18%) represents an oscillation between the Beaufort Sea–Canadian coast and the East Siberian– Laptev Seas. It is very similar to the east–west Arctic anomaly pattern and is well correlated with the winter AO index. It likely has a strong dynamic dependence because of this correlation. The third EOF (15%) represents an oscillation between the Chukchi–East Siberian Seas and the Svalbard region and is also correlated with the winter AO index, although not as strongly.

The model results were used to test potential ice thickness observation locations for their representation of the basinwide mean thickness. The single most effective location was found to be in the East Siberian Sea. Because a measurement program is already established at the North Pole, we determined a second location that would most aid in measuring the basinwide mean thickness and found it was in the Chukchi Sea. Between them they monitor 90% of the variance of the basinwide mean thickness. Because the two points explained so much of the variance, a third point was not sought.

The model results were also used to determine the optimal locations to monitor both the spatial and the temporal variability of the annual mean ice thickness. We introduce a new method to do this that uses the first 10 principal components to represent the variability and determines the locations that best explain it in multiple linear regression models. We test for the best single location, the best additional location given a measurement at the North Pole, and the best additional location given measurements at both the North Pole and in the Chukchi Sea. With the three stations, 57% of the spatial and temporal variability of the annual mean ice thickness in the Arctic Ocean is monitored.

While according to this analysis the optimal new location to monitor the basinwide spatial variability of the annual mean thickness is in the East Siberian Sea, other factors might also be important for selecting ice thickness monitoring locations. Among these is the monitoring of locations that show large differences in the simulations from different ice models so as to provide information about which models do well. Also, there is also little information about the mean and variability of the ice thickness in the region north of the Canadian Archipelago, a region that may be particularly important if it is one of the last bastions of multiyear ice in an ice-diminished Arctic Ocean. Of course, logistics requirements will always play an important role in selecting monitoring sites.

While we have shown how EOF analysis can be used to determine where best to monitor both the temporal and the spatial variability of the ice cover, there are two important caveats. First, the method is based on model simulations that may be flawed and second, perhaps more important, it is based on an assumption of a stationary time series: the patterns of variability in the past will be representative of those in the future. As we saw with the EOFs, the patterns of variability are very different in the most recent 16 yr when compared with the whole 51-yr period. There is abundant evidence that the ice cover in the Arctic is changing rapidly from one of predominantly multiyear ice to one with much more extensive summer open water and winter first-year ice. This dramatic regime change makes analyses of the historic ice cover difficult to interpret when planning future monitoring efforts.

## Acknowledgments

The submarine ice draft data were obtained from the National Snow and Ice Data Center, the NCEP–NCAR reanalysis data were from the National Center for Atmospheric Research, the Gice and HadISST ice concentration data were from the British Atmospheric Data Centre, the buoy data were from the International Arctic Buoy Program, and the SSM/I ice motion data were from the Polar Remote Sensing Group at the Jet Propulsion Laboratory. We gratefully acknowledge helpful discussions with M. Steele and A. Schweiger and the helpful comments of the reviewers. This work was supported by the NASA Cryospheric Sciences, the Radiation Sciences, and the Global Modeling and Analysis Programs, the NOAA Arctic Research Program, and the NSF Office of Polar Programs. This publication is partially funded by the Joint Institute for the Study of the Atmosphere and Ocean (JISAO) under NOAA Cooperative Agreement NA17RJ1232.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**.**

**.**

**.**

**,**

**,**

**,**

**,**

**,**

**.**

## Footnotes

*Corresponding author address:* Dr. R. Lindsay, Polar Science Center, University of Washington, 1013 N.E. 40th St., Seattle, WA 98105. Email: lindsay@apl.washington.edu

* Joint Institute for the Study of the Atmosphere and Ocean Contribution Number 1112.